电子科技大学信号处理实验4-数字滤波器设计及实现

电子科技大学实验报告一、实验室名称：数字信号处理实验室二、实验项目名称：数字滤波器设计及实现三、实验原理：1.数字滤波器设计步骤：（1）根据给定的滤波器设计要求，得到参数化描述，即通带，阻带截止频率p和s，通带阻带纹波p和s等数据。（2）找一个数字系统函数G(z)，使其频率响应逼近设计要求。（3）择合适的滤波器结构对满足要求的传递函数G(z)进行实现。2.数字滤波器设计中的注意事项：（1）设计要求的参数化：图4-1给出了一个典型的数字低通滤波器的幅频特性说明。理解每个参数的物理含义。（2）滤波器类型选择：在数字滤波器实现中可选择IIR滤波器和FIR滤波器两种。在实现相同幅频特性时，IIR滤波器的阶数会相对FIR滤波器的更低；而在实现中，对相同阶数的两种滤波器来看，对每个采样值所做的乘法数量，IIR约为FIR的两倍；另外，FIR还可以方便地设计成)(jeGP1P1sSPPS通带阻带过渡带图4-1典型的数字LPF幅频特性线性相位滤波器。总的来说，IIR滤波器除不能实现线性相位这一点外，由于阶数的原因，从计算复杂度上较FIR滤波器有很大的优势。根据以上这些区别，结合实际的设计要求，就可以选择一款合适的滤波器。（3）波器设计的方法：由于IIR滤波器和FIR滤波器各自的结构特点，所以它们的设计方法也不一样。在IIR滤波器的设计中，常用的方法是：先根据设计要求寻找一个合适的模拟原型滤波器)(sHa，然后根据一定的准则将此模拟原型滤波器转换为数字滤波器)(zG，即为我们需要设计的数字滤波器。在FIR滤波器设计中，一般使用比较直接的方法：根据设计的要求在时域对理想的冲击响应序列进行加窗逼近，或从频域对需要实现的频率响应特性进行采样逼近然后进行反FFT。（4）波器阶数估计：IIR滤波器的阶数就等于所选的模拟原型滤波器的阶数，所以其阶数确定主要是在模拟原型滤波器设计中进行的。FIR滤波器阶数估计可以根据很多工程中的经验公式，这些公式可以直接从设计的参数要求中估计滤波器阶数。例如，对FIR低通滤波器，已知通带截止频率p，阻带截止频率s，最大通带纹波p和最大最带纹波s，则可以使用下面的公式估计其阶数：2/)(6.1413)(log2010psspN3.数字滤波器的设计方法：（1）IIR滤波器设计方法：(a)冲击响应不变法：A.满足设计要求的模拟原型滤波器)(sHa进行部分分式展开为：B.由于，可以得到：NkkkassAsH1)()0)(Re(maxks)()(nThngaNkTskzeAzGk111)((b)双线性变换法：A.设计要求中给出的边界频率进行预畸处理，然后用得到的频率进行模拟滤波器设计，得到模拟原型滤波器)(sHa。B.用双线性变换法求出数字滤波器：1111|)()(zzzasHzG。（2）FIR滤波器设计方法：窗函数法：A.根据指标的要求选择合适的窗函数)(nw，然后根据此窗计算阶数等参数N。B.写出冲击响应序列的表达式：)()()(nwnhnhNd，其中，)(nhd为理想的冲击响应序列，一般为无限长的，)(nwN为长度为N的窗函数。C.计算所得冲击响应序列)(nh的DTFT，然后验证其是否满足设计要求。4.滤波器的实现结构（1）FIR滤波器：直接型实现结构级联结构并联结构线性相位型结构（2）IIR滤波器：直接型实现结构：I型和II型级联结构并联结构5.用MATLAB进行滤波器设计：1）IIR滤波器设计函数：butter,buttord,chebwin,cheb1ord,cheb2ord,cheby1,cheby2,ellip,ellipord等2）FIR滤波器设计函数：fir1,fir2,remez,remezord,kaiser,kaiserord,hanning,hamming,blackman等3）滤波器设计辅助设计软件在命令窗口中键入“fdatool”即可启动滤波器设计辅助设计软件，界面如图4-2所示。在此界面中填写需要设计的滤波器参数，即可设计出需要的滤波器。还可以通过此工具提供的幅度，相位观察窗口观察设计出来的滤波器的幅度，相位特性等，并可以将设计好的滤波器冲激响应系数导出用于实现。图4-2MATLAB中滤波器辅助设计软件界面四、实验目的：熟练掌握数字滤波器的双线性变换法（IIR）和窗函数法（FIR）两种设计方法；加深对数字滤波器的常用指标、设计过程及实现的理解。五、实验内容：本实验要求学生运用MATLAB编程完成IIR和FIR滤波器设计，并与用fdatool设计的滤波器进行对比；用设计出的滤波器系数对离散时间信号进行滤波，分析并观察输出结果。六、实验器材（设备、元器件）：安装了MATLAB软件的PC机一台。七、实验步骤：（1）给定输入信号：FSK信号（输入的二进制待调信号为随机信号，码元频率为100Hz，两个载频分别为2kHz和4kHz，采样频率为20kHz，）。利用MATLAB编程产生该信号，画出其时域和频域图。（2）利用MATLAB编程设计一个数字低通滤波器，指标要求如下：通带截止频率:2.2pfkHz;阻带截止频率:3.5sfkHz;采样频率20TfkHz；通带峰值起伏：1[]pdB；最小阻带衰减：40[]SdB。分别用MATLAB中的IIR和FIR设计命令进行滤波器设计，得出需要的滤波器系数。（3）用MATLAB滤波器辅助设计软件对上述滤波器进行设计，并将得到的滤波器系数对输入信号进行滤波，观察滤波结果。（4）将得到的滤波器系数在MATLAB中编程进行实现（选择直接型实现结果进行实现），对（1）中的输入信号进行滤波（分别用FIR和IIR滤波器进行），观察滤波结果，画出时域和频域图像。（5）（拓展要求）修改需要设计的滤波器的指标要求，比如：将通带截止频率修改为2kHz,或者将最小阻带衰减改为20[]SdB，这时再重复（3）和（5）的步骤，观察所得到的滤波器效果，并对这一结果进行解释。（6）（拓展要求）利用“fdatool”设计IIR级联和并联实现的滤波器系数，分析有限字长的影响，并与实际的滤波结果进行对比。八、实验数据及结果分析：需要的程序代码和结果：（1）产生FSK信号，画出其时域和频域图。%FSK信号的产生程序T=0.1;formatlongfw_b=100;fw1=2000;fw2=4000;fs=20000;ts=1/fs;tw=1/fw_b;tb=0:tw:T;N=length(tb);t=0:ts:T;random=rand(1,N);x_b=zeros(1,N);forn=1:Nif(random(n)=0.5)x_b(n)=1;elsex_b(n)=0;endendx=zeros(1,length(t));m=1;forn=2:Nfork=m:length(t)if(t(k)tb(n))if(x_b(n))x(k)=cos(2*pi*fw1*t(k));elsex(k)=cos(2*pi*fw2*t(k));endm=k+1;endendendfigure(1);subplot(2,1,1);plot(t,x);title('FSK信号的时域图');xlabel('时间/s');ylabel('y');nn=length(x);r=fft(x)/n;r=fftshift(r);f=linspace(-fs/2,fs/2,nn);subplot(2,1,2);plot(f,abs(r));title('FSK信号的频域图');xlabel('频率/Hz');ylabel('幅值/A');（2）用MATLAB中滤波器设计相关函数编程完成IIR和FIR滤波器设计，得出需要的滤波器系数。%IIR滤波器的设计程序Ffs=20000;Wp=2200*2/Ffs;Ws=3500*2/Ffs;Rp=1;Rs=40;[N,Wn]=buttord(Wp,Ws,Rp,Rs);[num,den]=butter(N,Wn);freqz(num,den)%FIR滤波器的设计rp=1;rs=40;Ffs=20000;Ff=[2200,3500];a=[1,0];dev=[(10^(rp/20)-1)/(10^(rp/20)+1),10^(-rs/20)];[N,fo,ao,w]=firpmord(Ff,a,dev,Ffs);Wn=2200/Ffs*2;h_fir=fir1(N,Wn);freqz(h_fir,1)（3）用fdatool完成IIR和FIR滤波器设计，得出需要的滤波器系数。（4）用直接型实现前面的滤波器系数，分别用FIR和IIR滤波器对FSK信号滤波，观察滤波结果，画出时域和频域图像。%IIR滤波器的实现T=0.1;formatlongfw_b=100;fw1=2000;fw2=4000;fs=20000;ts=1/fs;tw=1/fw_b;tb=0:tw:T;N=length(tb);t=0:ts:T;random=rand(1,N);x_b=zeros(1,N);forn=1:Nif(random(n)=0.5)x_b(n)=1;elsex_b(n)=0;endendx=zeros(1,length(t));m=1;forn=2:Nfork=m:length(t)if(t(k)tb(n))if(x_b(n))x(k)=cos(2*pi*fw1*t(k));elsex(k)=cos(2*pi*fw2*t(k));endm=k+1;endendendfigure(1);subplot(2,2,1);plot(t,x);title('FSK信号的时域图');xlabel('时间/s');ylabel('y');nn=length(x);r=fft(x)/n;r=fftshift(r);f=linspace(-fs/2,fs/2,nn);subplot(2,2,2);plot(f,abs(r));title('FSK信号的频域图');xlabel('频率/Hz');ylabel('幅值/A');Ffs=20000;Wp=2200*2/Ffs;Ws=3500*2/Ffs;Rp=1;Rs=40;[N,Wn]=buttord(Wp,Ws,Rp,Rs);[num,den]=butter(N,Wn);y=filter(num,den,x);subplot(2,2,3);plot(t,y);title('IIR滤波器滤波后的时域图');xlabel('时间/s');ylabel('y');r_f=fft(y)/n;r_f=fftshift(r_f);subplot(2,2,4);plot(f,abs(r_f));title('IIR滤波器滤波后的频域图');xlabel('频率/Hz');ylabel('幅值/A');%FIR滤波器的实现T=0.1;formatlongfw_b=100;fw1=2000;fw2=4000;fs=20000;ts=1/fs;tw=1/fw_b;tb=0:tw:T;N=length(tb);t=0:ts:T;random=rand(1,N);x_b=zeros(1,N);forn=1:Nif(random(n)=0.5)x_b(n)=1;elsex_b(n)=0;endendx=zeros(1,length(t));m=1;forn=2:Nfork=m:length(t)if(t(k)tb(n))if(x_b(n))x(k)=cos(2*pi*fw1*t(k));elsex(k)=cos(2*pi*fw2*t(k));endm=k+1;endendendfigure(1);subplot(2,2,1);plot(t,x);title('FSK信号的时域图');xlabel('时间/s');ylabel('y');nn=length(x);r=fft(x)/n;r=fftshift(r);f=linspace(-fs/2,fs/2,nn);subplot(2,2