数字信号处理实验三、四.docx

实验三、IIR数字滤波器的设计
测控1201 侍小青 12054109
一、 思考题：
1、双线性变换法中Ω和ω之间的关系是非线性的，在实验中你注意到这种非线性关系了吗？从哪几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？  
答：在双线性变化法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟滤波器经双线性变换后，得到的数字滤波器不再保持原有的线性相位了，在每一幅使用了双线性变换的图中，可以看到在采样频率一半处，幅度为零，这显然不是线性变换能够产生的，这是由于双线性变换将模拟域中的无穷远点映射到了改点处。  
采用双线性变化法设计的巴特沃斯型和切比雪夫数字型滤波器，可以观察到这种非线性关系。
2、能否利用公式Hz=Hs|s=1TlnZ完成脉冲响应不变法的数字滤波器设计？为什么？ 
答：IIR数字滤波器的设计实际上是求解滤波器的系数,它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。如果在S平面上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。但是它的缺点是，存在频率混叠效应，故只适用于阻带的模拟滤波器。
3、谈谈双线性变换法的特点，简述双线性变换法设计滤波器的全过程。
答：与脉冲响应不变法相比，双线性变换的主要优点：靠频率的严重非线性关系得到S平面与Z平面的单值一一对应关系，整个jΩ轴单值对应于单位圆一周，其中ω和Ω为非线性关系。在零频率附近，Ω～ω接近于线性关系，Ω进一步增加时，ω增长变得缓慢，(ω终止于折叠频率处)，所以双线性变换不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象。
双线性变换法的缺点：Ω与ω的非线性关系，导致数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变，(使数字滤波器与模拟滤波器在响应与频率的对应关系上发生畸变)。另外，一个线性相位的模拟滤波器经双线性变换后，滤波器就不再有线性相位特性。
双线性变换法设计滤波器的全过程：1.确定数字滤波器的性能指标：通带临界频率cf、阻带临界频率rf、通带波动d、阻带内的最小衰减At、采样周期T、采样频率fs； 
2.确定相应的数字角频率；
3.计算经过预畸的相应模拟低通原型的频率
4.根据Ωc和Ωr计算模拟低通原型滤波器的阶数N，并求得低通原型的传递函数Ha（s）；
5.用上面的双线性变换公式代入Ha（s），求出所设计的传递函数H(Z);
6.分析滤波器特性，检查其指标是否满足要求。
二、实验原理
 1． 脉冲响应不变法 
用数字滤波器的单位脉冲响应序列h（a）模仿模拟滤波器的冲激响应ha（t）,让h（a）正好等于ha（t）的采样值，即h（n）=ha（nt），其中T为采样间隔，如果以Ha（s）及H(z)分别表示ha（t）的拉式换及h（n）的Z变换，则-¥==+=p
 2．双线性变换法 
S平面与z平面之间满足以下映射关系： 
s平面的虚轴单值地映射于z平面的单位圆上，s平面的左半平面完全映射到z平面的单位圆内。双线性变换不存在混叠问题。双线性变换是一种非线性变换，这种非线性引起的幅频特性畸变可通过预畸而得到校正。
三、实验内容
1、wc=2*1000*tan(2*pi*300/(2*1000));
wt=2*1000*tan(2*pi*200/(2*1000)); %对临界频率进行预畸变；
[N,wn]=cheb1ord(wc,wt,0.8,20,'s'); %给定通带边界频率wc和阻带边界频率wr、通带波动0.8dB和最小阻带衰减20dB，求满足指标的模拟滤波器的最低阶数N和通带边界频率wn；
[B,A]=cheby1(N,0.8,wn,'high','s'); %给定模拟滤波器的最低阶数N、通带边界频率wn和通带波动0.8dB，设计N阶模拟巴特沃斯高通滤波器，B和A分别表示系统函数的分子和分母多项式的系数；
[num,den]=bilinear(B,A,1000); %给定模拟滤波器系统函数H(s)=B(s)/A(s)和采样频率1000Hz，根据双线性变换法求出数字滤波器的系统函数H(z)=B(z)/A(z)，num和den分别为数字滤波器的分子分母多项式系数；
[h,w]=freqz(num,den); %计算以(num,den)为参数的滤波器的频率向量w和复频率响应向量h；
f=w/pi*500;
plot(f,20*log10(abs(h)));
axis([0,500,-80,10]);
grid;
xlabel('频率/Hz');
ylabel('幅度/dB');

