python - 使用 python 从音频文件中提取 1/3 个八度音阶
问题描述
我是信号处理的初学者,我想在 mp3 或 wav 文件上应用三倍频程带通滤波器(产生大约 30 个新的滤波时间序列)中心频率:39 Hz、50 Hz、63 Hz、79 Hz、99 Hz , 125 赫兹, 157 赫兹, 198 赫兹, 250 赫兹, 315 赫兹, 397 赫兹, 500 赫兹, ...</p>
第一种方式...
在我读完 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)
解决方案
根据ANSI S1.11 中的公式 (5) 和 (6):八度音程、半八度音程和第三八度音程滤波器组规范,对于 1/3 八度音程,每个频带的低频和高频由下式给出:
factor = np.power(G, 1.0/6.0)
lowerCutoffFrequency_Hz=centerFrequency_Hz/factor;
upperCutoffFrequency_Hz=centerFrequency_Hz*factor;
whereG
是 2(根据指定的 base-2 规则np.power(10, 0.3)
设计过滤器时),或(根据指定的 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)
推荐阅读
- python-3.x - 如何在 Django 模型中创建 VIEW 表的数据库?
- java - 如何将图像/文件上传到 Tomcat 服务器
- flutter - Flutter:设备方向未在启动画面中显示
- regex - 在正则表达式中实现非贪婪组
- google-apps-script - 使用谷歌工作表添加日历事件,但如果事件已经存在,则不想创建
- ios - 在 AppStore 中打开 App 的正确方法是什么
- json - 如何从 interface{} 转换为 struct 实例
- css - .nav-item.hover 不会改变导航链接的颜色
- c++ - 访问动态分配列表中的元素
- c# - 在 SendMailAsync 方法中出现异常时,任务被取消异常