*The Article is based on a real project we’ve built recently. A more detailed case study coming soon. *

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 using Python**. It sounds like a nightmare for a programmer without specific knowledge about sound and acoustics, doesn't it? Fortunately, 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.

Any comment you want to share with us? Drop us a line: hello@start-up.house. * *