数字信号处理第三版高西泉上机实验完整版含子程序副本.docx

1、数字信号处理第三版高西泉上机实验完整版含子程序副本实验一：系统响应及系统稳定性一、实验原理与方法1、在时域求系统响应的方法有两种，第一种是通过解差分方程求得系统输出； 第二种是已知系统的单位脉冲响应，通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。2、检验系统的稳定性，其方法是在输入端加入单位阶跃序列，观察输出波形，如果波形稳定在一个常数值上，系统稳定，否则不稳定。二、实验结果及程序清单（含分析及函数用法）1、%-（2）调用filter解差分方程以及单位脉冲响应-close all;clear allA=1,-0.9;B=0.05,0.05; %系统差分方程系数向量B和Ax1n=one

2、s(1,8) zeros(1,25); %产生信号x1(n)=R8(n)，用zeros用来加点的个数x2n=ones(1,30); %产生信号x2(n)=u(n)hn=impz(B,A,25); %求系统单位脉冲响应h(n)subplot(3,1,1);stem(hn); %调用函数stem绘图title(a) 系统单位脉冲响应h(n);y1n=filter(B,A,x1n); %求系统对x1(n)的响应y1(n)subplot(3,1,2);stem(y1n);title(b) 系统对R8(n)的响应y1(n);y2n=filter(B,A,x2n); %求系统对x2(n)的响应y2(n)s

3、ubplot(3,1,3);stem(y2n);title(c) 系统对u(n)的响应y2(n);分析： （a）25个点数和程序所写一致。Filter函数实现线性常系数差分方程的递推求解。调用格式如下： Y=filter(B,A,x) 计算系统对输入信号x的零状态响应输出信号向量Y，B,A是差分方程的系数向量。即 B=a1,a2am A=b1,b2bn 2、%-（3）调用conv函数计算卷积-x1n=ones(1,8); %产生信号x1(n)=R8(n)h1n=ones(1,10) zeros(1,10);h2n=1 2.5 2.5 1 zeros(1,10);y21n=conv(h1n,x1

4、n);y22n=conv(h2n,x1n);subplot(2,2,1);stem(h1n);title(d) 系统单位脉冲响应h1(n);subplot(2,2,3);stem(y21n);title(e) h1(n)与R8(n)的卷积y21(n);subplot(2,2,2);stem(h2n);title(f) 系统单位脉冲响应h2(n);subplot(2,2,4);stem(y22n);title(g) h2(n)与R8(n)的卷积y22(n);分析： （d）（f）单位脉冲响应点数与程序要求一致 （e）（g）卷积点数满足M+N-1的要求，图形也满足要求。 Conv函数用于计算两个有限

5、长序列的卷积 C=conv（A,B）计算两个有限长序列向量A和B的卷积3、%-（4）实验方法检查系统是否稳定-close all;clear allun=ones(1,256); %产生信号u(n)n=0:255;xsin=sin(0.014*n)+sin(0.4*n); %产生正弦信号A=1,-1.8237,0.9801;B=1/100.49,0,-1/100.49; %系统差分方程系数向量B和Ay1n=filter(B,A,un); %谐振器对u(n)的响应y31(n)y2n=filter(B,A,xsin); %谐振器对u(n)的响应y31(n)subplot(2,1,1);stem(y

6、1n);title(h) 谐振器对u(n)的响应y31(n);subplot(2,1,2);stem(y2n);title(i) 谐振器对正弦信号的响应y32(n);分析： （h）中由：在系统的输入端加入单位阶跃序列，如果系统的输出趋近一个常数（包括零），就可以断定系统是稳定的 输出显然趋近于零，所以是稳定的 （i）中谐振器具有对某个频率进行谐振的性质，本实验中的谐振器的谐振频率是0.4 rad,因此稳定波形为sin(0.4n)。三、思考题简要分析(1) 如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应? 如何求？答： 如果输入信号为无限长序列，系统的单位

