Backend

Analiza dźwięku i czasu pogłosu w pomieszczeniach z wykorzystaniem Python

4 kwietnia 2019 Tomasz Trzos

Czy kiedykolwiek pomyślałeś, że jako web developer będziesz zajmował się analizą dźwięku? Ja też nie. Jednak, gdy na horyzoncie pojawił się projekt dla polskiego lidera w produkcji mebli akustycznych, stwierdziłem — dobra, spróbuję. Cel — stworzyć aplikację internetową, która pozwoli nagrać i przeanalizować dźwięk, tak aby obliczyć czas pogłosu pomieszczenia.

Tomasz Trzos. Programista backend w Startup Development House. Jego ulubionymi językami jest ruby i python. Pracuje również z technologiami: Ruby on Rails, Docker i Kubernetes oraz rozwija swoje umiejętności z zakresu devops na platformie Google Cloud Platform. Wolny czas poświęca na ściganie się zdalnie sterowanymi samochodami i budowanie dronów.

Innymi słowy, jako klient chce wiedzieć czy wstawić małą, czy dużą ściankę pochłaniającą dźwięk oraz gdzie najlepiej ją postawić. Spore wyzwanie jeśli akurat nie jesteś akustykiem, czyż nie? Na szczęście z pomocą przyszedł Python i sporo bibliotek, dzięki którym udało mi się zmierzyć z tym wyzwaniem.

Chcesz nauczyć się jak analizować dźwięk? W tym artykule przeprowadzę Cię przez cały proces — od odczytania danych z pliku do ich wizualizacji i obliczania czasu pogłosu.

Zobacz też

Backend Co jest najtrudniejsze w byciu developerem? Seniorzy odpowiadają (cz.3)

Spis treści

Odczytanie pliku WAV

Pierwsze co musimy zrobić to odczytanie pliku i uzyskanie danych, z którymi będziemy pracowali. W języku Python nie jest to zbyt trudne. Jedyne czego potrzebujemy to biblioteka “SciPy”, która jest całkiem obszerna, ale skupmy się tylko na pakiecie “scipy.io.wavfile”. Aby porównać wyniki z tymi opisanymi w artykule możesz pobrać plik dźwiękowy z nagranym silnym sygnałem (wykorzystywanym do obliczania czasu pogłosu):
download (uważaj na poziom głośności swoich słuchawek).

from scipy.io import wavfile
sample_rate, data = wavfile.read('file-name.wav')

1 2	from scipy.io import wavfile sample_rate, data = wavfile.read('file-name.wav')

Powyższa metoda zwraca dwie wartości: częstotliwość próbkowania (ang. sample rate) i tablicę danych, która zawiera wartości różnego typu i w różnych przedziałach (w zależności od formatu pliku WAV).

W przypadku naszego pliku jest to typ int16 NumPy dtype. Poniżej zaprezentowałem fragment kodu, który w łatwy sposób pozwoli sprawdzić podstawowe informacje na temat pliku WAV.

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

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

Wizualizacja dźwięku: amplituda sygnału cyfrowego

Co możemy zrobić z danymi wyciągniętymi w poprzednim kroku? Możemy skorzystać z kolejnej ciekawej biblioteki wykorzystywanej do generowania wielu różnych wykresów: matplotlib. To, czego szukamy, to wykres przebiegu sygnału cyfrowego. Możesz porównać swój wynik z wykresem, który wygenerowałem z wykorzystaniem edytora audio: 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()

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()

Otrzymaliśmy dane w zakresie od -2,15 do 2,15-1. W łatwy sposób możemy je przeskalować do wartości z zakresu od -1 do 1. Wystarczy podzielić wszystkie wartości przez 215. Następnie musimy przygotować tablicę, której długość będzie taka sama jak tablica z danymi pochodzącymi z pliku. Elementy tej tablicy będą zawierały dane dotyczące czasu (w sekundach), co pozwoli nam na przedstawienie przebiegu sygnału w czasie. Do przygotowania takiej tablicy wykorzystamy bibliotekę Numpy, która bardzo ułatwia obliczenia na tablicach. Kolejny krok to wykorzystanie funkcji “plot” z biblioteki matplotlib i voilá — mamy gotowy wykres.

