一、实验目的
1、加深对数字滤波器的基本设计方法的理解。
2、掌握用模拟滤波器原型设计IIR数字滤波器的方法。
3、掌握对窗函数法设计FIR数字滤波器的基本原理的理解。
4、了解MATLAB有关数字滤波器设计的子函数的调用方法。
二、课后习题
注意：我的第三题代码写的有点问题，改的方法治标不治本，等有时间了去问问我们这个课的老师再改。实验的文档在我的上传资源里面。
(1) 分析信号S=3+0.6sin(2pi90t+pi/2)+ sin(2pi50*t-pi/6)的频谱，并用脉冲响应不变法设计低通滤波器，滤除90Hz的正弦波信号。要求画出：
① 信号滤波前的时域图及频谱图（坐标轴是实际幅度及实际频率）
② 滤波器的脉冲响应、滤波器的幅频响应曲线和相频响应曲线。
③ 信号滤波后的时域图及频谱图（坐标轴是实际幅度及实际频率）
代码：
close all clc clear df=1;N=256; Ts=1/N/df;L=N*Ts;t=[0:Ts:L]; %采样时刻 %信号 S=3+0.6*sin(2*pi*90*t+pi/2)+ sin(2*pi*50*t-pi/6)+0.7*sin(2*pi*80*t-pi/6); %显示原始信号 plot(S); title('原始信号');grid on; figure; Y = fft(S,N); %做FFT变换 Ayy = (abs(Y)); %取模 Ayy=Ayy/(N/2); %换算成实际的幅度 Ayy(1)=Ayy(1)/2; F=([1:N]-1)/Ts/N; %换算成实际的频率值 plot(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果 title('幅度-频率曲线图'); grid on; figure; Pyy=[1:N/2]; for i=1:N/2 Pyy(i)=phase(Y(i)); %计算相位 Pyy(i)=Pyy(i)*180/pi; %换算为角度 end plot(F(1:N/2),Pyy(1:N/2)); %显示相位图 title('相位-频率曲线图'); grid on; figure; %%%%%%%%%%%%%%%%%%%%%%%数字滤波器设计6步 %(1)给定的数字滤波器的设计指标 wp=2*pi*50*Ts; %滤波器的通带截止频率 ws=2*pi*80*Ts; %滤波器的阻带截止频率 Rp=1;As=15; %滤波器的通阻带衰减指标 ripple=10^(-Rp/20); %滤波器的通带衰减对应的幅度值 Attn=10^(-As/20); %滤波器的阻带衰减对应的幅度值 %(2)转换为模拟滤波器的技术指标 Fs=2000;T=1/Fs;Omgp=wp*Fs;Omgs=ws*Fs; %(3)确定模拟滤波器的最小阶数N和截止频率； [n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s') %计算阶数n和截止频率 %(4) 根据最小阶数N计算模拟低通原型滤波器的系统传递函数 [z0,p0,k0]=buttap(n); %设计归一化的巴特沃思模拟滤波器原型 [ba1,aa1]=zp2tf(z0,p0,k0); %由零极点得到传递函数 %(5)利用模拟域频率变换法将模拟低通原型滤波器转换为实际模拟滤波器的系统传递函数； [ba,aa]=lp2lp(ba1,aa1,Omgc); %变换为模拟低通滤波器 %(6)用脉冲响应不变法将模拟滤波器转换为数字滤波器 [bd,ad]=impinvar(ba,aa,Fs) %脉冲响应不变法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%数字滤波器特性 [H,w]=freqz(bd,ad); %求数字系统的频率特性 dbH=20*log10((abs(H)+eps)/max(abs(H))); subplot(2,2,1);plot(w/pi,abs(H)); ylabel('|H|');title('幅度响应');axis([0,1,0,1.1]); set(gca,'XTickMode','manual','XTick',[0,wp,ws,2*ws]); set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid subplot(2,2,2);plot(w/pi,angle(H)/pi); ylabel('\phi');title('相位响应');axis([0,1,-1,1]); set(gca,'XTickMode','manual','XTick',[0,wp,ws,2*ws]); set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid subplot(2,2,3);plot(w/pi,dbH);title('幅度响应(dB)'); ylabel('dB');xlabel('频率(\pi)');axis([0,1,-40,5]); set(gca,'XTickMode','manual','XTick',[0,wp,ws,12*ws]); set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);grid subplot(2,2,4);zplane(bd,ad);axis([-1.1,1.1,-1.1,1.1]);title('零极点图'); S1=filter(bd,ad,S); figure; plot(S1); title('原始信号');grid on; figure; Y = fft(S1,N); %做FFT变换 Ayy = (abs(Y)); %取模 Ayy=Ayy/(N/2); %换算成实际的幅度 Ayy(1)=Ayy(1)/2; F=([1:N]-1)/Ts/N; %换算成实际的频率值 plot(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果 title('幅度-频率曲线图'); grid on; figure; Pyy=[1:N/2]; for i=1:N/2 Pyy(i)=phase(Y(i)); %计算相位 Pyy(i)=Pyy(i)*180/pi; %换算为角度 end plot(F(1:N/2),Pyy(1:N/2)); %显示相位图 title('相位-频率曲线图'); grid on;
2)选择合适的窗函数设计FIR数字低通滤波器，要求：ωp=0.2 ，Rp=0.05dB；ωs=0.3 ，As=45dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。
代码：
主程序：
close all N=64;beta=7.865;n=1:N; wbo=boxcar(N); wtr=triang(N); whn=hanning(N); whm=hamming(N); wbl=blackman(N); wka=kaiser(N,beta); plot(n,wbo,'-',n,wtr,'*',n,whn,'+',n,whm,'.',n,wbl,'o',n,wka,'d'); axis([0,N,0,1.1]); legend('矩形','三角形','汉宁','哈明','布莱克曼','凯塞') figure wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp; N0=ceil(6.2*pi/deltaw); %6.6*pi是哈明窗过渡带宽，ceil函数取括号中数的最大整数 N=N0+mod(N0+1,2) %为实现FIR类型1偶对称滤波器，应确保N为奇数 windows=(hanning(N))';wc=(ws+wp)/2; %(2)由数字滤波器的理想频率响应H(ejω)求出其单位冲激响应hd(n)。 hd=ideal_lp(wc,N); %(3)计算数字滤波器的单位冲激响应h(n)=w(n)hd(n)。 b=hd.*windows; %(4) 检查设计的滤波器是否满足技术指标 [db,mag,pha,grd,w]=freqz_m(b,1);% 在0-pi区间的db相对幅度,mag绝对幅度,pha相位 %响应,grd群延时,w为501个采样点 n=0:N-1;dw=2*pi/1000; w1=1:wp/dw+1;%通带内的频率采样点 Rp=-(min(db(w1))) %计算通带内的波动 w2= ws/dw+1:501;%阻带内的频率采样点 As=-round(max(db(w2))) %计算最小阻带衰减，round(n)的意思是四舍五入 subplot(2,2,1);stem(n,b);axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应'); xlabel('n');ylabel('h(n)');subplot(2,2,2);stem(n,windows); axis([0,N,0,1.1]);title('窗函数特性');xlabel('n');ylabel('wd(n)'); subplot(2,2,3);plot(w/pi,db);axis([0,1,-80,10]);title('幅度频率响应'); xlabel('频率（单位：\pi）');ylabel('H(e^{j\omega})'); set(gca,'XTickMode','manual','XTick',[0,wp/pi,ws/pi,1]); set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid subplot(2,2,4);plot(w/pi,pha);axis([0,1,-4,4]);title('相位频率响应'); xlabel('频率（单位：\pi）');ylabel('\phi(\omega)'); set(gca,'XTickMode','manual','XTick',[0,wp/pi,ws/pi,1]); set(gca,'YTickMode','manual','YTick',[-3.1416,0,3.1416,4]);grid
子程序一 :idel_lp.midel_lp.m：
function hd=ideal_lp(wc,N) %点0到N-1之间的理想脉冲响应 %wc=截止频率（弧度）%N=理想滤波器的长度 tao=(N-1)/2; n=[0:(N-1)]; m=n-tao+eps; %加一个小数以避免0作除数 hd=sin(wc*m)./(pi*m);
子程序二:freqz_m``.m：
function [db,mag,pha,grd,w]=freqz_m(b,a) [H,w]=freqz(b,a,1000,'whole'); H=(H(1:501)); w=(w(1:501)); mag=abs(H); db=20*log10((mag+eps)/max(mag)); pha=angle(H); grd=grpdelay(b,a,w);
子程序三:freqz_m2``.m：
function [db,mag,pha,w]=freqz_m2(b,a) [H,w]=freqz(b,a,1000,'whole'); H=(H(1:1:501)); w=(w(1:1:501)); mag=abs(H); db=20*log10((mag+eps)/max(mag)); pha=angle(H); grd=grpdelay;
(3)分析信号S=3sin(2pi20t)+ 2.5sin(2pi60t)+ 2sin(2pi90t)，并用窗口函数法设计带阻滤波器，滤除60Hz的正弦波信号。要求画出：
① 信号滤波前的时域图及频谱图（坐标轴是实际幅度及实际频率）
② 双线性带阻滤波器的脉冲响应图，幅频响应曲线和相频响应曲线
③ 信号滤波后的时域图及频谱图（坐标轴是实际幅度及实际频率）
close all clc clear df=1;N=512; Ts=1/N/df;L=N*Ts;t=[0:Ts:L]; %采样时刻 %信号 S=3*sin(2*pi*20*t)+ 2.5*sin(2*pi*60*t)+ 2*sin(2*pi*90*t); %显示原始信号 plot(S); title('原始信号');grid on; figure; Y = fft(S,N); %做FFT变换 Ayy = (abs(Y)); %取模 Ayy=Ayy/(N/2); %换算成实际的幅度 Ayy(1)=Ayy(1)/2; F=([1:N]-1)/Ts/N; %换算成实际的频率值 plot(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果 title('幅度-频率曲线图'); grid on; figure; Pyy=[1:N/2]; for i=1:N/2 Pyy(i)=phase(Y(i)); %计算相位 Pyy(i)=Pyy(i)*180/pi; %换算为角度 end plot(F(1:N/2),Pyy(1:N/2)); %显示相位图 title('相位-频率曲线图'); grid on; figure; %(1)根据过渡带和阻带衰减指标选择窗函数的类型，估算滤波器的阶数N。 wp1=2*pi*20*Ts;wp2=2*70*pi*Ts; ws1=2*30*Ts*pi;ws2=2*80*Ts*pi; wp=[wp1,wp2];ws=[ws1,ws2]; deltaw=ws1-wp1; N0=ceil(6.2*pi/deltaw);N=N0+mod(N0+1,2); windows=(hanning(N))'; wc1=(ws1+wp1)/2;wc2=(ws2+wp2)/2; %(2)由数字滤波器的理想频率响应H(ejω)求出其单位冲激响应hd(n)。 hd=ideal_lp(wc1,N)+ideal_lp(pi,N)-ideal_lp(wc2,N);%建立理想带阻 %(3)计算数字滤波器的单位冲激响应h(n)=w(n)hd(n)。 b=hd.*windows; %(4) 检查设计的滤波器是否满足技术指标 [db,mag,pha,grd,w]=freqz_m(b,1); n=0:N-1;dw=2*pi/1000; w1=[1:wp1/dw+1,wp2/dw+1:501];%建立通带频率样点数组 Rp=-(min(db(floor(w1)))); %检验通带波动 w2= ws1/dw+1:ws2/dw+1; %建立阻带频率样点数组 As=-round(max(db(floor(w2)))); %检验最小阻带衰减 subplot(2,2,1);stem(n,b);axis([0,N,1.1*min(b),1.1*max(b)]); title('实际脉冲响应');xlabel('n');ylabel('h(n)'); subplot(2,2,2);stem(n,windows); axis([0,N,0,1.1]);title('窗函数特性');xlabel('n');ylabel('wd(n)'); subplot(2,2,3);plot(w/pi,db);axis([0,1,-150,10]); title('幅度频率响应');xlabel('频率（单位：\pi）');ylabel('H(e^{j\omega})'); %set(gca,'XTickMode','manual','XTick',[0,wp1/pi,ws1/pi,ws2/pi,wp2/pi,1]); %set(gca,'YTickMode','manual','YTick',[-150,-40,-3,0]);grid subplot(2,2,4);plot(w/pi,pha);axis([0,1,-4,4]); title('相位频率响应');xlabel('频率（单位：\pi）');ylabel('\phi(\omega)'); %set(gca,'XTickMode','manual','XTick',[0,wp1/pi,ws1/pi,ws2/pi,wp2/pi,1]); %set(gca,'YTickMode','manual','YTick',[-pi,0,pi]);grid S1=filter(b,1,S); figure; plot(S1); title('原始信号');grid on; figure; Y = fft(S1,N); %做FFT变换 Ayy = (abs(Y)); %取模 Ayy=Ayy/(N/2); %换算成实际的幅度 Ayy(1)=Ayy(1)/2; F=([1:N]-1)/Ts/N; %换算成实际的频率值 plot(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果 title('幅度-频率曲线图'); grid on; figure; Pyy=[1:N/2]; for i=1:N/2 Pyy(i)=phase(Y(i)); %计算相位 Pyy(i)=Pyy(i)*180/pi; %换算为角度 end plot(F(1:N/2),Pyy(1:N/2)); %显示相位图 title('相位-频率曲线图'); grid on;
三． 总结
1.数字滤波器的设计步骤。
答：设计数字滤波器一般分为3步；
（1）根据任务要求，确定滤波器性能指标
（2）用因果稳定的离散线性时不变系统函数逼近性能要求，即滤波器设计
设计滤波器的任务就是寻求一个物理上可实现的系统函数H(z)，使其频率响应H(z)满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止频率、通带衰减系数和阻带衰减系数
（3）用软硬件方法实现数字滤波器，用硬件时，主要用延迟器，加法器和乘法器等；软件实现时，用专用处理设备以及软件，或者通用计算机编写程序。
注意：数字滤波器频率响应以2π为周期，图1中给出0到2π的一个周期的特性，其余按周期延括规律重复。低通通带中心在2π的整数倍，高通通带中心在π的奇数倍。
2.脉冲响应不变法的原理及缺点。
答：脉冲响应不变法又称冲激响应不变法，使数字滤波器的单位脉冲响应序列h(n)等于模拟滤波器的单位冲激响应和ha(t)的采样值，即： ，其中，T为采样周期。在脉冲响应不变法中，数字角频率和模拟角频率的变换关系为： ，可见，Ω和ω之间的变换关系为线性的。
MATLAB中用函数impinvar实现从模拟滤波器到数字滤波器的脉冲响应不变映射，调用格式为：
[B,A]=impinvar(b,a,fs) [B,A]=impinvar(b,a)
其中，b、a分别为模拟滤波器的分子和分母多项式系数向量；fs为采样频率（Hz），缺省值fs=1Hz；B、A分别为数字滤波器分子和分母多项式系数向量。
脉冲响应不变法是将系统从s平面映射到z平面的一种变换方法，使数字滤波器的单位脉冲响应h(n)模仿模拟滤波器的冲激响应ha(t)，其变换关系式为z=esT。由于esT是一个周期函数，因而s平面虚轴上每一段2π/T的线段都映射到z平面单位圆上一周。由于重叠映射，因而冲激响应不变法是一种多值映射关系。数字滤波器的频率响应是原模拟滤波器的频率响应的周期延拓（如图3）所示。只有当模拟滤波器的频率响应是有限带宽，且频带宽度|Ω|≤(π/T)=ΩS/2，才能避免数字滤波器的频率响应发生混叠现象。脉冲响应不变法只适用于限带的模拟滤波器，对于高频区幅频特性不等于零的高通和带阻滤波器不适用。
3.窗函数设计FIR数字滤波器的原理。
答：略，自己在我上传的资源文档里面找吧。