7、脉冲响应是有限长序列，可用分段线性卷积法求系统的响应 （2）如果信号经过低通滤波器，把信号的高频分量滤掉，时域信号会有何变化，用前面 第一个实验结果进行分析说明。 答：如果信号经过低通滤波器， 则信号的高频分量被滤掉， 时域信号的变化减缓，在有阶跃处附近产生过渡带。因此，当输入矩形序列时，输出序列的开始和终了都产生了明显的过渡带，见第一个实验结果的波形。四、总结总结即在实验原理中说明的两点：1、在时域求系统响应的方法有两种，第一种是通过解差分方程求得系统输出； 第二种是已知系统的单位脉冲响应，通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。2、检验系统的稳定性，其方法是在输入端加入单位

8、阶跃序列，观察输出波形，如果波形稳定在一个常数值上，系统稳定，否则不稳定。实验二：时域采样与频域采样一、实验原理与方法1、时域采样定理：a)对模拟信号以间隔T进行时域等间隔理想采样，形成的采样信号的频谱是原模拟信号频谱以采样角频率（）为周期进行周期延拓。公式为： b)采样频率必须大于等于模拟信号最高频率的两倍以上，才能使采样信号的 频谱不产生频谱混叠。2、频域采样定理：公式为： 由公式可知，频域采样点数N必须大于等于时域离散信号的长度M(即NM)，才能使时域不产生混叠，则N点IDFT得到的序列就是原序列x(n),即=x(n)。二、实验结果及程序清单（含分析及函数用法）1、% -（1）时域采样理

9、论验证程序-Tp=64/1000; %观察时间Tp=64毫秒%产生M长采样序列x(n)Fs=1000; T=1/Fs;M=Tp*Fs; n=0:M-1;A=444.128; alph=pi*50*20.5 ;omega=pi*50*20.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M); %M点FFTxnt)yn=xa(nT);subplot(3,2,1);stem(xnt); %调用自编绘图函数stem绘制序列图box on;title(a) Fs=1000Hz);k=0:M-1;fk=k/Tp;subplot(3,2,2);ste

