1/3 октавы из аудио файла с питоном - PullRequest
Новая
       26

1/3 октавы из аудио файла с питоном

1 голос
/ 27 июня 2019

Я новичок в обработке сигналов, и я хотел бы применить полосовые фильтры третьей октавы к mp3 или wav-файлу (дает около 30 новых отфильтрованных временных рядов) центральные частоты: 39 Гц, 50 Гц, 63 Гц, 79 Гц99 Гц, 125 Гц, 157 Гц, 198 Гц, 250 Гц, 315 Гц, 397 Гц, 500 Гц,…

Первый способ ...

После того, как я прочитал mp3 файл, яполучил стерео сигнал.Затем я разбил 1 аудиофайл на 4096 сэмплов.Затем я разделяю его на левый и правый канал.Теперь у меня есть массивы данных для каждого канала.Затем я применяю быстрое преобразование Фурье к образцу файла.Проблема в том, что мне нужно применить полосовые фильтры третьей октавы к этим выборкам.Мне нужно предложение, как мне поступить, так как я совсем не понимаю библиотеку акустики.

Другой способ ...

Я нашел какой-то сайт, связанный с моими ожиданиями, но он использует полосовой фильтр октав.Я использую код из ответа Михаила https://www.dsprelated.com/thread/7036/octave-bandpass-filter-on-audio-wav-files Поэтому я хотел бы применить этот код к третьей октаве. 

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
import math

sampleRate = 44100.0
nyquistRate = sampleRate/2.0

#center = [39, 50, 63, 79, 99, 125, 157, 198, 250, 315, 397, 500, 630, 
794, 1000, 1260, 1588, 2000, 2520, 3176, 4000, 5040, 6352, 8000, 10080, 
12704, 16000, 20160, 2508, 32000]
centerFrequency_Hz = 480.0;
lowerCutoffFrequency_Hz=centerFrequency_Hz/math.sqrt(2);
upperCutoffFrequenc_Hz=centerFrequency_Hz*math.sqrt(2);

# Determine numerator (b) and denominator (a) coefficients of the digital 
# Infinite Impulse Response (IIR) filter.
b, a = signal.butter( N=4, Wn=np.array([ lowerCutoffFrequency_Hz, 
upperCutoffFrequenc_Hz])/nyquistRate, btype='bandpass', analog=False, 
output='ba');

# Compute frequency response of the filter.
w, h = signal.freqz(b, a)

fig = plt.figure()
plt.title('Digital filter frequency response')
ax1 = fig.add_subplot(111)

plt.plot(w, 20 * np.log10(abs(h)), 'b')
plt.ylabel('Amplitude [dB]', color='b')
plt.xlabel('Frequency [rad/sample]')

ax2 = ax1.twinx()
angles = np.unwrap(np.angle(h))
plt.plot(w, angles, 'g')
plt.ylabel('Angle (radians)', color='g')
plt.grid()
plt.axis('tight')
plt.show()

fs, speech = wavfile.read(filename='segmented/atb30.wav');
speech = speech[:, 0]
fig=plt.figure()
plt.title('Speech Signal')
plt.plot(speech)

filteredSpeech=signal.filtfilt(b, a, speech)
fig=plt.figure()
plt.title('480 Hz Octave-band Filtered Speech')
plt.plot(filteredSpeech)

1 Ответ

0 голосов
/ 29 июня 2019

В соответствии с уравнениями (5) и (6) из ANSI S1.11: Спецификация для наборов полосовых фильтров октавы, полуоктавы и третьей октавы , для 1/3 октавы - нижняя и верхняяЧастоты каждой полосы определяются следующим образом: 

factor = np.power(G, 1.0/6.0)
lowerCutoffFrequency_Hz=centerFrequency_Hz/factor;
upperCutoffFrequency_Hz=centerFrequency_Hz*factor;

Где G равно либо 2 (при разработке фильтров в соответствии с указанными правилами base-2), либо np.power(10, 0.3) (при разработке фильтров в соответствии с указаннымиправила базы-10).В вашем случае предоставленные вами центральные частоты близки к значениям, полученным с использованием правил base-2, поэтому вам также следует G = 2 быть последовательным.

Обратите внимание, что для вашего заданного массива центральной частоты есть несколько значенийиз верхних частот будет больше, чем частота Найквиста (половина вашей частоты дискретизации).Соответственно, при попытке создать фильтр с scipy.signal.butter они будут давать недопустимые верхние значения нормализованной частоты.Чтобы избежать этого, вы должны ограничить массив центральной частоты значениями, меньшими, чем ~ 19644 Гц: 

centerFrequency_Hz = np.array([39, 50, 63, 79, 99, 125, 157, 198, 250, 315, 
397, 500, 630, 794, 1000, 1260, 1588, 2000, 2520, 3176, 4000, 5040, 6352, 8000, 10080, 
12704, 16000])

Кроме того, scipy.signal.butter может обрабатывать один набор низких и верхних частот одновременно,поэтому вы должны зациклить массивы нижней и верхней частот для разработки каждого полосового фильтра: 

for lower,upper in zip(lowerCutoffFrequency_Hz, upperCutoffFrequency_Hz):
    # Determine numerator (b) and denominator (a) coefficients of the digital 
    # Infinite Impulse Response (IIR) filter.
    b, a = signal.butter( N=4, Wn=np.array([ lower, 
    upper])/nyquistRate, btype='bandpass', analog=False, 
    output='ba');

    # Compute frequency response of the filter.
    w, h = signal.freqz(b, a)

    plt.plot(w, 20 * np.log10(abs(h)), 'b')

    # Filter signal
    filteredSpeech = signal.filtfilt(b, a, speech)

Это должно дать вам график, подобный следующему для амплитудных откликов:

enter image description here

В этот момент вы можете заметить некоторые признаки нестабильности в нижней полосе.Чтобы избежать этой проблемы, вы можете переключиться на представление sos: 

for lower,upper in zip(lowerCutoffFrequency_Hz, upperCutoffFrequency_Hz):
    # Design filter
    sos = signal.butter( N=4, Wn=np.array([ lower, 
    upper])/nyquistRate, btype='bandpass', analog=False, 
    output='sos');

    # Compute frequency response of the filter.
    w, h = signal.sosfreqz(sos)

    plt.plot(w, 20 * np.log10(abs(h)), 'b')

    # Filter signal
    filteredSpeech = signal.sosfiltfilt(sos, speech)
...