Wizualizacja dźwięku: spektrogram

Do obliczenia czasu pogłosu potrzebujemy poziomu decybeli w określonym czasie. Możesz próbować rozwiązać ten problem krok po kroku, ale ja wykorzystałem metodę “specgram” z biblioteki matplotlib, która zwraca wszystko czego potrzebujemy. Skupmy się na tym jak uzyskać spektrogram i w następnym kroku zajmiemy się obliczaniem czasu pogłosu. Tak jak w poprzedniej części — możesz porównać swój spektrogram ze spektrogramem wygenerowanym przez program 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()

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()

Nie chciałem skupiać się na analizie funkcji ‘specgram’, więc jeśli chcesz dowiedzieć więcej na ten temat to zachęcam Cię do sprawdzenia dokumentacji funkcji: matplotlib.pyplot.specgram.

Wartość decybeli dla określonej częstotliwości

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()

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()

Czas pogłosu pomieszczenia

Powinniśmy wiedzieć, że czas pogłosu to całkiem ważna kwestia, ponieważ rozmowa w pomieszczeniu gdzie czas pogłosu jest zbyt wysoki jest bardzo uciążliwa, a dźwięk słyszymy tak jakbyśmy rozmawiali wewnątrz studni. Po przeczytaniu dalszej części artykułu i wykonaniu poniższych kroków będziesz w stanie sprawdzić czas pogłosu pomieszczenia.

1. Znajdź maksymalną wartość decybeli w tablicy

2. “Przytnij” tablice tak, aby przycięta tablica z danymi dźwięku zaczynała się od maksymalnej wartości

3. Znajdź wartość max-5dB

4. “Przytnij” tablice tak, aby przycięta tablica z danymi dźwięku zaczynała się od max-5dB

5. Znajdź wartość max-25dB

6. Oblicz, po którym amplituda spadnie z max-5dB do max-25dB (RT20)

7. Pomnóż swój czas przez 3 i w rezultacie otrzymasz RT60

Zdaję sobie sprawę, że może to być niejasne dlatego z pomocą przyjdzie nam poniższy kod:

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))

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))

zielone kółko: wartość maksymalna,
żółte kółko: max-5dB,
czerwone kółko: max-25dB.

W ten sposób otrzymaliśmy algorytm, który działa na ‘surowym’ sygnale i nie jest odporny na zakłócenia dźwięku. Aby poprawić nasze wyniki powinniśmy obliczyć czas pogłosu na danych, które zostaną w jakiś sposób znormalizowane. Najlepszym wyborem dla tego problemu jest regresja liniowa.

Regresja liniowa

# 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

# 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

Ulepszony algorytm do obliczenia czasu pogłosu pomieszczeń

Następny krok to zmiana naszego algorytmu do obliczania czasu pogłosu, aby wyliczał go w taki sam sposób jak poprzednio, ale korzystając z regresji liniowej. Prawie skończyliśmy, ale spójrzmy jeszcze na poniższy kod:

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

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

W porządku, to wszystko! Właśnie napisałeś skrypt, który pozwala sprawdzić czas pogłosu w Twoim pokoju. Rozwiązanie i kod zaprezentowany w artykule był nieco uproszczony w porównaniu z rozwiązaniem wykorzystanym w projekcie dla naszego klienta, ale wystarcza do sprawdzenia czasu pogłosu we własnym pokoju.

Dziękuję wszystkim za czas poświęcony na artykuł. Tak jak wspomniałem na wstępie — nie jestem specjalistą z dziedziny dźwięku i akustyki oraz na co dzień programuję w ruby, dlatego jeśli ktoś ma jakieś ciekawe spostrzeżenia to zapraszam do kontaktu!

Zobacz też