10、m(fk,abs(Xk);title(a) T*FTxa(nT),Fs=1000Hz);xlabel(f(Hz);ylabel(幅度);axis(0,Fs,0,1.2*max(abs(Xk);% Fs=300Hz和 Fs=200Hz的程序与上面Fs=1000Hz完全相同。Tp=64/1000; %观察时间Tp=64毫秒%产生M长采样序列x(n)Fs=300; T=1/Fs;M=Tp*Fs; n=0:M-1;A=444.128; alph=pi*50*20.5 ;omega=pi*50*20.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,

11、M); %M点FFTxnt)yn=xa(nT);subplot(3,2,3);stem(xnt); %调用自编绘图函数stem绘制序列图box on;title(b) Fs=300Hz);k=0:M-1;fk=k/Tp;subplot(3,2,4);stem(fk,abs(Xk);title(b) T*FTxa(nT),Fs=300Hz);xlabel(f(Hz);ylabel(幅度);axis(0,Fs,0,1.2*max(abs(Xk);Tp=64/1000; %观察时间Tp=64毫秒%产生M长采样序列x(n)Fs=200; T=1/Fs;M=Tp*Fs; n=0:M-1;A=444.12

12、8; alph=pi*50*20.5 ;omega=pi*50*20.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M); %M点FFTxnt)yn=xa(nT);subplot(3,2,5);stem(xnt); %调用自编绘图函数stem绘制序列图box on;title(c) Fs=200Hz);k=0:M-1;fk=k/Tp;subplot(3,2,6);stem(fk,abs(Xk);title(c) T*FTxa(nT),Fs=200Hz);xlabel(f(Hz);ylabel(幅度);axis(0,Fs,0,1.2*ma

13、x(abs(Xk);分析：，当采样频率为1000 Hz时，频谱混叠很小：当采样频率为300 Hz时，频谱混叠很严重；当采样频率为200 Hz时，频谱混叠更很严重。 Fft函数的调用格式： Xk=fft(xn,N)调用参数xn为被交换的时域序列向量，N是DFT变换的区间长度，当N大于xn的长度时，fft函数自动在xn后面补零。当N小于xn的长度时，fft函数计算xn的前面N个元素构成的N长序列的N点DFT,忽略xn后面的元素。2、%-(2)频域采样理论验证-M=27;N=32;n=0:M;%产生M长三角波序列x(n)xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; x

14、n=xa,xb;Xk=fft(xn,1024); %1024点FFTx(n), 用于近似序列x(n)的TFX32k=fft(xn,32) ;%32点FFTx(n)x32n=ifft(X32k); %32点IFFTX32(k)得到x32(n)X16k=X32k(1:2:N); %隔点抽取X32k得到X16(K)x16n=ifft(X16k,N/2); %16点IFFTX16(k)得到x16(n)subplot(3,2,2);stem(n,xn);box ontitle(b) 三角波序列x(n);xlabel(n);ylabel(x(n);axis(0,32,0,20)k=0:1023;wk=2*

15、k/1024; %subplot(3,2,1);plot(wk,abs(Xk);title(a)FTx(n);xlabel(omega/pi);ylabel(|X(ejomega)|);axis(0,1,0,200)k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k);box ontitle(c) 16点频域采样);xlabel(k);ylabel(|X_1_6(k)|);axis(0,8,0,200)n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n);box ontitle(d)16点IDFTX_1_6(k);xlabel(n);yl

16、abel(x_1_6(n);axis(0,32,0,20)k=0:N-1;subplot(3,2,5);stem(k,abs(X32k);box ontitle(e) 32点频域采样);xlabel(k);ylabel(|X_3_2(k)|);axis(0,16,0,200)n1=0:N-1;subplot(3,2,6);stem(n1,x32n);box ontitle(f)32点IDFTX_3_2(k);xlabel(n);ylabel(x_3_2(n);axis(0,32,0,20)分析; 该图验证了频域采样理论和频域采样定理。 对信号x(n)的频谱函数X(ej)在0， 2上等间隔采样N

17、=16时， N点IDFTXN(k)得到的序列正是原序列x(n)以16为周期进行周期延拓后的主值区序列当N=16时，由于NM，频域采样定理，所以不存在时域混叠失真 Ifft函数用法同fft函数三、思考题简要分析1、如果序列x(n)的长度为M，希望得到其频谱在上的N点等间隔采样，当NM时， 如何用一次最少点数的DFT得到该频谱采样？答：先对原序列x(n)以N为周期进行周期延拓后取主值区序列，再计算N点DFT则得到N点频域采样：实验三：用FFT对信号作频谱分析一、实验原理与方法1、周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量

18、选择信号的观察时间长一些。 2、对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。 3、频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是，因此要求。二、实验结果及程序清单（含分析及函数用法）1、clear all;close all%-实验内容(1)-x1n=ones(1,4); %产生序列向量x1(n)=R4(n)M=8; xa=1:(M/2); xb=(M/2):-1:1; x2n=xa,xb; %产生长度为8的三角波序列x2(n)x3n=xb,xa;X1

19、k8=fft(x1n,8); %计算x1n的8点DFTX1k16=fft(x1n,16); %计算x1n的16点DFTX2k8=fft(x2n,8); %计算x1n的8点DFTX2k16=fft(x2n,16); %计算x1n的16点DFTX3k8=fft(x3n,8); %计算x1n的8点DFTX3k16=fft(x3n,16); %计算x1n的16点DFT%以下绘制幅频特性曲线subplot(3,2,1);mstem(X1k8); %绘制8点DFT的幅频特性图title(1a) 8点DFTx_1(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(

20、X1k8)subplot(3,2,2);mstem(X1k16); %绘制16点DFT的幅频特性图title(1b)16点DFTx_1(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X1k16)subplot(3,2,3);mstem(X2k8); %绘制8点DFT的幅频特性图title(2a) 8点DFTx_2(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X2k8)subplot(3,2,4);mstem(X2k16); %绘制16点DFT的幅频特性图title(2b)16点DFTx_2(n);

21、xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X2k16)subplot(3,2,5);mstem(X3k8); %绘制8点DFT的幅频特性图title(3a) 8点DFTx_3(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X3k8)subplot(3,2,6);mstem(X3k16); %绘制16点DFT的幅频特性图title(3b)16点DFTx_3(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X3k16)2、%-实验内容（2）-N=8;n=0:N

22、-1; %FFT的变换区间N=8x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k8=fft(x4n); %计算x4n的8点DFTX5k8=fft(x5n); %计算x5n的8点DFTN=16;n=0:N-1; %FFT的变换区间N=16x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k16=fft(x4n); %计算x4n的16点DFTX5k16=fft(x5n); %计算x5n的16点DFTsubplot(2,2,1);mstem(X4k8); %绘制8点DFT的幅频特性图title(a) 8点DFTx

23、_4(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X4k8)subplot(2,2,3);mstem(X4k16); %绘制16点DFT的幅频特性图title(b)16点DFTx_4(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X4k16)subplot(2,2,2);mstem(X5k8); %绘制8点DFT的幅频特性图title(a) 8点DFTx_5(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X5k8)subplot(2,2,4);mst

24、em(X5k16); %绘制16点DFT的幅频特性图title(b)16点DFTx_5(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X5k16)分析：的周期为8，所以N=8和N=16均是其周期的整数倍，得到正确的单一频率正弦波的频谱，仅在0.25处有1根单一谱线。如图（a）和（b）所示。的周期为16，所以N=8不是其周期的整数倍，得到的频谱不正确，如图（a）所示。N=16是其一个周期，得到正确的频谱，仅在0.25和0.125处有2根单一谱线, 如图（b）所示。3、%-实验内容（3）-Fs=64;T=1/Fs;N=16;n=0:N-1; %FFT

25、的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)16点采样X6k16=fft(x6nT); %计算x6nT的16点DFTX6k16=fftshift(X6k16); %将零频率移到频谱中心Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F; %产生16点DFT对应的采样点频率（以零频率为中心）subplot(3,1,1);stem(fk,abs(X6k16),.);box on %绘制8点DFT的幅频特性图title(6a) 16点|DFTx_6(nT)|);xlabel(f(H

26、z);ylabel(幅度);axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16)N=32;n=0:N-1; %FFT的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)32点采样X6k32=fft(x6nT); %计算x6nT的32点DFTX6k32=fftshift(X6k32); %将零频率移到频谱中心Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F; %产生16点DFT对应的采样点频率（以零频率为中心）subplot(3,1,2);stem(f

27、k,abs(X6k32),.);box on %绘制8点DFT的幅频特性图title(6b) 32点|DFTx_6(nT)|);xlabel(f(Hz);ylabel(幅度);axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32)N=64;n=0:N-1; %FFT的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)64点采样X6k64=fft(x6nT); %计算x6nT的64点DFTX6k64=fftshift(X6k64); %将零频率移到频谱中心Tp=N*T;F=1/Tp; %

