Matlab 与 Python 基于窗函数的滤波器设计对比 之 凯瑟窗

Kaiser 窗近似于扁长椭圆形窗，它使主瓣能量与旁瓣能量之比最大。对于特定长度的 Kaiser 窗，参数 β 控制相对旁瓣衰减。对于给定的 β，相对旁瓣衰减相对于窗长度是固定的。随着 β 的增加，相对旁瓣衰减降低，主瓣宽度增加。(以上内容来自: Kaiser 窗)

凯瑟窗的参数估计

Matlab与Python都提供了kaiserord函数, 用于估算特定频率响应条件下的Kaiser窗的阶数(阶数为窗长度-1)和$\beta$参数.

python中的kaiserord特指scipy.signal模块中的kaiserord函数.

Matlab 中的 kaiserord

Matlab中的函数语法形式如下:

1. [n,Wn,beta,ftype] = kaiserord(f,a,dev)
2. [n,Wn,beta,ftype] = kaiserord(f,a,dev,fs)
3. c = kaiserord(f,a,dev,fs,‘cell’)

参数:

• f : 通带边界频率和阻带截止频率构成的向量, 对于低通和高通滤波器, f的长度为2, 对于带通和带阻滤波器, f的长度为4; 形式1中的f为数字归一化频率.
• a: 频带幅度. 低通为[1, 0], 高通为[0, 1], 带通为[0,1,0], 带阻为[1,0,1]
• dev: 最大允许偏差(线性值).
• fs: 采样率

返回值:

• n: 滤波器阶数. n = 窗长度 - 1
• Wn: 数字归一化频率
• bata: $\beta$参数
• ftype: 滤波器类型. 取值: 'low' | 'bandpass' | 'high' | 'stop' | 'DC-0' | 'DC-1'

Python 中的 kaiserord

Python中的函数定义如下:

numtaps, beta = scipy.signal.kaiserord(ripple, width)

参数:

• ripple: 最大允许偏差(dB)
• width: 数字归一化的过渡带宽, 等于通带边界频率与阻带截止频率的差.
  返回值:
• numtaps: 窗的长度.
• beta: $\beta$参数

区别

因为线性相位特性为Type II(当win_len为偶数时)的滤波器无法实现高通滤波器和带阻滤波器(这里不做论证, 感兴趣的自行查找). 因此使用Matlab中的kaiserord 估算高通和带阻滤波器的阶数时, 阶数n必为偶数(win_len为奇数), 但python中的scipy.signal.kaiserord并不知道要实现何种滤波器, 因此numtaps可能为偶数 (参见:高通滤波器设计示例和带阻滤波器设计示例).

基于窗函数的滤波器设计函数

Matlab中使用fir1函数实现基于窗函数的滤波器设计, python中使用scipy.signal.firwin函数来实现. 这两个函数都能够实现Type I(当win_len为奇数时) 和 Type II(当win_len为偶数时)的滤波器设计.

Matlab 中的 fir1

Matlab中的fir1的语法形式如下:

• b = fir1(n,Wn,ftype, window, scaleopt)

参数:

• n: 滤波器阶数
• Wn: 归一化数字频率
• ftype: 滤波器类型. 取值: ‘low’ | ‘bandpass’ | ‘high’ | ‘stop’ | ‘DC-0’ | ‘DC-1’
• window: 窗向量
• scaleopt: 缩放选项. 取值: ‘scale’ (默认) | ‘noscale’
  • ‘scale’: 对系数进行缩放, 使得通带的幅频响应为 0 dB.
  • ‘noscale’: 不对系数进行缩放.

返回值:

b: 滤波器系数

Python 中的 firwin

python中的firwin的语法形式如下:

scipy.signal.firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True, nyq=None, fs=None)

参数:

• numtaps: 窗长度.
• cutoff: 截止频率. 对于低通和高通, 该参数为一个数值;对于带通和带阻, 该参数为一个二元元组.
• width: 过渡带带宽.
• window: 窗的名字, 或者包含窗名和参数的元组, 如: ('kaiser', beta).
• scale: 是否缩放滤波器系数, 使得通带中心频率的幅度响应为0 dB.
• pass_zero: 表示直流信号(频率为0)是否通过. 取值: True 或 False, 也可直接给滤波器的类型: { ‘bandpass’, ‘lowpass’, ‘highpass’, ‘bandstop’}. 对于低通和带阻, 该参数设为: True, 对于高通和带通, 该参数设为False.

返回值:

• h: 滤波器系数

区别

Matlab中的fir1函数同python中的scipy.signal.firwin函数完成相同的功能. 就凯瑟窗而言, 主要区别如下:

• Matlab的fir1函数不需要指定滤波器的类型(低通, 高通, 带通或带阻), 因为已经在kaiserord已经指定了滤波器类型. 而python中的firwin函数需要通过pass_zero参数指定滤波器类型.
• fir1函数使用滤波器的阶数n, 而firwin函数使用窗长度numtaps

