基于频率抽样法的FIR数字低通滤波器的设计

基于频率抽样法的FIR数字低通滤波器的设计1设计目的熟悉频率采样法的理论及其应用；掌握频率采样法设计FIR数字滤波器的方法。了解FIR数字滤波器的频率特性和相位特性，观察过渡带取样点对滤波器幅频特性的影响。掌握用频率采样法设计线性相位FIR低通数字滤波器的方法，并掌握该方法的matlab编程和仿真。2FIR数字滤波器设计的原理2.1频率抽样设计法FIR低通滤波器的设计一般方法有两种，即频率抽样法和窗函数法，频率抽样法设计不同于窗函数法，窗函数是从时域出发，把理想的()dhn用一定形状得窗函数截取成有限长的()hn，以此()hn来近似理想的()dhn，这样得到的频率响应()jwHe逼近于所要求的理想的频率响应()jwdHe。频率抽样法则是从频域出发，把给定的理想频率响应()jwdHe加以等间隔抽样，即2()|()jwddwkNHeHk然后以此()dHk作为实际FIR数字滤波器的频率特性的抽样值()Hk，即令2()()()|0,1,,1jwddwkNHkHkHekN，知道()Hk后，由DFT定义，可以用频域的这N个抽样值()Hk来唯一确定有限长序列()hn，而由()Xz的内插公式知道，利用这N个频域抽样值()Hk同样可求得FIR滤波器的系统函数()Hz及频率响应()jwHe。这个()Hz或()jwHe将逼近()dHz或()jwdHe，()Hz和()jwHe的内插公式为1101()()1NNkkNzHkHzNWz（2.1）102()()()NjwkHeHkwkN（2.2）其中()w是内插函数1()2sin()12()sin()2NjwwNwewN（2.3）将式（2.3）代入（2.2）式，化简后可得11()20sin()12()()sin()2kNNjjwjwNkwNHeeHkewkNN（2.4）即11(1)()20sin[()]12()()sin()2kNNjNjwjwNkwkNNHeeHkewkNN（2.5）从内插公式（2.2）看到，在各频率抽样点上，滤波器的实际频率响应是严格地和理想频率响应数值相等，即22()()()()jkjkNNddHeHkHkHe。但是在抽样点之间的频率响应则是由各抽样点的加权内插函数的延伸叠加而形成的，因而有一定的逼近误差，误差大小取决于理想频率响应曲线形状，理想频率响应特性变化越平缓，则内插值越接近理想值，逼近误差越小，如图2.1梯形理想频率特性所示。反之，如果抽样点之间的理想频率特性变化越陡，则内插值与理想值之误差就越大，因而在理想频率特性的不连续点附近，就会产生肩峰和波纹，如图2.2矩形理想频率特性所示。图2.1梯形理想频率特性图2.2矩形理想频率特性2.2线性相位的约束对于第二类线性相位FIR滤波器，由于()hn偶对称、N为偶数时,公式是1()2()()NjwjwHeHwe，其中幅度函数()Hw应为奇对称的，()(2)HwHw，如果抽样值2()()jkNHkHe也用幅值kH（纯标量）与相角kQ表示，则为：22()()()KKjkjQjQNkHkHeHkeHeN（2.6）其中kQ必须为：121()(1)2kNQkkNN（2.7）kH必须满足奇对称，即kNkHH。2.3线性相位第一种频率抽样210()()NjnkNnHkhne（2.8）当()hn为实数时，满足**()(())()()NNHkHNkRkHNk，由此得出()()HkHNk，()()QkQNk，也就是说，()Hk的模()Hk以2kN为对称中心呈偶对称，()Hk的相角()Qk以2kN为对称中心呈奇对称。再利用线性相位的条件1()2jwNQew，即可得到（N为偶数）：210,,122()0221()1,,122NNkkNNQkkNNNkkNN（2.9）21()221()()2()0,,12()02()1,,12NjkNNjNkNNHkekNHkkNHNkekN（2.10）12(0)sin2sin2NjjNHHeeN121sinsin22sinsin22NkkkNNHkNNkkNNN（2.11）2.4过渡带抽样的优化设计为了提高逼近质量，使逼近误差更小，也就是减小在通带边缘由于抽样点的陡然变化而引起的起伏振荡,这种起伏振荡使阻带内最小衰减变小。这里是使理想频率响应的不连续点的边缘加上一些过渡的抽样点,从而增加过渡带，减小频率边缘的突变，也就减小了起伏振荡，增大了阻带最小衰减。这些周扬点上的取值不同，效果也就不同，如果精心设计过渡带的抽样值，就有可能使它的有用频带（通带、阻带）的纹波得以减小。从而设计出较好的滤波器。一般过渡带取一、二、三点抽样值即可得到满意的结果。在低通设计中，不加过渡抽样点时，阻带最小衰减为－20dB，一点过渡抽样的最优设计，阻带最小衰减可提高到－40dB到－54dB左右，二点过渡抽样的最优设计可达－60dB到－75dB左右，而加三点过渡抽样的最优设计则可达－80dB到－95dB左右。加过渡抽样点的示意图如图2.3所示。图2.3加过渡抽样点3FIR数字滤波器设计方法3.1设计任务本次课程设计的任务是利用频率抽样法设计一个低通FIR数字低通滤波器，其理想频率特性是矩形的，即（3.1）给定抽样频率为421.510(/sec)srad，通带截止频率的321.610(/sec)prad，阻带起始频率为323.110(/sec)strad，通带波动11dB，阻带衰减250dB。幅度特性如图3.1所示。图3.1()Hj幅度特性曲线3.2参数计算通带的截止频率为3421.610220.21321.510pppsswf阻带的起始频率为3423.110220.41321.510stststsswf理想低通截止频率31()22.3510(/sec)2cpstrad其对应的数字频率3422.3510220.31321.510ccws过渡带带宽为0.2stp，设抽样点数为N，由于23wN，故算得抽样点数N等于30。3.3设计方法根据指标，可画出频率抽样后的()Hk序列，由于()Hk是对称于w的，我们又只对0w即015k的区间感兴趣，故可将2w即1630k的图略去不画。截止频率0.313cw，截止频率抽样点的位置应为：0.3133024.6954，按第一种频率抽样方式来设计，N=30，则10Int42()0Int151522ccNkHkNNk（3.2）[]Int表示取整数部分，将这些值代入式(2-11)可得：515130sin30sin30sin230230230sin30sin30sin2230230jjkkkHeekk（3.3）按此式计算1020log()jwHe的结果如图3.2所示。由图看出，过渡带宽为230，而最小阻带衰减约为－20dB。这一衰减在大多数情况下是不令人满意的，也不符合这次课程设计的要求.图3.2幅度频率特性曲线为了改善频率特性，以满足指标要求，可在通带和阻带交界处安排一个或者几个不等于1的抽样值。在本次课设中用优化算法使K＝5处值为0.5886，即(5)0.5886H，和K＝6处值为0.1065，即(6)0.1065H。则得到图3.3所示的结果，过渡带带宽为230.230，最小阻带衰减约为－60dB左右，如图3.4所示。图3.3增加过渡点后的()Hk图3.4过渡带抽样优化后的幅频特性综上所述，采用第一种频率抽样方式来设计，抽样点数取N＝30，优化抽样过渡带采用在通带和阻带交界处安排两个不等于1的点，即(5)0.5886H和(6)0.1065H，最后设计出来的基于频率抽样法FIR数字低通滤波器满足设计要求。4程序设计4.1两种频率抽样方法幅频特性对比程序closeall;clear;N=30;H=[ones(1,4),zeros(1,22),ones(1,4)];k=0:(N/2-1);k1=(N/2+1):(N-1);k2=0;A=[exp(-j*pi*k*(N-1)/N),exp(-j*pi*k2*(N-1)/N),exp(j*pi*(N-k1)*(N-1)/N)];HK=H.*A;hn=ifft(HK);fs=15000;[c,f3]=freqz(hn,1);f3=f3/pi*fs/2;figure(1)subplot(2,1,1)plot(f3,20*log10(abs(c)));title('频谱特性');xlabel('频率/HZ');ylabel('衰减/dB');grid;N=30;H=[ones(1,4),zeros(1,22),ones(1,4)];H(1,5)=0.5886;H(1,26)=0.5886;H(1,6)=0.1065;H(1,25)=0.1065;k=0:(N/2-1);k1=(N/2+1):(N-1);k2=0;A=[exp(-j*pi*k*(N-1)/N),exp(-j*pi*k2*(N-1)/N),exp(j*pi*(N-k1)*(N-1)/N)];HK=H.*A;hn=ifft(HK);fs=15000;[c,f3]=freqz(hn,1);f3=f3/pi*fs/2;subplot(2,1,2)plot(f3,20*log10(abs(c)));title('频谱特性');xlabel('频率/HZ');ylabel('衰减/dB');grid;4.2FIR数字低通滤波器程序closeall;clear;N=30;H=[ones(1,4),zeros(1,22),ones(1,4)];H(1,5)=0.5886;H(1,26)=0.5886;H(1,6)=0.1065;H(1,25)=0.1065;k=0:(N/2-1);k1=(N/2+1):(N-1);k2=0;A=[exp(-j*pi*k*(N-1)/N),exp(-j*pi*k2*(N-1)/N),exp(j*pi*(N-k1)*(N-1)/N)];HK=H.*A;hn=ifft(HK);fs=15000;[c,f3]=freqz(hn,1);f3=f3/pi*fs/2;figure(1);plot(f3,20*log10(abs(c)));title('频谱特性');xlabel('频率/HZ');ylabel('衰减/dB');grid;figure(2);stem(real(hn),'.');line([0,35],[0,0]);xlabel('n');ylabel('Real(h(n))');t=(0:100)/fs;W=sin(2*pi*t*750)+sin(2*pi*t*3000)+sin(2*pi*t*6500);q=filter(hn,1,W);[a,f1]=freqz(W);f1=f1/pi*fs/2;[b,f2]=freqz(q);f2=f2/pi*fs/2;figure(4);plot(f1,abs(a));title('输入波形频谱图