在Python中复制matlab过滤器

Question

我正在将 matlab 过滤器脚本重写为 python，并得到截然不同的输出。有谁知道如何使这些功能相似？

进行了一些初步计算以确定等效纹波和其他因素。

Python 脚本：

from scipy import signal
import math as m
import matplotlib.pyplot as plt    

# Variable Values:
total_ripple = -20*m.log10(0.001) #get decibels
fsamp = 100
Nyqu = fsamp/2
trans_width = (0.7 - 0.1)/Nyqu
cut = (0.7+0.1)/2*(Nyqu)

# Script:
numtaps, beta = signal.kaiserord(total_ripple, trans_width)
taps = signal.firwin(numtaps, cutoff=cut, width = trans_width,  fs = fsamp, window=('kaiser',numtaps+1, beta), scale=False, pass_zero=True)
filtered_data = (signal.lfilter(taps, 1, data1, axis = 0)).tolist()

Matlab 脚本：

mags = [1, 0]
devs = [0.001, 0.1]
fcuts = [0.1 0.7]
fsamp = 100
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fsamp);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
filteredData = filter(hh,1,data);

（编辑：添加变量值） 我知道 matlab 和 python 之间的 fcut 会有所不同。 Python 只接受以分贝为单位的单个 fcuts 值作为阻带。简单地绘制 taps & hh 会给出不同的图。最终的图与实际数据也不同。

Answer 1

进行了一些猜测和测试，但我找到了 Kaiser 窗口的等效项（至少从 Scipy 到 MATLAB）。如果您的过渡宽度为 5 的 40 Hz Kaiser 窗口的 Python 代码是：

freq = 40.
trans_width = 5.
sampling_rate = 2000.
nyquist = sampling_rate/2
ntaps40, beta40 = kaiserord(freq, trans_width/nyquist)

MATLAB 等效项是：

[ntaps40,W,beta40] = kaiserord([freq-trans_width/2 freq+trans_width/2],[1 0],2*trans_width/nyquist*[1 1],sampling_rate)

注意：我不知道为什么 2*trans_width 需要 2 倍！我可以告诉你 scipy 和 MATLAB 中的 ntaps 和 beta 是相同的。

Answer 2

之前的答案对我不起作用，因为 scipy 的 kaiserord 将可接受的纹波（以 dB 为单位）作为第一个参数（而不是以 Hz 为单位的频率）。我对其进行了如下调整，它适用于我的情况（具有 200Hz 通带的高通滤波器）：

import numpy as np
from scipy.signal import kaiserord

data = np.array([]) # your data here
fs = 48000
nyquist = fs / 2
passband = 200

stopband_atten = 60 # the default value for stopband attenuation in Matlab
steepness = 0.85 # the default steepness value in Matlab

passband_normalised = passband / nyquist
tw_percentage = -0.98 * steepness + 0.99 # values from Matlab
tw = tw_percentage * passband_normalised
w_stop_normalised = passband_normalised - tw
w_stop = w_stop_normalised * nyquist
trans_width = passband - w_stop

numtaps, beta = kaiserord(stopband_atten, trans_width/(0.5*fs))
mid_freq = passband - (trans_width / 2)
taps = firwin(numtaps, mid_freq, width=trans_width, scale=True, fs=fs, pass_zero='highpass')
delay = floor(numtaps / 2)
temp_data = np.concatenate([data, np.zeros(delay)])
fltrd = lfilter(taps, 1, temp_data)[delay:]

在 Matlab 中给出与以下类似的值：

[filtered_data, fltr] = highpass(data,200,fs);