滤波器设计举例

以下设计示例中, Matlab 和 Python生成的滤波器系数完全相同.

低通滤波器设计

基于凯瑟窗, 设计一个低通滤波器, 采样频率为1000 Hz, 通带边界频率为 150 Hz, 阻带截止频率为200 Hz, 通带最大允许偏差为0.05, 阻带最大允许偏差为 0.01.

• 使用Matlab 设计此滤波器

figure;
fs = 1000;
fcuts = [150, 200];
mags = [1, 0];
devs = [0.05, 0.01];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(hh,1,1024,fs)
title('低通滤波器')

• 使用Python设计此滤波器

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

fs = 1000
fcuts = np.array([150, 200]) / (fs / 2)
devs = np.array([0.05, 0.01])
ripple = np.min(20 * np.log10(devs))
width = fcuts[1] - fcuts[0]
numtaps, beta = signal.kaiserord(ripple, width)
taps = signal.firwin(numtaps, np.mean(fcuts), window=('kaiser', beta), 
                     scale=False)
w, h = signal.freqz(taps, fs=fs)
plt.plot(w, np.log10(np.abs(h)) * 20)
plt.grid(True)

高通滤波器设计

基于凯瑟窗, 设计一个高通滤波器, 采样频率为1000 Hz, 通带边界频率为 350 Hz, 阻带截止频率为300 Hz, 阻带最大允许偏差为0.05, 通带最大允许偏差为 0.01.

• 使用Matlab 设计此滤波器

figure;
fs = 1000;
fcuts = [300, 350];
mags = [0, 1];
devs = [0.05, 0.01];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(hh,1,1024,fs)
title('高通滤波器')

• 使用Python 设计此滤波器

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

fs = 1000
fcuts = np.array([300, 350]) / (fs / 2)
devs = np.array([0.05, 0.01])
ripple = np.min(20 * np.log10(devs))
width = fcuts[1] - fcuts[0]
numtaps, beta = signal.kaiserord(ripple, width)
if numtaps % 2 == 0:
    numtaps += 1
taps = signal.firwin(numtaps, np.mean(fcuts), window=('kaiser', beta), 
                     scale=False, pass_zero=False)
w, h = signal.freqz(taps, fs=fs)
plt.plot(w, np.log10(np.abs(h)) * 20)
plt.grid(True)

带通滤波器设计

基于凯瑟窗, 设计一个带通滤波器, 采样频率为1000 Hz, 通带边界频率分别为 200 Hz 和 300 Hz, 阻带截止频率分别为150 Hz 和 350 Hz, 通带最大允许偏差为 0.05, 阻带最大允许偏差为 0.01.

• 使用Matlab 设计此滤波器

figure;
fs = 1000;
fcuts = [150, 200, 300, 350];
mags = [0, 1, 0];
devs = [0.01, 0.05, 0.01];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(hh,1,1024,fs)
title('带通滤波器')

• 使用Python 设计此滤波器

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

fs = 1000
fcuts = np.array([150, 200, 300, 350]) / (fs / 2)
devs = np.array([0.01, 0.05, 0.01])
ripple = np.min(20 * np.log10(devs))
width = np.min(fcuts[1::2] - fcuts[0::2])
numtaps, beta = signal.kaiserord(ripple, width)
taps = signal.firwin(numtaps, (fcuts[1] - width / 2, fcuts[2] + width / 2), 
                     window=('kaiser', beta), scale=False, pass_zero=False)
w, h = signal.freqz(taps, fs=fs)
plt.plot(w, np.log10(np.abs(h)) * 20)
plt.grid(True)

带阻滤波器设计

基于凯瑟窗, 设计一个带阻滤波器, 采样频率为1000 Hz, 通带边界频率分别为150 Hz 和 350 Hz, 阻带截止频率分别为 200 Hz 和 300 Hz, 阻带最大允许偏差为 0.05, 通带最大允许偏差为 0.01.

• 使用Matlab 设计此滤波器

figure;
fs = 1000;
fcuts = [150, 200, 300, 350];
mags = [1, 0, 1];
devs = [0.01, 0.05, 0.01];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(hh,1,1024,fs)
title('带阻滤波器')

• 使用Python 设计此滤波器

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

fs = 1000
fcuts = np.array([150, 200, 300, 350]) / (fs / 2)
devs = np.array([0.01, 0.05, 0.01])
ripple = np.min(20 * np.log10(devs))
width = np.min(fcuts[1::2] - fcuts[0::2])
numtaps, beta = signal.kaiserord(ripple, width)
if numtaps % 2 == 0:
    numtaps += 1
taps = signal.firwin(numtaps, (fcuts[1] - width / 2, fcuts[2] + width / 2), 
                     window=('kaiser', beta), scale=False)
w, h = signal.freqz(taps, fs=fs)
plt.plot(w, np.log10(np.abs(h)) * 20)
plt.grid(True)

以上内容为个人的经验总结. 有不当之处, 欢迎留言指正!!!