python - 来自带有python的音频文件的1/3倍频程

Question

我是信号处理的初学者，我想在 mp3 或 wav 文件上应用 third-octave 带通滤波器(产生大约 30 个新的过滤时间序列)中心频率:39 Hz、50 Hz、63 Hz、79 Hz、99 Hz、125 Hz、157 Hz、198 Hz、250 Hz、315 Hz、397 Hz、500 Hz……

第一种方式...

在我读取 mp3 文件后，我得到了一个立体声信号。然后我将 1 个音频文件分割成 4096 个样本。然后我将它分成左右声道。现在我有每个 channel 的数据数组。接下来，我将快速傅里叶变换应用于示例文件。问题是我需要对这些样本应用三分之一倍频程带通滤波器。我需要关于我应该怎么做的建议，因为我不太了解声学库。

另一种方式...

我找到了一些与我的期望相关的网站，但他使用的是 Octave 带通滤波器。我使用迈克尔在 https://www.dsprelated.com/thread/7036/octave-bandpass-filter-on-audio-wav-files 上的回复中的代码所以我想将此代码应用于第三个 Octave 音阶。

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)

Answer 1

根据 ANSI S1.11: Specification for Octave, Half-Octave, and Third Octave Band Filter Sets 中的等式 (5) 和 (6) ，对于 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)(根据指定的 base-2 规则设计过滤器时)到指定的 base-10 规则)。在您的情况下，您提供的中心频率接近使用 base-2 规则获得的值，因此您还应该 G = 2 保持一致。

请注意，对于给定的中心频率阵列，较高频率的一些值将大于奈奎斯特频率(采样率的一半)。当尝试使用 scipy.signal.butter 设计滤波器时，这些将相应地产生无效的上归一化频率输入。 .为避免这种情况，您应该将中心频率阵列的值限制在小于 ~19644Hz 的值:

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)

这应该为您提供类似于以下幅值响应的图:

此时您可能会注意到最低波段的一些不稳定迹象。为避免此问题，您可以切换到 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)