设计多频带的FIR滤波器

有时在信号处理中需要对信号进行多频带的数字滤波，以往的办法是单独设计一个个频带的滤波器，把信号输入到每一个滤波器（并联），再把它们的输出信号叠加在一起。是否能设计一个多频带的滤波器，一次输入/输出把多频带的滤波都处理了。答案是可以做到，下面将介绍用fir2函数和等波纹法设计FIR的多频带滤波器。

用频率采样法设计FIR滤波器
函数：fir2
功能：用频率采样法设计FIR滤波器
调用格式：
B=fir2(n,f,a,window)

说明：其中n是滤波器阶数，f和a是频率和滤波器幅值矢量（类似于kaiserord中的f和a参数），当设计低通滤波器时，f=[fpfs]，a=[1o]；当设计带通滤波器时，f=[fs1 fp1 fp2 fs2]，a=[0 1 0]，对高通和带通滤波器可按低通和带通给出。f的长度是等于a的长度，f必须是归一化频率：0≤f≤1。当然f和a的值不限于低通、高通、带通和带阻滤波器，也可以表示出任意的幅值响应。在fir2函数中用频率采样法设计FIR滤波器后还结合窗函数法优化滤波器响应，所以可以指定窗函数，默认设置为海明窗。除海明窗外，还可用boxcar、hann、bartlett、blackman、kaiser和chebwin等窗函数。B是滤波器系数。

等波纹FIR滤波器的设计
函数：firpm(remez)
功能：用Parks-McClellan方法设计等波纹FIR滤波器
调用格式：
b =firpm(n,f,a)
b = firpm(n,f,a,w)
b = firpm(n,f,a,'ftype')
b = firpm(n,f,a,w,'ftype')

说明：在MATLAB旧版本中曾用remez函数，但在新版本中用firpm替代了，但功能是一样的。在调用格式中，其中n是滤波器阶数，f和a是频率和滤波器幅值矢量，f的长度是2*length（a）-2,f必须是归一化频率：0≤f≤1。w是计权矢量，ftype是滤波器类型。利用本函数可以设计微分器或Hilbert变换器。

各参数和fir2函数中的参数一样。

案例、设待处理信号的采样频率是2000Hz,而需要信号的频率区间为260~340Hz、640~680Hz、860~1000Hz等，这样有3个频带，其中2个是带通滤波器，1个是高通滤波器。要求滤波器的通带波纹为1dB，阻带衰减不低于40dB。设置滤波器参数如下：fcl=220，fc2=260，fc3=340，fc4=380，fc5=520，fc6=560，fc7=640，fc8=680，fc9=820，fc10=860。其中，阻带频率为fcl，fc4，fc5，fc8，fc9，通带频率有fc2，fc3，fc6，fc7，fc10。用fir2函数（频率采样法和窗函数法结合的方法）和等波纹法设计FIR滤波器。

程序如下：

clear all; clc; close all;
 
Fs=2000;                               % 采样频率
Fs2=Fs/2;                              % 奈奎斯特频率
% 滤波器各个频率点值
fc1=220; fc2=260; fc3=340; fc4=380; fc5=520;
fc6=560; fc7=640; fc8=680; fc9=820; fc10=860;
rp=1; as=40;                           % 通带波纹和阻带衰减 
% 归一化各频点值
Fc=[fc1 fc2 fc3 fc4 fc5 fc6 fc7 fc8 fc9 fc10]/Fs2;
% fir2法
f=[0 Fc 1];                            % 构成理想滤波器的频率矢量
a=[0 0 1 1 0 0 1 1 0 0 1 1];           % 构成理想滤波器的幅值矢量
dw=(fc3-fc2)*pi/Fs2;                   % 归一化的过渡带值
N1=ceil(6.6*pi/dw);                    % 计算海明窗时滤波器的阶数
N1=N1+mod(N1,2);                       % 保证滤波器阶数为偶数
b=fir2(N1,f,a);                        % 用fir2函数求出滤波器系数
[db1,mag1,pha1,grd1,w1]=freqz_m(b,1);  % 求出滤波器频域响应
% 等波纹法
A=[0 1 0 1 0 1];                       % 构成幅值矢量
devp=(10^(rp/20)-1)/(10^(rp/20)+1);    % 求出通带的偏差值
devs=10^(-as/20);                      % 求出阻带的偏差值
dev=[devs devp devs devp devs devp];   % 构成滤波器设计中的偏差矢量
[N2,F0,A0,W]=firpmord(Fc,A,dev);       % 用firpmord求出滤波器的阶数
N2=N2+mod(N2,2);                       % 保证滤波器阶数为偶数
df=Fs2/500;                            % 频域分辨率
ns1=ceil(fc1/df)+1;                    % fc1对应的样点值
wlis=1:ns1;                            % 阻带样点区间
Acs=1;                                 % Acs初始值
while Acs>-as                          % 阻带衰减不满足条件将循环
    h=firpm(N2,F0,A0,W);               % 用firpm函数设计滤波器
    [db2, mag2, pha2, grd2,w2]=freqz_m(h,1);  % 计算滤波器频域响应
    Acs=max(db2(wlis));                % 求阻带衰减值
    fprintf('N=%4d   As=%5.2f\n',N2,-Acs); % 显示滤波器阶数和衰减值
    N2=N2+2;                           % 阶数加2,保证为第1类滤波器
end
N2=N2-2;                               % 修正N2值
% 作图
subplot 211; plot(w1/pi*Fs2,db1,'k','linewidth',2)
grid; axis([0 1000 -80 10]);
set(gca,'XTickMode','manual','XTick',[0 220 260 340 380 520 560 640 680 820 860 1000])
set(gca,'YTickMode','manual','YTick',[-40,0])
title('(a)fir2函数设计滤波器幅值响应'); 
xlabel('频率/kHz'); ylabel('幅值/dB')
subplot 212; plot(w2/pi*Fs2,db2,'k','linewidth',2)
grid; axis([0 1000 -80 10]);
set(gca,'XTickMode','manual','XTick',[0 220 260 340 380 520 560 640 680 820 860 1000])
set(gca,'YTickMode','manual','YTick',[-40,0])
title('(b)等波纹法设计滤波器幅值响应'); 
xlabel('频率/kHz'); ylabel('幅值/dB')
set(gcf,'color','w');

运行结果如下：

①用fir2函数设计中用了海明窗函数，以海明窗计算了滤波器的阶数N1，但从图3-13-13（a）中可以看到，过渡带较宽，在设置的阻带频率点上（220Hz、360Hz、520Hz、680Hz、820Hz）并没有满足设计要求，即衰减没有达到40dB。也就是说，fir2+海明窗函数设计中计算滤波器阶数的方法对于多频带不完全适用。

②用等波纹设计FIR多频带滤波器，虽然计算出来的阶数N2也不满足阻带的衰减要求，但经过几次增加阶数的循环，最后还是获得了满意的结果。在本例中fir2函数求得滤波器阶数N1为84，等波纹法求得滤波器阶数N2为88，两者较为接近，但在幅值响应方面显然等波纹法要优于fir2函数法。

参考文献：MATLAB数字信号处理85个实用案例精讲——入门到进阶；宋知用（编著）