2、T = 0.001;fs = 1000;fc = 200;fr = 300; 
wp1 = 2*pi*fc;wr1 = 2*pi*fr; 
[N1,wn1] = buttord(wp1,wr1,1,25,'s') 
[B1,A1] = butter(N1,wn1,'s'); 
[num1,den1] = impinvar(B1,A1,fs);%脉冲响应不变法 
[h1,w] = freqz(num1,den1); 
wp2 = 2*fs*tan(2*pi*fc/(2*fs)) 
wr2 = 2*fs*tan(2*pi*fr/(2*fs)) 
[N2,wn2] = buttord(wp2,wr2,1,25,'s') 
[B2,A2] = butter(N2,wn2,'s'); 
[num2,den2] = bilinear(B2,A2,fs);%双线性变换法 
[h2,w] = freqz(num2,den2); 
f = w/(2*pi)*fs; 
plot(f,20*log10(abs(h1)),'-.',f,20*log10(abs(h2)),'-'); 
axis([0,500,-80,10]);
grid;xlabel('频率/Hz ');
ylabel('幅度/dB') 

从图中可见，通带边界和阻带边界分别为200hz,300hz,衰减量也满足为25Db。 
总结： 脉冲响应不变法的一个重要特点是频率坐标的变换是线性的，ω＝ΩΤ，ω与Ω是线性关系：在某些场合，要求数字滤波器在时域上能模仿模拟滤波器的功能时，如要实现时域冲激响应的模仿，一般使用脉冲响应不变法。脉冲响应不变法的最大缺点：有频谱周期延拓效应，因此只能用于带限的频响特性，如衰减特性很好的低通或带通,而高频衰减越大，频响的混淆效应越小，至于高通和带阻滤波器,由于它们在高频部分不衰减，因此将完全混淆在低频响应中，此时可增加一保护滤波器，滤掉高于 /2sW 的频带，再用脉冲响应不变法转换为数字滤波器，这会增加设计的复杂性和滤波器阶数，只有在一定要满足频率线性关系或保持网络瞬态响应时才采用。双线性变换法的主要优点是S平面与Z平面一一单值对应，s平面的虚轴(整个jΩ)对应于Z平面单位圆的一周，S平面的Ω=0处对应于Z平面的ω=0处， Ω= ∞处对应于Z平面的ω= π处,即数字滤波器的频率响应终止于折叠频率处，所以双线性变换不存在混迭效应。 
双线性变换缺点: Ω与ω成非线性关系，导致： 
a. 数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变，(使数字滤波器与模拟滤波器在响应与频率的对应关系上发生畸变)。      
b. 线性相位模拟滤波器经双线性变换后，得到的数字滤波器为非线性相位。
c.要求模拟滤波器的幅频响应必须是分段恒定的，故双线性变换只能用于设计低通、高通、带通、带阻等选频滤波器
3、
分析：
N=9，为九阶巴特沃思低通滤波器，从图中可以看出通带波动和阻带衰减都满足设计要求。
N=5，为五阶切比雪夫低通滤波器，从图中可以看出通带波动和阻带衰减都满足设计要求。
N=4，为四阶椭圆型数字低通滤波器，从图中可以看出通带波动和阻带衰减都满足设计要求。
4、wp1 = 2*pi*2000;wp2 = 2*pi*3000; 
ws1= 2*pi*1500;ws2= 2*pi*6000;
[N1,wn1] = buttord([wp1 wp2],[ws1 ws2],3, 20 ,'s') 
[B1,A1] = butter(N1,wn1,'s'); 
[num1,den1] = impinvar(B1,A1,30000) 
[h1,w] = freqz(num1,den1); 
w1=2*30000*tan(2*pi*2000/(2*30000)); 
w2=2*30000*tan(2*pi*3000/(2*30000)); 
wr1=2*30000*tan(2*pi*1500/(2*30000)); 
wr2=2*30000*tan(2*pi*6000/(2*30000)); 
[N,wn]=buttord([w1 w2],[wr1 wr2],3,20,'s'); 
[B,A]=butter(N,wn,'s'); 
[num,den]=bilinear(B,A,30000); 
[h2,w]=freqz(num,den); 
f=w/(2*pi)*30000; 
plot(f,20*log10(abs(h1)),'-.',f,20*log10(abs(h2)),'-'); 
axis([0,15000,-60,10]); 
grid;xlabel('频率/Hz');ylabel('幅度/dB') 

5、w1=2*10000*tan(2*pi*1000/(2*10000)); 
w2=2*10000*tan(2*pi*2000/(2*10000)); 
wr1=2*10000*tan(2*pi*500/(2*10000)); 
wr2=2*10000*tan(2*pi*3000/(2*10000)); 
[N,wn]=cheb1ord([wr1 wr2],[w1 w2],3,18,'s'); 
[B,A]=cheby1(N,3,wn,'stop','s'); 
[num,den]=bilinear(B,A,10000) 
[h,w]=freqz(num,den); 
f=w/(2*pi)*10000; 
plot(f,20*log10(abs(h))); 
axis([0,5000,-120,10]); 
grid;xlabel('频率/Hz');ylabel('幅度/dB')

