Question

MATLAB中FIR滤波器如何设计

Answer 1

一般滤波的要求主要是通带边界频率、阻带边界频率、通带最大波纹及阻带最小衰减。

而由FIR滤波器的窗函数基本参数，可以知道，最小阻带衰减只由窗形状决定，不受窗宽N的影响；而过渡带的宽度则既与窗形状有关，且随窗宽N的增加而减小。

这样的话，设计一个FIR滤波器，主要是由阻带最小衰减来确定窗形状，再根据过渡带宽的要求来确定窗宽N。有一个窗函数基本参数表，可以对照着选。然后用MATLAB中fir1函数来设计，其语法格式为：b=fir1(N,wn,'ftype',window)。需简单计算N，wn

例题：
设计一个低通数字滤波器，给定抽样频率为fs=5000Hz，通带截止频率wp=500Hz，阻带起始频率ws=800Hz，阻带衰减不小于-50dB。

解答：
由于阻带衰减为50dB，查表，可选海明窗，其阻带最小衰减为53dB，过渡带宽度为6.6π/N。

MATLAB程序如下：
wp=500*2/5000;% 频率归一化
ws=800*2/5000;
wdel=ws-wp;% 过渡带宽
wn=0.5*(wp+ws);% 近似计算截止频率
N=ceil(6.6*pi/wdel);% 根据过渡带宽度求滤波器阶数
window=hamming(N+1);% 海明窗
b=fir1(N,wn,window);% FIR滤波器设计
freqz(b,1,512);% 查看滤波器幅频及相频特性

追问

非常感谢，你写的比较清楚。那这个通带最大波纹怎么选啊？不能像IIR滤波器那样根据给定的通带波纹设计么？

追答

一般来说窗函数选定了后，N增加只会减小过渡带宽，而不会改变肩峰的相对值，这就是Gibbs效应。对于矩形窗情况下，最大相对肩峰总是保持8.95%不变，对于其他窗，第一旁瓣相对主瓣的幅度衰减也是常数。

因此，在窗函数设计FIR滤波器中，通带最大波纹就不是那么重要了。要查看通带最大波纹，可以加一段话：
omega=linspace(0,wp,500);
h=freqz(b,1,omega);% 计算在通带内的幅频
Rp=20*log10(max(abs(h))/min(abs(h))) % 计算通带最大波纹

另外，除了fir1用窗函数设计，还有fir2函数，对应的是频率抽样法设计FIR滤波器，可以查看相应的。

Answer 2

说起来不太方便，基本就是根据指标进行设计。给你举一个例子。
如果用Hamming窗设计满足指标：Wp=0.6*pi;Ws=0.4*pi;Ap=1;As=45的线性相位FIR高通滤波器。
利用I型线性相位滤波器设计FIR高通滤波器。
clear;
%确定滤波器阶数，并使滤波器为I型
Wp=0.6*pi;Ws=0.4*pi;Ap=1;As=45;
N=ceil(7*pi/(Wp-Ws));
N=mod(N+1,2)+N;
M=N-1;
%hamming窗长度N满足条件
w=hamming(N)';
%理想低通截频，计算截频
Wc=(Wp+Ws)/2;
k=0:M;
hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);hd(0.5*M+1)=hd(0.5*M+1)+1;
h=hd.*w;
omega=linspace(0,pi,512);
mag=freqz(h,[1],omega);
plot(omega/pi,20*log10(abs(mag)));grid;%画出增益响应
omega1=linspace(0,Wp,512);
h1=freqz(h,[1],omega1);
omega2=linspace(0,Ws,512);
h2=freqz(h,[1],omega2);
fprintf('Ap=%.4f\n',-20*log10(max(abs(h1))));
fprintf('As=%.4f\n',-20*log10(max(abs(h2))));%计算AS和AP

在举一个例子，是我以前做过的：
用BLACKMAN窗设计满足Wp=0.4*pi;Ws=0.6*pi;Ap=0.5;As=45;的低通线性相位FIR滤波器：
clear;
Wp=0.4*pi;Ws=0.6*pi;Ap=0.5;As=45;
N=ceil(11.4*pi/(Ws-Wp));
N=mod(N+1,2)+N;
M=N-1;
w=blackman(N)';
Wc=(Wp+Ws)/2;
k=0:M;
hd=(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);
h=hd.*w;
omega=linspace(0,pi,512);
mag=freqz(h,[1],omega);
plot(omega/pi,20*log10(abs(mag)));
grid;

希望对你有帮助~欢迎追问~

追问

谢谢你的回答。不过还是有点看不懂。第一段程序里hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi)算的是什么啊？

追答

这是计算IDTFT得到的结果，对于线性相位FIR类型。
由于有幅度函数A=1(|w|<=wc),A=0(其他) 角度-0.5Mw
可以计算IDTFT得到
hd[k]=w/pi*Sa[wc(k-0.5M)]
然后再进行窗截断。
由于这个公式适用于I和II型线性相位系统，所以编程时一般作为固定语句使用。实际上就是计算其频率响应。