数字信号处理实验指导书

1、数字信号处理实验指导书2实验一 离散时间与系统的傅立叶分析一、实验目的用傅立叶变换对信号和系统进行频域分析。二、实验原理对信号进行频域分析就是对信号进行傅立叶变换。对系统进行频域分析即对它的单位脉冲响应进行傅立叶变换，得到系统的传输函数。也可以由差分方程经过；傅立叶变换直接 求它的传输函数。传输函数代表的就是系统的频率响应特性。但传输函数是的连续函数，计算机只能计算出有限个离散频率点的传输函数值，因此得到传输函数以后，应该在 02 之间取许多点，计算这些点的传输函数的值，并取它们的包络，该包络才是需要的频率特性。当然，点数取得多一些，该包络才能接近真正的频率特性。 注意：非周期信号的频率特性是

2、的连续函数，而周期信号的频率特性是离散谱，它们的计算公式不一样，响应的波形也不一样。 三、实验内容 1已知系统用下面差分方程描述： y(n)=x(n)十ay(n一1) 试在a095和a一0.5两种情况下用傅立叶变换分析系统的频率特性。要求写出系统的传输函数和幅度响应，并打印|H(ejw)|曲线。 2已知两系统分别用下面差分方程描述： y1(n)x(n)+x(n一1) y2(n)x(n)一x(n一1) 试分别写出它们的传输函数和幅度响应，并分别打印|H(ejw)|曲线。 3已知信号x(n)R3(n)，试分析它的频域特性，要求打印|H(ejw)|曲线。 4假设x(n)a(n)，将x(n)以2为周期

3、进行周期延拓，得到x(n)，试分析它的频率特性，并画出它的幅频特性。 四、实验用MATLAB函数介绍 1abs 功能：求绝对值(复数的模)。 yabs(x)：计算实数x的绝对值。当x为复数时得到x的模(幅度值)。 当x为向量时，计算其每个元素的模，返回模向量y。 2angle 功能：求相角。 Ph=angle(x)：计算复向量x的相角(rad)。Ph值介于 - 和+之间. 3freqz： 计算数字滤波器H(z)的频率响应。 Hfreqz(B，A，w)：计算由向量w指定的数字频率点上数字滤波器H(z)的频率响应H(ejw)，结果存于H向量中.向量B和A分别为数字滤波器系统函数H(z)的分子和分母

4、多项式系数。 H，wfreqz(B，A，M，whole)：计算出M个频率点上酌频率响应，存放在，H向量中，M个频率存放在向量w中。厅reqz函数自动将这M个频点均匀设置在频率范围 0,2 上。缺省whole时，M个频点均匀设置在频率范围（0，)上。调用参数B和A与系统函数的关系由下式给出： 缺省w和M时，freqz自动选取512个频率点计算。不带输出向量的freqz函数将自动绘出幅频和相频曲线。其它几种调用格式可用命令help查阅。Freqz函数用于计算模拟滤波器的频率响应函数，详细功能及调用格式用help查看。例如，八阶梳状滤波器系统函数为 H(z) =B(z)=1-z -8 用下面的简单程

5、序绘出 H(z)的幅频和相频特性曲线如图10.2.1 所示。 %Example for Freqz B=1 0 0 0 0 0 0 0 -1; A=1; freqz(B,A)五、实验报告要求1对各实验内容进行理论分析和推导。2分析各实验内容，并和理论分析推导结果进行对比 。13实验二 快速变换及其应用一、实验目的(1)在理论学习的基础上。通过本实验，加深对FFT的理解，熟悉FFT子程序；(2)熟悉应用FFT对典型信号进行频谱分析的方法；(3)了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT；(4)熟悉应用FFT实现两个序列的线性卷积的方法；(5)对DCT变换用作数

