实验1离散系统的时域分析

1、数字信号处理实验 学号： 姓名：实验1 离散系统的时域分析一、实验目的：加深对离散系统的差分方程、单位抽样响应和卷积分析方法的理解。 二、实验原理：离散系统Discrete-time systmex(n)y(n)其输入、输出关系可用以下差分方程描述：输入信号分解为冲激信号，系统单位抽样序列h（n），则系统响应为如下的卷积计算式：当时，h(n)是有限长度的（n：0，M），称系统为FIR系统；反之，称系统为IIR系统。三 、实验内容编制程序求解下列两个系统的单位抽样响应，并绘出其图形。（1）（2）（1）源程序：N=21;b=1 -1;a=1 0.75 0.125;x=1 zeros(1,N-1);

2、n=0:1:N-1;y=filter(b,a,x);stem(n,y);xlabel('n');ylabel('幅度');title('单位抽样响应');图形：（2）源程序：N=20;d=0 0.25 0.25 0.25 0.25;c=1;x=1 zeros(1,N-1);n=0:1:N-1;y=filter(d,c,x);stem(n,y);xlabel('n');ylabel('幅度');title('单位抽样响应');图形如下： 实验2 离散系统的频率响应分析和零、极点分布一、实验目的：加深对

3、离散系统的频率响应分析和零、极点分布的概念理解。二、实验原理：离散系统的时域方程为其变换域分析方法如下：系统函数为 分解因式 ，其中 和 称为零、极点。 在MATLAB中，可以用函数z，p，K=tf2zp（num，den）求得有理分式形式的系统函数的零、极点，用函数zplane（z，p）绘出零、极点分布图；也可以用函数zplane（num，den）直接绘出有理分式形式的系统函数的零、极点分布图。使h=freqz(num,den,w)函数可求系统的频率响应，w是频率的计算点，如w=0:pi/255:pi, h是复数，abs(h)为幅度响应，angle(h)为相位响应。另外，在MATLAB中，可以

