有限冲激响应滤波器设计

有限冲激响应滤波器设计实验目的：通过实验深入理解有限冲击响应数字滤波器的概念、基本原理和设计的基本方法。掌握信号处理中滤波器设计的原则和需要注意的问题。能够根据不同要求设计不同的滤波器。了解信号处理中几种不同的滤波器设计方法的优缺点。实验原理：FIR原理及特性：有限冲击响应滤波器（FIR）的系统函数可表示为：1010MnMMnHzbkzbbzbz系统的单位冲击响应序列0,1,2,hhhhM即为012,,,,Mbbbb并且当nM时，0hn。此时系统函数可写成0MnnHzhnz其输入输出关系变成卷积的形式0MnYkhnxkn经过分析，可知FIR滤波器最重要的特点便是具有线性相位。其具体可分为两大类：第一类情况：2MhnhMn第二类情况：22MhnhMnFIR设计：从根本上，FIR滤波器的设计就在于寻找一个系统函数jHe取逼近需要的具有理想响应滤波器jdHe。相应的在时域上就是设法用有限冲击响应序列hn去逼近理想响应的单位冲击响应dhn，而2012jjnddnjjnddHehnehnHeed窗函数法：在时域进行设计，实际情况下hn必然是有限长的，即用一个有限长的窗函数n去截取dhn，得到：=dhnhnn要保证所得到的是线性相位FIR滤波器，必须保证所截取的冲击响应序列是偶对成或者奇对称的。窗函数法的关键在于选择合适的窗函数n，其准则在于过渡带宽和阻带衰减大小的问题，即尽可能的增大主瓣，减少旁瓣，使能量尽可能集中在主瓣。窗口主瓣尽可能窄，使过渡带陡峭。频率抽样法：窗函数法是在时域进行设计，但是滤波器的设计指标往往是从频率给出。频率抽样法则是直接从频域出发，对理想频域响应jdHe在02之间做等间隔采样N点，获得频域上的序列dHk。22==0,1,,1jkjNdddkNHkHeHekN再对长度为N的序列dHk做离散傅里叶逆变换，得hn2110,1,,1NjnkNdkhnHkenNNhn就是频率抽样法设计的滤波器的单位冲击响应。其幅度函数Hk和相位函数（序列）k计算公式如下：当N为奇数时：110,1,,21HkHNkNNkkkNNNkkN当N为偶数时：020,1,,1121HkHNkNHNkNkkNNNkkN采样点理想特性越平缓，其内插值就越接近理想值，逼近误差越小。采样点之间变化激烈，内插逼近误差很大，在理想特性的不连续点附件差生明显的肩峰变化。理想滤波器的边界频率出都是阶跃式变化，频率响应起伏震荡大，使阻带衰减程度降低。为此，在突变采样点之间加入一、二个采样过渡点，缓和突变，减少响应的起伏，加大阻带衰减。契比雪夫最佳一致逼近法：为了弥补窗函数法和频率抽样法的通带不够平，阻带衰减不够大，过渡带过宽，频率边缘不能精确指定的问题，产生了契比雪夫最佳一致逼近法。寻找px使得在区间,ab上误差Expxfx较均匀一致，并且使误差的极限值达到最小。在n阶多项式集合中，寻找多项式ˆpx，使其相对于所有其他多项式px对fx的偏差最小，令maxnaxbEpxfxExpxfx则px是fx的最佳一致逼近的充要条件是Ex在,ab上至少存在2n个交错点，122naxxxb，使得1,1,2,,2,1,2,,1iniiExEinExExin所以122,,,nxxx是Ex的极值点，他们构成了一个交错点组。而契比雪夫多项式cosarccos01Cnnxx在0,1上恰好满足此关系。由此便可以利用契比雪夫多项式设计（过程略）01cos,2MjgnNHeannM来逼近理想jdHe。而阶次N可根据经验公式12211log310Nb来进行估计，其中1为通带纹波，2为阻带纹波，2spb。实验内容：参数初始化，设置通带频率为8001000Hz，采样频率3.3sfkHz，采样点数33N，相应代码如下：clc;clearall;fs=3300;N=33;f1=800;f2=1000;Wn=[2*f1/fs,2*f2/fs];窗函数法设计FIR数字带通滤波器，分别采用矩形窗，hamming窗和blackman窗，并画出其幅频特性和相频特性曲线，相应代码如下：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%窗函数法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%b_box=fir1(N,Wn,boxcar(N+1));%矩形窗b_ham=fir1(N,Wn,hamming(N+1));%汉明窗b_black=fir1(N,Wn,blackman(N+1));%布莱克曼窗[H_box,f_box]=freqz(b_box);[H_ham,f_ham]=freqz(b_ham);[H_black,f_black]=freqz(b_black);figure(1);%画幅频特性曲线holdon;plot(f_box,20*log10(abs(H_box)),'r');plot(f_ham,20*log10(abs(H_ham)),'g');plot(f_black,20*log10(abs(H_black)),'b');gridon;figure(2);%画相频特性曲线holdon;plot(f_box,unwrap(angle(H_box)),'r');plot(f_ham,unwrap(angle(H_ham)),'g');plot(f_black,unwrap(angle(H_black)),'b');gridon;频率抽样法设计FIR数字带通滤波器，为了改善阻带衰减，加入一个过渡点=0.5H，分别求出加入过渡点和未加过渡点情况下hamming窗以及矩形窗的幅频特性与相频特性，并画图，响应代码如下：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%频率抽样法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%F=0:1/N:1;M=zeros(N+1,3);M(2*f1/fs*N+1:2*f2/fs*N+1,:)=1;M(2*f1/fs*N+1,2)=0.5;M(2*f2/fs*N+1,2)=0.5;%Hw(2*f1/fs*N+1,3)=2/3;%Hw(2*f1/fs*N,3)=1/3;%Hw(2*f2/fs*N+1,3)=2/3;%Hw(2*f2/fs*N+2,3)=1/3;bf1_box=fir2(N,F,M(:,1)',boxcar(N+1));bf1_ham=fir2(N,F,M(:,1)',hamming(N+1));bf2_box=fir2(N,F,M(:,2)',boxcar(N+1));bf2_ham=fir2(N,F,M(:,2)',hamming(N+1));[H_fbox1,f_fbox1]=freqz(bf1_box);[H_fham1,f_fham1]=freqz(bf1_ham);[H_fbox2,f_fbox2]=freqz(bf2_box);[H_fham2,f_fham2]=freqz(bf2_ham);figure(3);%画幅频特性曲线holdon;plot(f_fbox1,20*log10(abs(H_fbox1)),'r');plot(f_fbox2,20*log10(abs(H_fbox2)),'g');plot(f_fham1,20*log10(abs(H_fham1)),'b');plot(f_fham2,20*log10(abs(H_fham2)),'y');gridon;figure(4);%画相频特性曲线holdon;plot(f_fbox1,unwrap(angle(H_fbox1)),'r');plot(f_fbox2,unwrap(angle(H_fbox2)),'g');plot(f_fham1,unwrap(angle(H_fham1)),'b');plot(f_fham2,unwrap(angle(H_fham2)),'y');gridon;采用契比雪夫最佳一致逼近法设计FIR数字带通滤波器，分别取通带阻带相同权重，以及加大通带权重，加大阻带权重，画出幅频特性和相频特性曲线，相应代码如下：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Chebyshev最佳一致逼近%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fc=[0,14/33,16/33,20/33,22/33,1];mc=[0,0,1,1,0,0];w=[1,1,1;1,10,1;10,1,10];b_c1=firpm(N,fc,mc,w(1,:));b_c2=firpm(N,fc,mc,w(2,:));b_c3=firpm(N,fc,mc,w(3,:));[H_c1,f_c1]=freqz(b_c1);[H_c2,f_c2]=freqz(b_c2);[H_c3,f_c3]=freqz(b_c3);figure(5);%画幅频特性曲线holdon;plot(f_c1,20*log10(abs(H_c1)),'r');plot(f_c2,20*log10(abs(H_c2)),'g');plot(f_c3,20*log10(abs(H_c3)),'b');gridon;figure(6);%画相频特性曲线holdon;plot(f_c1,unwrap(angle(H_c1)),'r');plot(f_c2,unwrap(angle(H_c2)),'g');plot(f_c3,unwrap(angle(H_c3)),'b');gridon;结果分析：窗函数法设计的FIR数字带通滤波器其幅频特性曲线（20logH，以下幅频特性均为此单位）和相频特性曲线分别入图1和图2所示。可以看出从幅频特性上3dB带宽，矩形窗最好，blackman窗最差，即矩形窗过渡带最陡峭；而主瓣能量集中情况blackman窗最好，矩形窗最差，即blackman窗阻带衰减最高。而hanming窗的特性介于两者之间。从相频特性来分析，三种窗函数都具有着几乎一样较好的线性相位，但是blackman窗在零点处产生了的相移。对比矩形窗，hamming窗以及blackman窗可知000,0minffHblackmanHhamg窗，矩形窗、窗由于00所以blackman窗在零点处可能的相移。图1窗函数法幅频特性图2窗函数法相频特性频率抽样法设计的FIR数字带通滤波器幅频特性和相频特性曲线分别如图3和图4所示。不同窗函数之间的比较与窗函数法相同，不在赘述，不同的是hamming窗的通带最大值出衰减不是0，这是由于一共抽取134N个点，而这样正好使得归一化幅值为1的点没有被采样到，带来了一定的衰减。而相较来说加入过渡点之后无论是3dB带宽还是主瓣能量集中度都有较好的性能。相频特性来说，都具有着较好的特性。图3频率抽样法幅频特性图4频率抽样法相频特性契比雪夫最佳一致逼近法，幅频特性曲线和相频特性曲线分别如图5和图6所示。可以看出阻带加权具有较好的阻带衰减，但是通带上却带来幅值过高的问题。而通带加权有着较好的通带特性，但是阻带衰减却不够大。无加权介于两者之间。相频特性都比较好。图5契