6、据压缩有初步的认识。二、实验内容实验中用到的信号序列a)余弦序列 b)衰减正弦序列c)三角波序列d)反三角波序列上机实验内容：(1) xn=1, 2, 3, 4, 5, hn=-3, -2, -1, 0, 1, 2, 3, 4, 5, 6. 借助FFT求它们的线性卷积。(2)观察余弦序列的时域和幅频特性，了解当长度N取不同值时，对信号序列的时域和幅频特性的影响。（i）信号频率f=500Hz。采样频率fs=1500 Hz。取（a）N=16 时；（b）N=160时；（c）N=1600时。观察出现的泄露现象。混叠是否同时出现？绘出泄露现象最严重时时域序列和幅频特性曲线。（ii）信号频率f=500Hz

7、。N=1600。取采样频率（a）fs=500Hz时；（b）fs=1000Hz时；（c）fs=2000Hz时。观察出现的混叠现象。绘出泄露现象最严重时时域序列和幅频特性曲线。 (3)观察衰减正弦序列Xb(n)的时域和幅频特性，a=0.1,f=0.0625,检查谱峰出现位置是否正确，注意频谱的形状，绘出时域序列和幅频特性曲线。改变f，使f分别等于0.4375和0.5625,观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄露现象？说明产生现象的原因。(4)观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列Xc(n)和Xd(n)的幅频特性，观察两者的序列形状和频谱曲线有什么

8、异同？绘出两时域序列及其幅频特性曲线。在Xc(n)和Xd(n)末尾补零，用N=16点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？两情况的FFT频谱还有相同之处吗？这些变化说明了什么？(5)一个连续信号含两个频率分量，经采样得已知N=16，f分别为1/16，1/64，观察其频谱；当N=128时，f不变，其结果有何不同，为什么？四、实验报告要求(1)简述实验目的及原理。(2)分析得到的波形，说明参数改变对时域和频域的影响，按实验要求附上实验信号序列和幅频特性曲线，(3)总结实验中的主要结论实验三 IIR数字滤波器的设计一、实验目的（1）掌握双线性变换法及脉冲响应不变法设计IIR数字

9、滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应设计低通、高通和带通IIR数字滤波器的计算机编程。（2）观察双线性变换及脉冲响应不变法设计的滤波器的频率特性，了解双线性变换法及脉冲响应不变法的特点。（3）熟悉Butterworth滤波器、Chebyshew滤波器和椭圆滤波器的频率特性。二、实验内容实验中有关变量的定义：fp 通带边界频率fr 阻带边界频率Ap 通带波动Ar 最小阻带衰减fs 采样频率T 采样周期（1）fp=0.3kHz,Ap=0.8dB,fr=0.2kHz,Ar=20dB,T=1ms。设计一Chebyshew高通滤波器；局部放大图形，观察其通带损耗和阻带衰减是否满足要

10、求。记录fp=0.3kHz和fr=0.2kHz处的通带波动和阻带衰减，绘制所设计滤波器的图形。（2）fp=0.2kHz,Ap=1dB,fr=0.3kHz,Ar=25dB,T=1ms。用双线性变换法设计一Butterworth数字低通滤波器，观察其通带损耗和阻带衰减是否满足要求。记录fp=0.3kHz和fr=0.2kHz处的通带波动和阻带衰减，绘制所设计滤波器的图形。（3）利用双线性变换法设计一Butterworth型数字带通滤波器，已知fs=30kHz,其等效的模拟滤波器指标为Ap<3dB,2kHz<f3kHz; Ar 5 dB,f 6kHz; Ar 20dB,f 1.5kHz。观

11、察其通带损耗和阻带衰减是否满足要求。记录通带波动和阻带衰减，绘制所设计滤波器的图形。三、思考题（1）双线性变换法中和之间的关系是非线性的，在实验中你注意到这样非线性关系了吗？从哪几种数字滤波器的幅频特性曲线中可以观察到这几种非线性关系？四、实验报告要求（1）简述实验目的及原理。（2）按实验要求附上所设计的滤波器传输函数及相应的幅频特性曲线，定性分析他们的性能，判断设计是否满足要求。（3）简要回答思考题。实验四 FIR数字滤波器的设计一、实验目的（1）掌握用窗函数法设计FIR滤波器的原理及方法，熟悉相应的计算机编程。（2）熟悉线性相位FIR滤波器的幅频特性和相频特性。二、实验内容（1）给定滤波器

