Week 4 [05-11.11.18] Sound analysis and reverberation time of rooms using Python

Today I would like to share with you my first article which I wrote after finish work on interesting project for company which is produces acoustic furnitures. I hope that this article will be quite interesting for most of you because most of us are programmers.

"Have you ever thought that as a web developer you would have to deal with sound analysis? Neither have I but I faced it during a project for a company which is the leading Polish acoustic furniture producer. The case was a web app with a solution to record sound, analyse it and calculate reverberation time. It sounds like a nightmare for a programmer without specific knowledge about sound and acoustics, doesn't it? Fortunately, we have Python and a lot of Python libraries which come in handy!

So let me walk you through the whole process. Would you like to learn how to analyse sound?

Read a wav file to measure reverberation time

First of all, you have to read a wav file and find out some details from it. In Python it is a simple task because you only need a library called the SciPy, which is quite extensive. But let's focus only on the package scipy.io.wavfile. You can download a sample sound file with strong signal (which is used to measure reverberation time) to compare your results with mine: download (watch out for level of your volume).

from scipy.io import wavfile

sample_rate, data = wavfile.read('file-name.wav')

The method above returns two variables: sample rate of wav file and data (an array containing numbers). The numbers can have different ranges and types depending on the different types of wav formats.

In our case we have a int16 NumPy dtype that is simple to check. Let's see what we can 

read from our file.

from scipy.io import wavfile

sample_rate, data = wavfile.read('STE-010_mono.wav')

amount_of_samples = len(data)

length_of_sound = amount_of_samples / sample_rate

print("Sample rate:", sample_rate)

print("Amount of samples:", amount_of_samples)

print("Length of sound:", round(length_of_sound, 2), "seconds")

print("Data type:", data.dtype)

Sample rate: 44100

Amount of samples: 112640

Length of sound: 2.55 seconds

Data type: int16

Visualisation of a sound: an amplitude of digital signal

What can you do with such data? You can use another interesting library which is widely used for visualising any kind of data and preparing different kinds of charts: matplotlib. What we want to achieve is a chart which shows us a waveform of a digital signal. You can compare your result with a chart I generated using an audio editor called Audacity.

import numpy as np

import matplotlib.pyplot as plt

from scipy.io import wavfile

sample_rate, data = wavfile.read('STE-010_mono.wav')

amount_of_samples = len(data)

# convert sound array to floating point values ranging from -1 to 1

scaled_data = data / (2.**15)

# “return evenly spaced values within a given interval”

time_array = np.arange(0, float(amount_of_samples), 1) / sample_rate

plt.plot(time_array, scaled_data, linewidth=0.3, alpha=0.7, color='#004bc6')

plt.xlabel('Time (s)')

plt.ylabel('Amplitude')

plt.show()

You have your sound data in a range from -2^15 to 2^15-1 so we can simply scale it to values in a range from -1 to 1 by dividing all the values by 2^15. Then you have to prepare a time array whose length will be the same as the array with our sound data. For that case, you should use a NumPy library, which helps a lot with calculations on arrays. The next step is to use a plot function from matplotlib and, voilá, we have a finished chart.

Visualisation of a sound: a spectrogram

To calculate a reverberation time you need a level of dB in a particular time. You can try to write a code step by step, but I used a specgram method from matplotlib which returns me everything I need. Let's focus on how to get a spectrogram and in the next part, we will go to reverberation time. As in the previous example, you can compare your result with a spectrogram from Audacity.

import numpy as np

import matplotlib.pyplot as plt

from scipy.io import wavfile

sample_rate, data = wavfile.read('STE-010_mono.wav')

spectrum, freqs, t, im = plt.specgram(data, Fs=sample_rate, NFFT=1024, cmap=plt.get_cmap('autumn_r'))

cbar = plt.colorbar(im)

plt.xlabel('Time (s)')

plt.ylabel('Frequency (Hz)')

cbar.set_label('Intensity (dB)')

plt.show()

I won't analyse a parameter for specgram function, but if you want to read more about it I encourage you to check the documentation of matplotlib.pyplot.specgram. Let’s see what this function returns.

Source: https://matplotlib.org/api/_as_gen/matplotlib.pyplot.specgram.html

Decibel values for a specific frequency

import numpy as np

import matplotlib.pyplot as plt

from scipy.io import wavfile

sample_rate, data = wavfile.read('STE-010_mono.wav')

spectrum, freqs, t, im = plt.specgram(data, Fs=sample_rate, NFFT=1024, cmap=plt.get_cmap('autumn_r'))

# you can choose a frequency which you want to check

# print(freqs)

index_of_frequency = np.where(freqs == 4005.17578125)[0][0]

# find a sound data for a particular frequency

data_for_frequency = spectrum[index_of_frequency]

# change a digital signal for a values in decibels

data_in_db = 10 * np.log10(data_for_frequency)

plt.plot(t, data_in_db, linewidth=1, alpha=0.7, color='#004bc6')

plt.xlabel('Time (s)')

plt.ylabel('Power (dB)')