28、频率分辨率Fk=-N/2:N/2-1;fk=k*F; %产生16点DFT对应的采样点频率（以零频率为中心）subplot(3,1,3);stem(fk,abs(X6k64),.); box on%绘制8点DFT的幅频特性图title(6a) 64点|DFTx_6(nT)|);xlabel(f(Hz);ylabel(幅度);axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64)分析： x6(t)有3个频率成分， f1=4 Hz, f2=8 Hz, f3=10 Hz。所以x6(t)的周期为0.5 s。 采样频率Fs=64 Hz=16f1=8f2=6.4f3。变换区间N

29、=16时，观察时间Tp=16T=0.25 s，不是x6(t)的整数倍周期， 所以所得频谱不正确，如图8.3.1（6a）所示。变换区间N=32，64 时，观察时间Tp=0.5 s，1 s，是x6(t)的整数周期，所以所得频谱正确补充：上述三个程序中用到的子程序：function mstem(Xk)M=length(Xk);k=0:M-1; wk=2*k/M;stem(wk,abs(Xk),.); box on;xlabel(/); ylabel(幅度);axis(0,2,0,1.2*max(abs(Xk);三、思考题简要分析（1）对于周期序列，如果周期不知道，如何用FFT进行谱分析？答：可先截取M点进行DFT,再将截取长度扩大1倍，比较两次的结果。如果二者的主谱差别满足分析误差要求，则以两者中的一个近似表示周期序列的频谱，否则，继续把截取长度加倍，重