N=4，为四阶切比雪夫带阻滤波器，由图知满足设计要求。
四、实验总结
通过这次实验，我掌握了双线性变换法及脉冲相应不变法设计IIR数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。
实验四、FIR数字滤波器的设计

一、思考题
（1）定性地说明用本实验程序设计的 FIR 滤波器的3dB截止频率在什么位置？它等于理想频率响应Hd(ejω)的截止频率吗？ 
答：从图形中看出，本实验设计的FIR滤波器的3dB截止频率在0.4π或0.3π和0.5π，基本等于理想频率响应的截止频率。  
（2）如果没有给定 h(n)的长度 N，而是给定了通带边缘截止频率ωc和阻带临界频率ωp，以及相应的衰减，能根据这些条件用窗函数法设计线性相位FIR低通滤波器吗？ 
答：可以；h(n)的长度对实验的影响不是很大，N的取值对h(n)的值有影响但对h(n)的幅频特性影响不大。
二、实验内容
1、clear all;
N=15; %N=15的汉明窗
h=fir1(N-1,[0.3 0.5],'bandpass',hanning(N));
figure(1);
freqz(h,1);
title('N=15,汉宁窗');
N=45; %N=45的汉明窗
h=fir1(N-1,[0.3 0.5],'bandpass',hanning(N));
figure(2);
freqz(h,1);
title('N=45,汉宁窗');

clear all; %矩形窗
N=15;h=fir1(N-1,[0.3 0.5],'bandpass',kaiser(N,0)); %N=15的矩形窗
[h1,w1]=freqz(h,1);
subplot(2,1,1);plot(w1/pi,20*log10(abs(h1))); %画出其幅频特性
axis([0,1,-80,10]);grid;
xlabel('归一化频率/\pi');ylabel('幅度/dB');title('N=15,矩形窗'); %为坐标标注x、y轴以及标题
N=45;h=fir1(N-1,[0.3 0.5],'bandpass',kaiser(N,0)); %N=45的矩形窗
[h1,w1]=freqz(h,1);
subplot(2,1,2);plot(w1/pi,20*log10(abs(h1)));
axis([0,1,-80,10]);grid;
xlabel('归一化频率/\pi');ylabel('幅度/dB');title('N=45,矩形窗');
clear all; %布莱克曼窗
N=15;h=fir1(N-1,[0.3 0.5],'bandpass',blackman(N)); %N=15的布莱克曼窗
[h1,w1]=freqz(h,1);
subplot(2,1,1);plot(w1/pi,20*log10(abs(h1)));
axis([0,1,-80,10]);grid;
xlabel('归一化频率/\pi');ylabel('幅度/dB');title('N=15,布莱克曼窗');
N=45;h=fir1(N-1,[0.3 0.5],'bandpass',blackman(N)); %N=45的布莱克曼窗
[h1,w1]=freqz(h,1);
subplot(2,1,2);plot(w1/pi,20*log10(abs(h1)));
axis([0,1,-80,10]);grid;
xlabel('归一化频率/\pi');ylabel('幅度/dB');title('N=45,布莱克曼窗');