plt.show()

Reverberation time of room

We should be conscious that reverberation time is an important issue because it's really inconvenient to stay in a room where it is too high. Now you will be able to check it yourself  following these steps:

1. Find a maximum value of dB's in our array

2. Slice data and time arrays to the point where you found a maximum value

3. Find a value which will be equal to max-5dB

4. Slice the data and time arrays to the point where you found a max-5dB

5. Find a value which will be equal to max-25dB

6. Calculate the time which takes amplitude to drop from max less 5dB to max less 25dB (RT20)

7. Multiply your time by 3 which gives you a RT60

Does it sound a bit tricky and hard to understand? Let's have a look on a code below.

import numpy as np

import matplotlib.pyplot as plt

from scipy.io import wavfile

sample_rate, data = wavfile.read('STE-010_mono.wav')

spectrum, freqs, t, im = plt.specgram(data, Fs=sample_rate, NFFT=1024, cmap=plt.get_cmap('autumn_r'))

# you can choose a frequency which you want to check

# print(freqs)

index_of_frequency = np.where(freqs == 4005.17578125)[0][0]

# find a sound data for a particular frequency

data_for_frequency = spectrum[index_of_frequency]

# change a digital signal for a values in decibels

data_in_db = 10 * np.log10(data_for_frequency)

plt.figure(2)

plt.plot(t, data_in_db, linewidth=1, alpha=0.7, color='#004bc6')

plt.xlabel('Time (s)')

plt.ylabel('Power (dB)')

# find a index of a max value

index_of_max = np.argmax(data_in_db)

value_of_max = data_in_db[index_of_max]

plt.plot(t[index_of_max], data_in_db[index_of_max], 'go')

# slice our array from a max value

sliced_array = data_in_db[index_of_max:]

value_of_max_less_5 = value_of_max - 5

# find a nearest value because it is a big chance that you won't find exactly a value_of_max_less_5 value

def find_nearest_value(array, value):

array = np.asarray(array)

idx = (np.abs(array - value)).argmin()

return array[idx]

value_of_max_less_5 = find_nearest_value(sliced_array, value_of_max_less_5)

index_of_max_less_5 = np.where(data_in_db == value_of_max_less_5)

plt.plot(t[index_of_max_less_5], data_in_db[index_of_max_less_5], 'yo')

# slice our array from a max-5dB

value_of_max_less_25 = value_of_max - 25

value_of_max_less_25 = find_nearest_value(sliced_array, value_of_max_less_25)

index_of_max_less_25 = np.where(data_in_db == value_of_max_less_25)

plt.plot(t[index_of_max_less_25], data_in_db[index_of_max_less_25], 'ro')

plt.show()

rt20 = (t[index_of_max_less_5] - t[index_of_max_less_25])[0]

rt60 = 3 * rt20

print(round(abs(rt60), 2))

green oval: max

yellow oval: max - 5dB

red oval: max - 25dB

Now you have an algorithm which works on a raw signal. To improve our results we should calculate a reverberation time on data, which was normalized in some way. The best choice for that problem is a linear regression.

Linear regression

# linear regression

from scipy import stats

# find a value which is 35dB less than our max

value_of_max_less_35 = value_of_max - 35

value_of_max_less_35 = find_nearest_value(sliced_array, value_of_max_less_35)

index_of_max_less_35 = np.where(data_in_db == value_of_max_less_35)[0][0]

# slice arrays to from max to max-35dB to calculate a linear regression for it

x = t[index_of_max:index_of_max_less_35]

y = data_in_db[index_of_max:index_of_max_less_35]

# you do not have to worry if the gap between min value in y array and max value is a bit more than 35 dB

slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)

plt.plot(x, intercept + slope*x, 'r', label='linear regression')

linregress = intercept + slope * x

An improved algorithm to calculate a reverberation time of room

The following step is changing your algorithm to calculate the same way as before but on a linear regression. Let's look at a code below. You almost finished!

linregress_data = intercept + slope * x

index_of_max = 0

value_of_max_less_20 = linregress_data[0] - 20

value_of_max_less_20 = find_nearest_value(linregress_data, value_of_max_less_20)

index_of_max_less_20 = np.where(linregress_data == value_of_max_less_20)[0][0]

rt20 = (x[index_of_max] - x[index_of_max_less_20])

rt60 = 3 * rt20

print(round(abs(rt60), 2)) # result = 0.98s

Ok, that’s it! You just wrote a script which allows you to check reverberation time in your room. The solution in this article was a bit simplified compared to the solution in the project for our client, but it is enough for checking reverberation time in your own room. Python libraries give us a possibility to deal with unusual use cases which would be more difficult to implement e.g. in Ruby."

1. Have you ever heard about reverberation time?

2. Do you have any experience with programming in python or ruby language? Maybe can you share with us any interesting libraries or plugins?

3. Do you have an idea which programming language could also deal with case as described in the article? Have you ever had such unusual cases in you job?

source: https://start-up.house/blog/articles/reverberation-time-sound-analysis-python