4、用函数 r，p，k=residuez（num，den）完成部分分式展开计算；可以用函数sos=zp2sos（z，p，K）完成将高阶系统分解为2阶系统的串联。三、实验内容练习1 求下列直接型系统函数的零、极点，并将它转换成二阶节形式 解 用MATLAB计算程序如下： num=1 -0.1 -0.3 -0.3 -0.2;den=1 0.1 0.2 0.2 0.5;z,p,k=tf2zp(num,den);disp('零点');disp(z);disp('极点');disp(p);disp('增益系数');disp(k);sos=zp2sos(z,p,

5、k);disp('二阶节');disp(real(sos);zplane(num,den)输入到“num”和“den”的分别为分子和分母多项式的系数。计算求得零、极点增益系数和二阶节的系数：零点 0.9615 -0.5730 -0.1443 + 0.5850i -0.1443 - 0.5850i  极点 0.5276 + 0.6997i 0.5276 - 0.6997i -0.5776 + 0.5635i -0.5776 - 0.5635i  增益系数 1  二阶节 1.0000 -0.3885 -0.5509 1.0000 1.1552 0.65

6、11 1.0000 0.2885 0.3630 1.0000 -1.0552 0.7679  系统函数的二阶节形式为： 极点图如右图。    2 差分方程 所对应的系统的频率响应。 解：差分方程所对应的系统函数为用MATLAB计算的程序如下：k=256;num=0.8 -0.44 0.36 0.02;den=1 0.7 -0.45 -0.6;w=0:pi/k:pi;h=freqz(num,den,w);subplot(2,2,1);plot(w/pi,real(h);gridtitle('实部')xlabel('omega/pi')

7、;ylabel('幅度')subplot(2,2,2);plot(w/pi,imag(h);gridtitle('虚部')xlabel('omega/pi');ylabel('Amplitude')subplot(2,2,3);plot(w/pi,abs(h);gridtitle('幅度谱')xlabel('omega/pi');ylabel('幅值')subplot(2,2,4);plot(w/pi,angle(h);gridtitle('相位谱')xlabel(&

8、#39;omega/pi');ylabel('弧度')                   3求系统 的零、极点和幅度频率响应和相位响应。要求：绘出零、极点分布图，幅度频率响应和相位响应曲线。 源程序：num=0.0528 0.0797 0.1295 0.1295 0.0797 0.0528;den=1 -1.8107 2.4947 -1.8801 0.9537 -0.2336;subplot(3,1,

9、1); z,p,k=tf2zp(num,den);disp('零点');disp(z);disp('极点');disp(p);disp('增益系数');disp(k);sos=zp2sos(z,p,k);disp('二阶节');disp(real(sos);zplane(num,den)title('零极点分布')k=256;m=0.0528 0.0797 0.1295 0.1295 0.0797 0.0528;n=1 -1.8107 2.4947 -1.8801 0.9537 -0.2336;w=0:pi/k:p

10、i;h=freqz(m,n,w);subplot(3,2,3);plot(w/pi,real(h);gridtitle('实部')xlabel('omega/pi');ylabel('幅度')subplot(3,2,4);plot(w/pi,imag(h);gridtitle('虚部')xlabel('omega/pi');ylabel('Amplitude')subplot(3,2,5);plot(w/pi,abs(h);gridtitle('幅度谱')xlabel('om

11、ega/pi');ylabel('幅值')subplot(3,2,6);plot(w/pi,angle(h);gridtitle('相位谱')xlabel('omega/pi');ylabel('弧度') 实验结果：零点 -1.0000 -0.3018 + 0.9534i -0.3018 - 0.9534i 0.0471 + 0.9989i 0.0471 - 0.9989i极点 0.2788 + 0.8973i 0.2788 - 0.8973i 0.3811 + 0.6274i 0.3811 - 0.6274i 0.491

12、0 增益系数 0.0528二阶节 0.0528 0.0528 0 1.0000 -0.4910 0 1.0000 0.6036 1.0000 1.0000 -0.7622 0.53891.0000 -0.0941 1.0000 1.0000 -0.5575 0.8829图形：实验3 FFT算法的应用一 、实验目的：加深对离散信号的DFT的理解及其FFT算法的运用。二 、实验原理：N点序列的DFT和IDFT变换定义式如下：， 利用旋转因子具有周期性，可以得到快速算法（FFT）。 在MATLAB中，可以用函数X=fft（x，N）和x=ifft（X，N）计算N点序列的DFT正、反变换。三、实验内容1

13、. 对连续的单一频率周期信号 按采样频率采样，截取长度N分别选N =20和N =16，观察其DFT结果的幅度谱。 解 此时离散序列： 即N=8。用MATLAB计算并作图，函数fft用于计算离散傅里叶变换DFT，程序如下： q=8;k1=0:1:19;xa1=sin(2*pi*k1/q);subplot(2,2,1)plot(k1,xa1)xlabel('k');ylabel('x1(k)');xm1=fft(xa1);xm1=abs(xm1);subplot(2,2,2)stem(k1,xm1)xlabel('m');ylabel('X1

14、(m)');k2=0:1:15;xa2=sin(2*pi*k2/q);subplot(2,2,3)plot(k2,xa2)xlabel('k');ylabel('x2(k)');xm2=fft(xa2);xm2=abs(xm2);subplot(2,2,4)stem(k2,xm2)xlabel('m');ylabel('X2(m)');(d)(c)(b)(a)         图2.1计算结果示于图2.1，(c) 和(d) 分别是N=16时的截

15、取信号和DFT结果，由于截取了两个整周期，得到单一谱线的频谱；而(a) 和(b)分别是N=20时的截取信号和DFT结果，由于截取了两个半周期，不再是单一谱线的频谱，而是含有丰富的谐波分量，与N=16时的差别很大；上述频谱的误差主要是由于有限长序列的离散傅里叶变换是把有限长序列作为周期序列的一个周期来表示的，隐含有周期性意义。 2. 分别计算16点序列 的16点和32点DFT,绘出幅度谱图形。源程序：q=16;n1=0:1:15;xn1=sin(5*pi*n1/q);subplot(2,2,1)plot(n1,xn1)xlabel('n');ylabel('x1(n)&#

16、39;);xm1=fft(xn1);xm1=abs(xm1);subplot(2,2,2)stem(n1,xm1)xlabel('m');ylabel('X1(m)');n2=0:1:31;xn2=sin(5*pi*n2/q);subplot(2,2,3)plot(n2,xn2)xlabel('n');ylabel('x2(n)');xm2=fft(xn2);xm2=abs(xm2);subplot(2,2,4)stem(n2,xm2)xlabel('m');ylabel('X2(m)');实验结果

17、：实验4 无限冲激响应数字滤波器设计一、实验目的掌握双线性变换法及冲激响应不变法设计IIR数字滤波器的具体设计方法及其原理；熟悉用双线性变换法及冲激响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。二、实验编程函数 在MATLAB中，可以用下列函数辅助设计IIR数字滤波器：1) 利用buttord和cheblord可以确定低通原型巴特沃斯和切比雪夫滤波器的阶数和截止频率；2) num,den=butter（N,Wn）(巴特沃斯)；3) num,den=cheby1（N,Wn）,num,den=cheby2（N,Wn）(切比雪夫1型和2型)可以进行滤波器的设计；4) lp2hp，lp2

18、bp，lp2bs可以完成低通滤波器到高通、带通、带阻滤波器的转换；5) 使用bilinear可以对模拟滤波器进行双线性变换，求得数字滤波器的传输函数系数；6) 利用impinvar可以完成冲激响应不变法的模拟滤波器到数字滤波器的转换。注： 双线性变换法通过将数字频率的取值范围从0到对应到模拟频率，也就对应于模拟域中所有可能的频率值。双线性变换法不会出现频率混叠，但非线性关系却导致数字滤波器的频率响应不能逼真地模仿模拟滤波器的频率响应。冲激响应不变法通过选择满足设计要求的模拟滤波器单位冲激响应h（t）的采样值的h（n），得到的被采样的冲激响应将给出与原模拟滤波器非常相近的滤波器形状。由于该方法不

19、可避免的要发生频率混叠现象，所以只适合设计低通和带通滤波器。 三、实验内容 1.设采样周期T=250s，用冲激响应不变法和双线性变换法设计一个三阶巴特沃斯滤波器，其3dB边界频率为fc =1kHz。 B,A=butter(3,2*pi*1000,'s'); num1,den1=impinvar(B,A,4000); h1,w=freqz(num1,den1); B,A=butter(3,2/0.00025,'s'); num2,den2=bilinear(B,A,4000); h2,w=freqz(num2,den2); f=w/pi*2000; plot(f,

20、abs(h1),'-.',f,abs(h2),'-'); grid; xlabel('频率/Hz ') ylabel('幅值/dB') 程序中第一个butter的边界频率2×1000，为冲激响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图1给出了这两种设计方法所得到的频响，虚线为冲激响应不变法的结果；实线为双线性变换法的结果。冲激响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。 图12. 设计一巴特沃斯带通滤

21、波器，其通带边界频率分别为f2=110kHz和f1=90kHz处的最大衰减小于3dB，在阻带边界频率分别为f3 =80kHz 和f4=120kHz处的最小衰减大于10dB，采样频率fs=400Hz。 w1=2*400*tan(2*pi*90/(2*400); w2=2*400*tan(2*pi*110/(2*400); wr1=2*400*tan(2*pi*80/(2*400); wr2=2*400*tan(2*pi*120/(2*400); N,wn=buttord(w1 w2, wr1 wr2,3,10,'s'); B,A=butter(N,wn,'s');

22、 num,den=bilinear(B,A,400); h,w=freqz(num,den); f=w/pi*200; plot(f,20*log10(abs(h); axis(40,160,-30,10); grid; xlabel('频率/kHz') ; ylabel('幅度/dB') 图23. 利用MATLAB编程设计一个数字带通滤波器，指标要求如下：通带边缘频率：，通带衰减。阻带边缘频率：，最小阻带衰减 。给出IIR数字滤波器参数，绘出幅度和相位频响曲线，讨论实现形式和特点。源程序：wp=0.45 0.65; ws=0.3 0.75; Ap=1;As=4

23、0;n1,wn1=buttord(wp,ws,Ap,As);num,den=butter(n1,wn1); w=0:pi/255:pi; h=freqz(num,den,w); g=20*log10(abs(h); pha=angle(h); subplot(1,2,1)plot(w/pi,g); grid xlabel('频率/kHz');ylabel('幅度');title('数字带通滤波器幅频曲线')subplot(1,2,2);plot(w/pi,pha);gridxlabel('频率/kHz');ylabel('

24、相位');title('数字带通滤波器相位曲线'); 实验结果：实现形式和特点分析：实验5 有限冲激响应数字滤波器设计一 、实验目的：（1） 掌握用窗函数法设计FIR数字滤波器的原理和方法。（2） 了解各种窗函数对滤波特性的影响。二、实验原理：如果所希望的滤波器的理想频率响应函数为 ,则其对应的单位抽样响应为窗函数设计法的基本原理是用有限长单位抽样响应序列h(n)逼近hd(n)，由于hd(n)往往是无限长序列，且是非因果的，所以用窗函数w(n)将hd(n)截断，并进行加权处理，得到:h(n)就作为实际设计的FIR数字滤波器的单位抽样响应序列，其频率响应函数为：式中，N为

25、所选窗函数的长度 在MATLAB中，可以用b=fir1(M,Wc,ftype,taper) 等函数辅助设计FIR数字滤波器。M代表滤波器阶数(单位抽样响应h（n）的长度N=M+1)；Wc代表滤波器的截止频率(对归一化频率)，当设计带通和带阻滤波器时，Wc为双元素相量；ftype代表滤波器类型，如high高通，stop带阻等；taper为窗函数类型，默认为海明窗；窗函数用blackman, hamming, hanning ， kaiser产生。三、实验内容1.设计一FIR低通滤波器，通带截止频率 ，通带衰减不大于0.25dB，阻带截止频率 ，阻带衰减不小于50dB。选择一个合适的窗函数，绘出所

26、设计的滤波器的频率响应图。wp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;M=ceil(6.6*pi/tr_width);n=0:1:M-1;wc=(ws+wp)/2;hd=ideal_lp(wc,M);w_ham=(hamming(M)'h=hd.*w_ham;H,w=freqz(h,1,1000,'whole');H=(H(1:1:501)'w=(w(1:1:501)'mag=abs(H);db=20*log10(mag+eps)/max(mag);delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/del

27、ta_w+1);disp('实际通带波动为',num2str(Rp);As=-round(max(db(ws/delta_w+1:1:501);disp('最小阻带衰减为',num2str(As); subplot(2,2,1);stem(n,hd);title('理想抽样响应');axis(0 M-1 -0.1 0.3);xlabel('n');ylabel('hd(n)');subplot(2,2,2);stem(n,w_ham);title('海明窗');axis(0 M-1 0 1.1);xlabel('n');ylabel('w(n)');subplot(2,2,3);stem(n,h);title('实际脉冲响应');axis(0 M-1 -0.1 0.3);xlabel('n');ylabel('h(n)');subplot(2,2,4);p