分析：1、矩形窗函数具有最窄的主瓣宽度，但有最大的旁瓣峰值；汉明窗函数的主瓣稍宽，而旁瓣较小；布莱克曼窗函数则更甚之。矩形窗设计的滤波器过渡带最窄，但是阻带最小衰减也最差；布莱克曼窗设计的滤波器阻带衰减最好，过度带最宽，约为矩形窗设计的的三倍。汉明窗设计的滤波器处于矩形窗和布莱克曼窗之间。
2、观察它们的实际3dB和20dB带宽，发现N=15时，其3dB带宽约为0.2pi,20dB带宽约为0.45pi；N=45时，其3dB带宽约为0.16pi,20dB带宽约为0.3pi;可见N增大，其3dB带宽和20dB带宽分别减小，滤波器特性变好，过渡带变陡，幅频曲线显示其通带较平缓，波动小，阻带衰减大。相频特性曲线显示其相位随频率变化也变大。
3、同一N值，分别用矩形窗，汉宁窗，汉明窗，布莱克曼窗设计滤波器时，主瓣宽度逐渐增大,过渡带变宽，但阻带衰减性能变好；N增加，主瓣变窄，旁瓣的分量增加，过渡带变陡，起伏震荡变密。   
加窗处理对滤波器的频率响应会产生以下主要影响：  
（1）使理想特性不连续的边沿加宽，形成一过渡带，过渡带的宽度取决于窗函数频谱的主瓣宽度。  
（2）在过渡带两旁产生肩峰和余振，它们取决于窗函数频谱的旁瓣；旁瓣越多，余振也越多；旁瓣相对值越大，肩峰则越强。  
（3）增加截断长度N，只能缩小窗函数频谱的主瓣宽度而不能改变旁瓣的相对值；旁瓣与主瓣的相对关系只决定于窗函数的形状。因此增加N，只能相对应减小过渡带宽。而不能改变肩峰值。肩峰值的大小直接决定通带内的平稳和阻带的衰减，对滤波器性能有很大关系。
4、clear all; 
N=40;  
f = [0 0.2 0.2 0.4 0.4 0.6 0.6 0.8 0.8 1]  
a = [0 0 1 1 0 0 1 1 0 0]  
beta=4; 
h = fir2(N-1,f,a,kaiser(N,beta)); 
[h1,w1]=freqz(h,1);  
subplot(3,1,1); 
plot(w1/pi,20*log10(abs(h1))); 
axis([0,1,-80,10]); 
grid;
xlabel('归一化频率/\pi');  
ylabel('幅度/dB'); 
title('beta=4 时凯塞窗专用线性相位滤波器'); 
beta=6; 
h = fir2(N-1,f,a,kaiser(N,beta)); 
[h1,w1]=freqz(h,1);  
subplot(3,1,2); 
plot(w1/pi,20*log10(abs(h1))); 
axis([0,1,-80,10]); 
grid;  
xlabel('归一化频率/\pi');  
ylabel('幅度/dB'); 
title('beta=6 时凯塞窗专用线性相位滤波器'); 
beta=10; 
h = fir2(N-1,f,a,kaiser(N,beta)); 
[h1,w1]=freqz(h,1);  
subplot(3,1,3); 
plot(w1/pi,20*log10(abs(h1))); 
axis([0,1,-80,10]); 
grid;  
xlabel('归一化频率/\pi');  
ylabel('幅度/dB'); 
title('beta=10 时凯塞窗专用线性相位滤波器');

分析：β越大，w(n)窗越窄，频谱的旁瓣越小，但主瓣宽度也相应增加，过渡带变宽，相位特性变好。
5、clear all; N=40; 
Hk=[zeros(1,3) 0.5 ones(1,5) 0.5 zeros(1,1) 0.5 ones(1,5) 0.5...
zeros(1,5) -0.5 -ones(1,5) -0.5 zeros(1,1) -ones(1,5) -0.5 zeros(1,3)]; 
k=0:N-1; 
hn=real(ifft(Hk.*exp(-j*pi*(N-1)*k/N))); 
[H w]=freqz(hn, 1);  
plot(w/pi, 20 *log10(abs(H)));  
axis([0 1 -80 10])  grid;  
xlabel('归一化频率/\pi')  
ylabel('幅度/dB') 
title('频率采样法设计专用线性相位滤波器'); 

6、clear all; N=40;  
f=[0 0.15 0.2 0.4 0.45 0.55 0.6 0.8 0.85 1]   
a=[0 0 1 1 0 0 1 1 0 0]   
wt=[2 1 2 1 2]  
b=remez(N-1,f,a,wt);  
[h,w]=freqz(b,1);
plot(w/pi,20*log10(abs(h)));
axis([0 1 -70 10])  
grid; 
xlabel('归一化频率/\pi')  
ylabel('幅度/dB') 
title('雷米兹交替算法设计专用线性相位滤波器');

7、clear all;  
fedge=[500 800]   
mval=[0 1]   
dev=[0.01 0.109]   
fs=5000; 
[N,fpts,mag,wt]=remezord(fedge,mval,dev,fs);  
b=remez(N,fpts,mag,wt);  
[h,w]=freqz(b,1);  
plot(w*2500/pi,20*log10(abs(h)));  
axis([0 2500 -80 10])  grid;  
xlabel('频率/Hz')  
ylabel('幅度/dB')
title('雷米兹交替算法设计线性相位高通FIR数字滤波器');

三、实验总结
1. 通过本实验掌握窗函数设计了FIR滤波器，掌握了其原理和方法。  
2. 了解了线性相位FIR滤波器的幅频特性和相位特性。  
3. 窗函数首先对滤波器提出性能指标，通过不同窗对滤波器的响应，充分理解了各个窗的功能和用法。会用窗函数设计FIR滤波器。
这部分知识比较欠缺，需要自学。上课介绍的比较少，所以懂得比较少，做的比较吃力，但以后会好好学习。