搜索
《数字信号处理》第三版，高西全、丁玉美，实验程序清单：实验0实验内容：1、实验用Matlab工具箱函数简介熟悉以下函数的功能、格式（1）abs（2）angle（3）conv（4）filter（5）freqz（6）impz（7）fft（8）ifft（9）plot（10）stem（11）subplot（12）figure2、常用序列的产生及其频谱分析（1）单位脉冲序列及其频谱n=0:50;%定义序列的长度是50x=[1,zeros(1,50)];%注意：MATLAB中数组下标从1开始%x(1)=1;closeall;subplot(3,1,1);stem(x);title('单位冲击信号序列');k=-25:25;X=x*(exp(-j*pi/25)).^(n'*k);magX=abs(X);%绘制x(n)的幅度谱subplot(3,1,2);stem(magX);title('单位冲击信号的幅度谱');angX=angle(X);%绘制x(n)的相位谱subplot(3,1,3);stem(angX);title('单位冲击信号的相位谱')（2）矩形序列及其频谱n=1:50x=sign(sign(10-n)+1);closeall;subplot(3,1,1);stem(x);title(‘矩形序列’);k=-25:25;X=x*(exp(-j*pi/25)).^(n’*k);magX=abs(X);%绘制x(n)的幅度谱subplot(3,1,2);stem(magX);title(‘矩形序列的幅度谱’);angX=angle(X);%绘制x(n)的相位谱subplot(3,1,3);stem(angX);title(‘矩形序列的相位谱’)(3)特定冲击串及其频谱%h2(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)n=1:50;%定义序列的长度是50x=zeros(1,50);%注意：MATLAB中数组下标从1开始x(1)=1;x(2)=2.5;x(3)=2.5;x(4)=1;closeall;subplot(3,1,1);stem(x);title(‘特定冲击串’);k=-25:25;X=x*(exp(-j*pi/12.5)).^(n’*k);magX=abs(X);%绘制x(n)的幅度谱subplot(3,1,2);stem(magX);title(‘特定冲击串的幅度谱’);angX=angle(X);%绘制x(n)的相位谱subplot(3,1,3);stem(angX);title(‘特定冲击串的相位谱’)（4）采样信号序列及其频谱%x(n)=Ae-anTsin(w0nT)u(nT)（0=n=50）;%A=444.128;a=50*sqrt(2.0)*pi;fs=1kHz(T=0.001);w0=50*sqrt(2.0)*pin=0:50;%定义序列的长度是50A=444.128;%设置信号有关的参数a=50*sqrt(2.0)*pi;T=0.001;%采样率，fs=1kHzw0=50*sqrt(2.0)*pi;x=A*exp(-a*n*T).*sin(w0*n*T);%pi是MATLAB定义的π，信号乘可采用“.*”closeall%清除已经绘制的x(n)图形subplot(3,1,1);stem(x);%绘制x(n)的图形title(‘理想采样信号序列’);k=-25:25;W=(pi/12.5)*k;X=x*(exp(-j*pi/12.5)).^(n’*k);magX=abs(X);%绘制x(n)的幅度谱subplot(3,1,2);stem(magX);title(‘理想采样信号序列的幅度谱’);angX=angle(X);%绘制x(n)的相位谱subplot(3,1,3);stem(angX);title(‘理想采样信号序列的相位谱’)010203040506000.51单位冲击信号序列010203040506000.51单位冲击信号的幅度谱0102030405060-101单位冲击信号的相位谱0510152025303540455000.51矩形序列01020304050600510矩形序列的幅度谱0102030405060-505矩形序列的相位谱05101520253035404550024特定冲击串01020304050600510特定冲击串的幅度谱0102030405060-505特定冲击串的相位谱0102030405060-2000200理想采样信号序列010203040506005001000理想采样信号序列的幅度谱0102030405060-505理想采样信号序列的相位谱实验一系统响应及系统稳定性%====内容1：调用filter解差分方程，由系统对u(n)的响应判断稳定性====closeall;clearallA=[1,-0.9];B=[0.05,0.05];%系统差分方程系数向量B和Ax1n=[11111111zeros(1,50)];%产生信号x1(n)=R8(n)x2n=ones(1,128);%产生信号x2(n)=u(n)hn=impz(B,A,58);%求系统单位脉冲响应h(n)subplot(3,1,1);stem(hn);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);%y='y2(n)';stem(y2n);title('(c)系统对u(n)的响应y2(n)');boxon%===内容2：调用conv函数计算卷积=====x1n=[11111111];%产生信号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);figure(2)subplot(2,2,1);stem(h1n);%调用函数stem绘图title('(d)系统单位脉冲响应h1(n)');boxonsubplot(2,2,2);stem(y21n);title('(e)h1(n)与R8(n)的卷积y21(n)');boxonsubplot(2,2,3);stem(h2n);title('(f)系统单位脉冲响应h2(n)');subplot(2,2,4);stem(y22n);title('(g)h2(n)与R8(n)的卷积y22(n)');%====内容3：谐振器分析=======un=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和Ay31n=filter(B,A,un);%谐振器对u(n)的响应y31(n)y32n=filter(B,A,xsin);%谐振器对u(n)的响应y31(n)figure(3)subplot(2,1,1);stem(y31n);title('(h)谐振器对u(n)的响应y31(n)');subplot(2,1,2);stem(y32n);title('(i)谐振器对正弦信号的响应y32(n)');010203040506000.050.1(a)系统单位脉冲响应h(n)010203040506000.51(b)系统对R8(n)的响应y1(n)02040608010012014000.51(c)系统对u(n)的响应y2(n)0510152000.51(d)系统单位脉冲响应h1(n)010203002468(e)h1(n)与R8(n)的卷积y21(n)0510150123(f)系统单位脉冲响应h2(n)010203002468(g)h2(n)与R8(n)的卷积y22(n)050100150200250300-0.0500.05(h)谐振器对u(n)的响应y31(n)050100150200250300-1-0.500.51(i)谐振器对正弦信号的响应y32(n)实验二时域采样与频域采样1时域采样理论的验证程序清单%时域采样理论验证程序exp2a.mTp=64/1000;%观察时间Tp=64微秒%产生M长采样序列x(n)%Fs=1000;T=1/Fs;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=fft(xnt,M);%M点FFT[xnt)]yn='xa(nT)';subplot(3,2,1);stem(xnt);title('(a)Fs=1000Hz');k=0:M-1;fk=k/Tp;subplot(3,2,2);plot(fk,abs(Xk));title('(b)FT[xa(nT)],Fs=1000Hz');xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))])%================================%Fs=300Hz和Fs=200Hz的程序与上面Fs=1000Hz完全相同。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=fft(xnt,M);%M点FFT[xnt)]yn='xa(nT)';subplot(3,2,3);stem(xnt);title('(c)Fs=300Hz');k=0:M-1;fk=k/Tp;subplot(3,2,4);plot(fk,abs(Xk));title('(d)T*FT[xa(nT)],Fs=300Hz');xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))])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);title('(e)Fs=200Hz');k=0:M-1;fk=k/Tp;subplot(3,2,6);plot(fk,abs(Xk));title('(f)T*FT[xa(nT)],Fs=200Hz');xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))])2频域采样理论的验证程序清单%频域采样理论验证程序exp2b.mM=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)的FTX32k=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