12、设计指标通带边界频率 0.25阻带边界频率 0.35通带波动 0.03db最小阻带衰减 50db 海明窗选择固定窗自行设计MATLAB程序，得到满足设计指标的滤波器。（2）给定滤波器设计指标通带边界频率 0.5阻带边界频率 0.7通带波动 0.01db最小阻带衰减 80db用Kaiser窗自行设计MATLAB程序，得到满足设计指标的滤波器。（3）绘制所设计滤波器的幅度和相位图形三、实验报告要求（1）简述实验目的及原理。（2）按实验要求附上所设计的滤波器的MATLAB程序。（3）绘出相应的幅频特性曲线（幅度和相位）。用海明窗设计一个FIR滤波器，其中Wp=0.2*pi，Ws=0.3*pi,通带衰

13、减不大于0.25dB,阻带衰减不小于50dB。代码：wp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;M=ceil(6.6*pi/tr_width)+1;n=0:M-1;wc=(ws+wp)/2;%ideal LPF cutoff frequencyhd = ideal_lp(wc,M);w_ham=(hamming(M)'h=hd.*w_ham;% .*db,mag,pha,grd,w=freqz_m(h,1);delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1:1:501) %Min stopband attenuati

14、on%Plotssubplot(2,2,1);stem(n,hd);title('Ideal Impulse Response')axis(0 M-1 -0.1 0.3);xlabel('n');ylabel('hd(n)');%subplot(2,2,2);stem(n,w_ham);title('Hamming Window')axis(0 M-1 0 1.1);xlabel('n');ylabel('w(n)');%subplot(2,2,3);stem(n,h);title('Act

15、ual Impulse Response')axis(0 M-1 -0.1 0.3);xlabel('n');ylabel('h(n)');%实验一离散时间与系统的频谱分析%1. y(n)=x(n)+ay(n-1)%1)a=0.95B=1A=1, -0.95AM,w=freqz(B,A,'whole');am=abs(AM);an=angle(AM);subplot(2,1,1)plot(w/pi,am)xlabel('omega/pi')ylabel('|H|')grid onsubplot(2,1,2)

16、plot(w/pi,an)grid onxlabel('omega/pi')ylabel('|H|')%a=-0.5B=1；A=1, 0.5；AM,w=freqz(B,A,'whole');am=abs(AM);an=angle(AM);subplot(2,1,1)plot(w/pi,am)xlabel('omega/pi')ylabel('|H|')grid onsubplot(2,1,2)plot(w/pi,an)grid onxlabel('omega/pi')ylabel('|H|&

17、#39;)%2.y1(n)=x(n)+x(n-1), y2(n)=x(n)-x(n-1)%1)y1(n)=x(n)+x(n-1)B=1 1；A=1；AM,w=freqz(B,A,'whole');am=abs(AM);an=angle(AM);subplot(2,1,1)plot(w/pi,am)xlabel('omega/pi')ylabel('|H|')grid onsubplot(2,1,2)plot(w/pi,an)grid onxlabel('omega/pi')ylabel('|H|')%y2(n)=x

18、(n)-x(n-1)B=1 -1；A=1；AM,w=freqz(B,A,'whole');am=abs(AM);an=angle(AM);subplot(2,1,1)plot(w/pi,am)xlabel('omega/pi')ylabel('|H|')grid onsubplot(2,1,2)plot(w/pi,an)grid onxlabel('omega/pi')ylabel('|H|')%x(n)=R3(N)B=1 1 1 1；A=1；AM,w=freqz(B,A,'whole');am=a

19、bs(AM);an=angle(AM);subplot(2,1,1)plot(w/pi,am)xlabel('omega/pi')ylabel('|H|')grid onsubplot(2,1,2)plot(w/pi,an)grid onxlabel('omega/pi')ylabel('|H|') 实验二 快速变换及其应用%（1）xn=1:5hn=-3:6xnd=length(xn)hnd=length(hn)ynd=xnd+hnd-1xnL=xn zeros(1,ynd-xnd)hnL=hn zeros(1,ynd-hnd)X

20、L=fft(xnL)HL=fft(hnL)YL=XL.*HLyn=ifft(YL)%conv(xn,hn);直接计算%(2)n=0:1600;omega=2*pi*500/1500x=sin(omega*n)y=fft(x);subplot(3,1,1)stem(x);title('序列')subplot(3,1,2)plot(n,abs(y);title('幅度特性')subplot(3,1,3)plot(n,angle(y);title('相位特性')%(3)n=0:15;a=input('a=');f=input('

21、f=');x=exp(-(a*n).*sin(2*pi*f*n);y=fft(x);subplot(3,1,1)stem(x);title('序列')subplot(3,1,2)plot(n,abs(y);title('幅度序列')subplot(3,1,3)plot(n,angle(y);title('相位特性')%(4)%1) N=8; n=0:7; x=0,1,2,3,4,3,2,1;y=fft(x);subplot(3,1,1)stem(n,x);title('序列')subplot(3,1,2)plot(n,ab

22、s(y);title('幅度特性')subplot(3,1,3)plot(n,angle(y);title('相位特性')%2)N=8;n=0:7;x=4,3,2,1,0,1,2,3;y=fft(x);subplot(3,1,1)stem(n,x)subplot(3,1,2)plot(n,abs(y);title('幅度特性')subplot(3,1,3)plot(n,angle(y);title('相位特性')%3) N=16; n=0:15; x=0,1,2,3,4,3,2,1,0,0,0,0,0,0,0,0;y=fft(x)

23、;subplot(3,1,1)stem(n,x);title('序列')subplot(3,1,2)plot(n,abs(y);title('幅度特性')subplot(3,1,3)plot(n,angle(y);title('相位特性')%4) N=16; n=0:15; x=4,3,2,1,0,1,2,3,0,0,0,0,0,0,0,0;y=fft(x);subplot(3,1,1)stem(n,x);title('序列')subplot(3,1,2)plot(n,abs(y);title('幅度特性')sub

24、plot(3,1,3)plot(n,angle(y);%(4)N=input('N=');n=0:N-1;df=input('df=');x=sin(2*pi*0.125*n)+cos(2*pi*(0.125+df)*n);y=fft(x);subplot(3,1,1)stem(x);title('序列')subplot(3,1,2)plot(n,abs(y);title('幅度特性')subplot(3,1,3)plot(n,angle(y);title('相位特性')实验三 IIR滤波器设计%(1)fp=0.3

25、k,Ap=0.8db;fr=0.2k,Ar=20db;T=1%ms.T=1wp=2*0.3*T;wr=2*0.2*T;Ap=0.8;Ar=20;N,wc=cheb1ord(wp,wr,Ap,Ar) b,a=cheby1(N,Ap,wc,'high') %滤波器差分方程：%y(n)=b0x(n)+ b1x(n-1) + b2x(n-2)+%+a1y(n-1)+ a2y(n-2)+sys1=tf(b,a,-1)% w=0:pi/512:pi;H, w=freqz(b,a,w);fg=20*log10(abs(H);plot(w/pi,fg); title('幅度特性'

26、;)grid%(2)fp=0.2k,Ap=1db,fr=0.3k,Ar=25db,T=1m%s%双线性变换法wp=2*0.2*1;wr=2*0.3*1;T=0.001Ap=1;Ar=25;N,wc=buttord(wp,wr,Ap,Ar,'s');z,p,k=buttap(N)WC=(2/T)*tan(wc*pi/2)num=kn=N:-1:0den=real(poly(p)den=den./(WC.n);sys=tf(num,den)b,a=bilinear(num,den,1/T)%滤波器差分方程：%y(n)=b0x(n)+ b1x(n-1) + b2x(n-2)+%+a1y(n-1)+ a2y(n-2)+% m=linspace(0,pi,512);H,w=freqz(b,a,m);H2=20*log10(abs(H);W2=angle(H);subplot(2,1,1)hold onplot(w/pi,H2,'r');grid onsubplot(2,1,2)hold onplot(w/pi,W2,'r');grid on%(3)banpass:fs=30k;Ap<3db,2k<=f<=3k;Ar>=5db,f>=6k;Ar>=20db,f<=1.5k.wp1=2*2