我正在尝试分析每个1/3倍频程频带频率的幅度，因此我使用了许多带通Butterworth滤波器。但是，它们只能在三阶时以50 Hz的频率工作。我想使用6阶，但是由于某些原因，在1 kHz以下没有得到任何结果。
[fs, x_raw] = wavfile.read('ruido_rosa.wav')
x_max=np.amax(np.abs(x_raw))
x=x_raw/x_max
L=len(x)
# Creates the vector with all frequencies
f_center=np.array([50.12, 63.10, 79.43, 100, 125.89, 158.49, 199.53, 251.19, 316.23, 398.11, 501.19, 630.96, 794.33, 1000, 1258.9, 1584.9, 1995.3, 2511.9, 3162.3, 3981.1, 5011.9, 6309.6, 7943.3, 10000, 12589.3, 15848.9])
f_low=np.array([44.7, 56.2, 70.8, 89.1, 112, 141, 178, 224, 282, 355, 447, 562, 708, 891, 1120, 1410, 1780, 2240, 2820, 3550, 4470, 5620, 7080, 8910, 11200, 14100])
f_high=np.array([56.2, 70.8, 89.1, 112, 141, 178, 224, 282, 355, 447, 562, 708, 891, 1120, 1410, 1780, 2240, 2820, 3550, 4470, 5620, 7080, 8910, 11200, 14100, 17800])
L2=len(f_center)
x_filtered=np.zeros((L,L2))
for n in range (L2):
order=6
nyq = 0.5*fs
low = f_low[n]/nyq
high = f_high[n]/nyq
b,a = butter(order,[low,high],btype='band')
x_filtered[:,n] = lfilter(b,a,x)
x_filtered_squared=np.power(x_filtered,2)
x_filtered_sum=np.sqrt(np.sum(x_filtered_squared,axis=0)/L)
pyplot.figure(2)
pyplot.semilogx(f_center,20*np.log10(np.abs(x_filtered_sum)))
pyplot.xlim((50,16000))
pyplot.xlabel('Frequência (Hz)')
pyplot.ylabel('Amplitude (dB)')
如何使用6阶Butterworth滤波器正确过滤50 Hz带通？
参考方案
IIR滤波器的多项式“ ba”表示很容易受到滤波器系数量化误差的影响，该误差会导致极点移动到单位圆之外，并相应地导致不稳定的滤波器。对于具有窄带宽的滤波器，这可能尤其成问题。
为了更好地了解发生了什么，我们可以将使用scipy.butter(..., output='zpk')获得的预期极点位置与通过计算反馈多项式的根（a系数）获得的有效极点位置进行比较。可以使用以下代码完成此操作：
fs = 44100
order = 6
bands = np.array([np.array([2240, 2820]),
np.array([1650,1.2589*1650]),
np.array([1300,1.2589*1300]),
np.array([1120,1410])])
plt.figure()
index = 1
for band in bands:
low = band[0]/(fs/2)
high = band[1]/(fs/2)
b,a=butter(order,[low,high],btype='band',output='ba')
z,p,k=butter(order,[low,high],btype='band',output='zpk')
r = np.roots(a)
t = np.linspace(-np.pi,np.pi,200)
plt.subplot(2,2,index)
index += 1
plt.title("low={:.0f}, high={:.0f} (BW={:.0f}Hz)".format(band[0], band[1], band[1]-band[0]))
plt.plot(np.cos(t),np.sin(t),'k:')
plt.plot(np.real(p),np.imag(p),'bx',label='poles')
plt.plot(np.real(r),np.imag(r),'rs',fillstyle='none',label='roots')
# Focus on the poles/roots at positive frequencies
x0 = np.cos(0.5*np.pi*(low+high))
y0 = np.sin(0.5*np.pi*(low+high))
delta = 0.05
plt.ylim(y0-delta,y0+delta)
plt.xlim(x0-delta,x0+delta)
plt.grid(True)
plt.legend()
这表明这两个集合是为带宽较大的滤波器并置的，但是位置随着带宽的减小而开始变化，直到误差足够大以至于根被推到单位圆之外为止：
为避免此问题，可以使用级联的二阶节过滤器实现：
sos = butter(order,[low,high],btype='band',output='sos')
x_filtered[:,n] = sosfilt(sos,x)
如下图所示，这应将输出扩展到较低的频率：
我在Angular工作，正在使用Http请求和响应。是否可以在“响应”中发送多个参数。角度文件：this.http.get("api/agent/applicationaware").subscribe((data:any)... python文件：def get(request): ... return Response(seriali…
Python exchangelib在子文件夹中读取邮件 - python我想从Outlook邮箱的子文件夹中读取邮件。Inbox ├──myfolder 我可以使用account.inbox.all()阅读收件箱，但我想阅读myfolder中的邮件我尝试了此页面folder部分中的内容，但无法正确完成https://pypi.python.org/pypi/exchangelib/ 参考方案 您需要首先掌握Folder的myfo…
R'relaimpo'软件包的Python端口 - python我需要计算Lindeman-Merenda-Gold（LMG）分数，以进行回归分析。我发现R语言的relaimpo包下有该文件。不幸的是，我对R没有任何经验。我检查了互联网，但找不到。这个程序包有python端口吗？如果不存在，是否可以通过python使用该包？ python参考方案 最近，我遇到了pingouin库。
Python ThreadPoolExecutor抑制异常 - pythonfrom concurrent.futures import ThreadPoolExecutor, wait, ALL_COMPLETED def div_zero(x): print('In div_zero') return x / 0 with ThreadPoolExecutor(max_workers=4) as execut…
如何用'-'解析字符串到节点js本地脚本？ - python我正在使用本地节点js脚本来处理字符串。我陷入了将'-'字符串解析为本地节点js脚本的问题。render.js：#! /usr/bin/env -S node -r esm let argv = require('yargs') .usage('$0 [string]') .argv; console.log(argv…