数字信号处理第三版高西泉上机实验完整版(含子程序)

实验一：系统响应及系统稳定性一、实验原理与方法1、在时域求系统响应的方法有两种，第一种是通过解差分方程求得系统输出；第二种是已知系统的单位脉冲响应，通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。2、检验系统的稳定性，其方法是在输入端加入单位阶跃序列，观察输出波形，如果波形稳定在一个常数值上，系统稳定，否则不稳定。二、实验结果及程序清单（含分析及函数用法）1、%-----（2）调用filter解差分方程以及单位脉冲响应------closeall;clearallA=[1,-0.9];B=[0.05,0.05];%系统差分方程系数向量B和Ax1n=[ones(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)subplot(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,a2……am]A=[b1,b2……bn]2、%-----（3）调用conv函数计算卷积-------x1n=[ones(1,8)];%产生信号x1(n)=R8(n)h1n=[ones(1,10)zeros(1,10)];h2n=[12.52.51zeros(1,10)];y21n=conv(h1n,x1n);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函数用于计算两个有限长序列的卷积C=conv（A,B）计算两个有限长序列向量A和B的卷积3、%-----（4）实验方法检查系统是否稳定-------closeall;clearallun=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(y1n);title('(h)谐振器对u(n)的响应y31(n)');subplot(2,1,2);stem(y2n);title('(i)谐振器对正弦信号的响应y32(n)');分析：（h）中由：在系统的输入端加入单位阶跃序列，如果系统的输出趋近一个常数（包括零），就可以断定系统是稳定的输出显然趋近于零，所以是稳定的（i）中谐振器具有对某个频率进行谐振的性质，本实验中的谐振器的谐振频率是0.4rad,因此稳定波形为sin(0.4n)。三、思考题简要分析(1)如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应?如何求？答：如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可用分段线性卷积法求系统的响应（2）如果信号经过低通滤波器，把信号的高频分量滤掉，时域信号会有何变化，用前面第一个实验结果进行分析说明。答：如果信号经过低通滤波器，则信号的高频分量被滤掉，时域信号的变化减缓，在有阶跃处附近产生过渡带。因此，当输入矩形序列时，输出序列的开始和终了都产生了明显的过渡带，见第一个实验结果的波形。四、总结心里总结：起初没有看书，什么都不懂，做实验的时候也是什么都不懂，混时间，非常的痛苦。说要认真的写实验报告，也是笑话，还要分析，更是痛苦。这次发时间把书本的第一章好好的看了几遍，包括原理性的东西，以及其中涉及的matlab的代码全看了。再则就是参考了matlab基础书籍以及网上的一些资料，最重要的是matlab的帮助文档作用很大，函数怎么用的都写好了。把实验一完成后，觉得也没有什么呀，当初那么的害怕做什么呢！！！！总结即在实验原理中说明的两点：1、在时域求系统响应的方法有两种，第一种是通过解差分方程求得系统输出；第二种是已知系统的单位脉冲响应，通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。2、检验系统的稳定性，其方法是在输入端加入单位阶跃序列，观察输出波形，如果波形稳定在一个常数值上，系统稳定，否则不稳定。实验二：时域采样与频域采样一、实验原理与方法1、时域采样定理：a)对模拟信号)(txa以间隔T进行时域等间隔理想采样，形成的采样信号的频谱)(ˆjX是原模拟信号频谱()aXj以采样角频率s（Ts/2）为周期进行周期延拓。公式为：)](ˆ[)(ˆtxFTjXaa)(1nsajnjXTb)采样频率s必须大于等于模拟信号最高频率的两倍以上，才能使采样信号的频谱不产生频谱混叠。2、频域采样定理：公式为：()IDFT[()][()]()NNNNixnXkxniNRn由公式可知，频域采样点数N必须大于等于时域离散信号的长度M(即N≥M)，才能使时域不产生混叠，则N点IDFT[()NXk]得到的序列()Nxn就是原序列x(n),即()Nxn=x(n)。二、实验结果及程序清单（含分析及函数用法）1、%---------（1）时域采样理论验证程序---------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*2^0.5;omega=pi*50*2^0.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M);%M点FFT[xnt)]yn='xa(nT)';subplot(3,2,1);stem(xnt);%调用自编绘图函数stem绘制序列图boxon;title('(a)Fs=1000Hz');k=0:M-1;fk=k/Tp;subplot(3,2,2);stem(fk,abs(Xk));title('(a)T*FT[xa(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*2^0.5;omega=pi*50*2^0.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M);%M点FFT[xnt)]yn='xa(nT)';subplot(3,2,3);stem(xnt);%调用自编绘图函数stem绘制序列图boxon;title('(b)Fs=300Hz');k=0:M-1;fk=k/Tp;subplot(3,2,4);stem(fk,abs(Xk));title('(b)T*FT[xa(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.128;alph=pi*50*2^0.5;omega=pi*50*2^0.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M);%M点FFT[xnt)]yn='xa(nT)';subplot(3,2,5);stem(xnt);%调用自编绘图函数stem绘制序列图boxon;title('(c)Fs=200Hz');k=0:M-1;fk=k/Tp;subplot(3,2,6);stem(fk,abs(Xk));title('(c)T*FT[xa(nT)],Fs=200Hz');xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))]);分析：，当采样频率为1000Hz时，频谱混叠很小：当采样频率为300Hz时，频谱混叠很严重；当采样频率为200Hz时，频谱混叠更很严重。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;xn=[xa,xb];Xk=fft(xn,1024);%1024点FFT[x(n)],用于近似序列x(n)的TFX32k=fft(xn,32);%32点FFT[x(n)]x32n=ifft(X32k);%32点IFFT[X32(k)]得到x32(n)X16k=X32k(1:2:N);%隔点抽取X32k得到X16(K)x16n=ifft(X16k,N/2);%16点IFFT[X16(k)]得到x16(n)subplot(3,2,2);stem(n,xn);boxontitle('(b)三角波序列x(n)');xlabel('n');ylabel('x(n)');axis([0,32,0,20])k=0:1023;wk=2*k/1024;%subplot(3,2,1);plot(wk,abs(Xk));title('(a)FT[x(n)]');xlabel('\omega/\pi');ylabel('|X(e^j^\omega)|');axis([0,1,0,200])k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k));boxontitle('(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);boxontitle('(d)16点IDFT[X_1_6(k)]');xlabel('n');ylabel('x_1_6(n)');axis([0,32,0,20])k=0:N-1;subplo