DSP实验Matlab程序实例.docx

1、DSP实验Matlab程序实例2离散时间信号（序列）的产生利用MATLAB产生和绘制下列有限长序列：单位脉冲序列，单位阶跃序列，矩形序列x=-10:10;y=; for i=1:21 if x(i)=0 y(i)=1; else y(i)=0; endendsubplot(3,1,1);%在子图画单位冲激序列stem(x,y,*);title(冲激); xlabel(n); hold on; for i=1:21 if x(i)0 y(i)=0; else y(i)=1; endendsubplot(3,1,2); %在子图画单位阶跃序列stem(x,y,*); title(阶跃); xlab

2、el(n);hold on; for i=1:21 if x(i)4 y(i)=0; else y(i)=1; end end end subplot(3,1,3);%在子图画矩形脉冲序列stem(x,y,*); title(矩形); xlabel(n); hold on在一幅图上绘出曲线，A=2, =0.5, f=2Hz，A1=1, A2=0.5, A3=0.2，f=2Hz。t=0:0.5:720*pi/180;x1=5*sin(pi*t/5+pi/3);x2=5*cos(t-pi/3);x3=2*sin(3*t).*cos(2*t-pi/3);subplot(3,1,1); %在子图中画出

3、一系列正弦余弦曲线plot(t,x1,-red,t,x2,-gr,t,x3,-bl);title(正弦余弦曲线);xlabel(t);hold on;A=2;a=0.5;f=2;x4=A*exp(-a*t).*sin(2*pi*f*t);subplot(3,1,2);%在子图中画出正弦衰减信号plot(t,x4)title(正弦衰减信号);xlabel(t);hold on;A1=1;A2=0.5;A3=0.2;f=2;x5=A1*sin(2*pi*f*t)+A2*sin(2*pi*2*f*t)+A3*sin(2*pi*3*f*t);subplot(3,1,3);%在子图中画出谐波信号plot

4、(t,x5)title(谐波信号);xlabel(t);hold on;3序列的运算生成下列序列：(1) 利用MATLAB编程完成上述两序列的卷积，并绘制运算后序列的波形。x=1 2 3 4 5;h=1 2 1 2;n1=0 1 2 3 4;n2=2 3 4 5;y=conv(x,h); %计算序列x与h的卷积和序列yn0=n1(1)+n2(1); %计算卷积和序列y的起点位置n3=length(x)+length(h)-2+n0; %计算卷积和序列y的终止位置n=n0:1:n3; %确定卷积和y非零样值的时间向量subplot(3,1,1);stem(n1,x,.); %在子图绘x(n)ti

5、tle(x(n);xlabel(n);ylabel(x(n);subplot(3,1,2);stem(n2,h,.); %在子图绘h(n)title(h(n);xlabel(n);ylabel(h(n);subplot(3,1,3);stem(n,y,.); %在子图绘卷积和y(n)title(y(n)=x(n)*h(n);xlabel(n);ylabel(y(n);（2）编写一个完成卷积的函数，输入是两个序列，输出是结果。function y,n=jsjuanji(x,n1,h,n2)%计算序列卷积和y(n)=x(n)*h(n)% y: 卷积和y(n)对应的非零样值向量% n：卷积和y(n)

6、对应的时间向量% x: x(n)对应的非零样值向量% n1: x(n)对应的时间向量% h: h(n)对应的非零样值向量% n2: h(n)对应的时间向量n3=n1(1)+n2(1); %卷积和序列y的起始位置n4=length(x)+length(h)-1;%卷积和序列y的长度m=0;for i=1:n4 %求卷积和y if ilength(x) x(i)=0; end for j=1:i if jlength(h) h(j)=0; end m=m+x(i+1-j)*h(j); end y(i)=m; m=0;endn1=n1(1):1:n1(1)+n4-1; n1: % x(n)对应的新的

7、时间向量n2=n2(1):1:n2(1)+n4-1; % h(n)对应的新的时间向量n=n3:1:n3+n4-1; %y(n)对应的时间向量subplot(3,1,1);stem(n1,x,.); %在子图绘x(n)title(x(n);xlabel(n);ylabel(x(n);subplot(3,1,2);stem(n2,h,.); %在子图绘h(n)title(h(n);xlabel(n);ylabel(h(n);subplot(3,1,3);stem(n,y,.); %在子图绘卷积和y(n)title(y(n)=x(n)*h(n);xlabel(n);ylabel(y(n);4采样定理

8、的研究分别令采样周期Ts为不同值，绘出不同采样周期下x(t) =sin(t)的频谱（直接用FFT函数），观察频谱混叠现象和防止混叠的采样周期。fs=2; %设定采样频率T=1/fs; %采样周期N=128; %采样点数n=0:N-1; t=n*T; %采样时间点x=sin(t); %生成正弦信号subplot(4,2,1);plot(t,x);%画正弦信号的时域波形xlabel(t);ylabel(x);title(正弦信号x=sin(t)时域波形); y=fft(x,N); %进行fft变换mag=sqrt(y.*conj(y); %求幅值f=(0:N-1)*fs/N; %进行对应的频率转换

9、subplot(4,2,2);plot(f,mag); %画频谱图axis(-10,50,0,80);xlabel(频率);ylabel(幅值);title(正弦信号x=sin(t)幅频谱图 fs=2);hold on;fs=10; %设定采样频率T=1/fs; %采样周期N=1024; %采样点数n=0:N-1;t=n*T; %采样时间点x=sin(t); %生成正弦信号subplot(4,2,3);plot(t,x);%画正弦信号的时域波形xlabel(t);ylabel(x);title(正弦信号x=sin(t)时域波形);y=fft(x,N); %进行fft变换mag=sqrt(y.*

10、conj(y); %求幅值f=(0:N-1)*fs/N; %进行对应的频率转换subplot(4,2,4);plot(f,mag); %画频谱图axis(-10,50,0,400);xlabel(频率);ylabel(幅值);title(正弦信号x=sin(t)幅频谱图 fs=10);hold on;fs=20; %设定采样频率T=1/fs; %采样周期N=1024; %采样点数n=0:N-1;t=n*T; %采样时间点x=sin(t); %生成正弦信号subplot(4,2,5);plot(t,x);%画正弦信号的时域波形xlabel(t);ylabel(x);title(正弦信号x=sin

11、(t)时域波形);y=fft(x,N); %进行fft变换mag=sqrt(y.*conj(y); %求幅值f=(0:N-1)*fs/N; %进行对应的频率转换subplot(4,2,6);plot(f,mag); %画频谱图axis(-10,50,0,400);xlabel(频率);ylabel(幅值);title(正弦信号x=sin(t)幅频谱图 fs=20);hold on;fs=30; %设定采样频率T=1/fs; %采样周期N=1024; %采样点数n=0:N-1; t=n*T; %采样时间点x=sin(t); %生成正弦信号subplot(4,2,7);plot(t,x);%画正弦

12、信号的时域波形xlabel(t);ylabel(x);title(正弦信号x=sin(t)时域波形); y=fft(x,N); %进行fft变换mag=sqrt(y.*conj(y); %求幅值f=(0:N-1)*fs/N; %进行对应的频率转换subplot(4,2,8);plot(f,mag); %画频谱图axis(-10,50,0,400);xlabel(频率);ylabel(幅值);title(正弦信号x=sin(t)幅频谱图 fs=30);五、思考题1如何产生方波信号序列和锯齿波信号序列？答：方波信号由square函数可得，锯齿波信号由循环语句产生。%方波T=0:0.001:2*pi

13、; %方波的时间向量y=square(6*T); %方波w=6,周期为T=2*pi/6subplot(2,1,1);plot(T,y); %画方波title(方波);axis(T(1)-1 T(end)+1 -2 2); %坐标轴区域hold on%锯齿波x=-10:10;y=;for i=1:21 %rem为取余函数，此循环使得y(i)交替取值1和-1if rem(i,2)=0 y(i)=-1;else y(i)=1;endendsubplot(2,1,2);plot(x,y);title(锯齿波);1编写一个 DFT 计算程序，要求该程序具有正变换及反变换的双重功能。令 x(n)=sin(

14、n)n(-1)n, n=0,1,.,31, 这是一个震荡信号，由它来验证所编程序是否正确（即：对 x(n)作DFT 得 X(k)，再对 X(k)作反变换，检验是否得到同样的 x(n)）。function X,x2=bianhuan(x,N)%DFT X=zeros(1,N);for m=1:N for n=1:N X(m)=X(m)+x(n)*exp(-j*2*pi*n*m/N); endend %IDFT x2=zeros(1,N);for m=1:N for n=1:N x2(m)=x2(m)+X(n)*exp(j*2*pi*n*m/N); end x2(m)=x2(m)/N;endcle

15、arN=32;x=;for p=1:N x(p)=sin(p)*p*(-1)pendX,x2=bianhuan(x,N) n=0:1:N-1;subplot(3,1,1);stem(n,x,*);title(original);subplot(3,1,2);stem(n,X,*);title(after DFT);subplot(3,1,3);stem(n,x2,*);title(after IDFT);2. 序列 x(n)的 N 点 DTFT 的物理意义是对 X(ej)的在0，2上进行 N 点等间隔采样clear;N1=1024;w1=(0:2*pi/1024:2*pi);X1=(1-exp

16、(-j*4*w1)./(1-exp(-j*w1); N2=8;for i=1:N2 X2(i)=X1(i*N1/N2); w2(i)=2*pi*i/N2;end N3=16;for i=1:N3 X3(i)=X1(i*N1/N3); w3(i)=2*pi*i/N3;end subplot(3,2,1); plot(w1,abs(X1);titlesubplot(3,2,2), plot(w1,angle(X1);titlesubplot(3,2,3); stem(w2,abs(X2);titlesubplot(3,2,4), stem(w2,angle(X2);titlesubplot(3,2

17、,5); stem(w3,abs(X3);titlesubplot(3,2,6), stem(w3,angle(X3);title 4. 混叠现象研究衰减正弦信号 x(t)=1+sin(7f0t)cos(2f0t), f0=120Hz, fs=200Hz,采样点数 N=64，用 DFT 计算信号的幅频和相频，观察混叠现象clear;f0=120;N=64; fs1=200;t1=(0:1:N-1)/fs1;x1=(1+sin(7*pi*f0*t1).*cos(2*pi*f0*t1);y1=fft(x1); fs2=1500;t2=(0:1:N-1)/fs2;x2=(1+sin(7*pi*f0*

18、t2).*cos(2*pi*f0*t2);y2=fft(x2); subplot(2,1,1);stem(abs(y1); title(fs=200);subplot(2,1,2);stem(abs(y2); title(fs=1500);5. 泄漏现象 令余弦信号 x(t)= cos(2f0t),f0=60Hz,fs=200Hz,采样点数 N=64要求：先用 DFT 计算采样信号 x(n)的频谱 X(k)，显示结果图形。然后将采样信号 x(n)乘上哈明窗函数：x4(t)=0.54-0.46cos(2n/N),n=0,1,2,N-1用 DFT 计算加窗后信号的频谱 XW(k)，分别在屏幕上显示

19、信号幅频 X(k)，窗函数幅频 X4(k)及加窗后信号的频谱 XW(k)，比较频谱 X(k)与 XW(k)的不同，并给出解释。clear;f0=60;fs=200;N=64;t=(0:1:N-1)/fs;x=cos(2*pi*f0*t);X=fft(x); n=0:1:N-1;x1=0.54-0.46.*cos(2*pi*n/N);X1=fft(x1); x2=x.*x1;X2=fft(x2); subplot(3,1,1);stem(abs(X);title；subplot(3,1,2);stem(abs(X1);title;subplot(3,1,3);stem(abs(X2);title

20、;6. 栅栏效应将内容1中的信号补零至L, L 分别为 16，32，64，计算 x(n)的频谱 X(k)，并与 X(k)比较，将频谱图 X(k)与 X(k)显示出来，观察补零的效果。clear;N=8; x=;for p=1:N x(p)=sin(p)*p*(-1)pendX=zeros(1,N);for p=1:N for q=1:N X(p)=X(p)+x(q)*exp(-j*2*pi*p*q/N); endend L1=16; x1=;for p=1:L1 if p=8 x1(p)=sin(p)*p*(-1)p else x1(p)=0 endend X1=zeros(1,L1); fo

21、r p=1:L1 for q=1:L1 X1(p)=X1(p)+x1(q)*exp(-j*2*pi*p*q/L1); endend L2=32; x2=;for p=1:L2 if p=8 x2(p)=sin(p)*p*(-1)p else x2(p)=0 endend X2=zeros(1,L2); for p=1:L2 for q=1:L2 X2(p)=X2(p)+x2(q)*exp(-j*2*pi*p*q/L2); endend L3=64; x3=;for p=1:L3 if p=8 x3(p)=sin(p)*p*(-1)p else x3(p)=0 endend X3=zeros(1

22、,L3); for p=1:L3 for q=1:L3 X3(p)=X3(p)+x3(q)*exp(-j*2*pi*p*q/L3); endend n=0:1:N-1;n1=0:1:L1-1;n2=0:1:L2-1;n3=0:1:L3-1;subplot(4,1,1);stem(n,X,*);title;subplot(4,1,2);stem(n1,X1,*);title(X1:L=16);subplot(4,1,3);stem(n2,X2,*);title(X2:L=32);subplot(4,1,4);stem(n3,X3,*);title(X3:L=64);1.基于 FFT 的卷积计算对

23、上述序列做基于 FFT 的卷积计算，做 FFT 时注意延拓序列长度，将结果显示出来，并和直接卷积结果比较.clear;N1=40;N2=30;N=41+31-1; x1=cos(0.025*pi*(0:N1).2);x=x1 zeros(1,N-N1); h1=sin(2.5*pi*(0:N2);h=h1 zeros(1,N-N2); y1=conv(x1,h1); X=fft(x); H=fft(h)Y=X.*H;y=ifft(Y); subplot(2,1,1);stem(y,.);title(基于FFT卷积计算);subplot(2,1,2);stem(y1,.);title(直接卷积计

24、算);2.基于 FFT 的相关计算请计算 h(n)和 x(n)的基于 FFT 的相关函数，将结果显示出来，并和直接相关结果比较。clear;n1=0:128;n2=129:192;n3=193:256;n4=0:31;N=257+32-1; x1=cos(0.25*pi*n1.2)+sin(0.5*pi*n1);x2=cos(pi*n2)+sin(0.5*pi*n2)/6;x3=sin(0.5*pi*n3)/6+(-1).n3;x=x1 x2 x3 zeros(1,(N-257); h1=cos(pi*n4)+sin(0.5*pi*n4)/6;h=h1 zeros(1,(N-32); r1=x

25、corr(x1 x2 x3,h1); X=fft(x); H=fft(h);R=X.*conj(H);r=ifft(R); subplot(2,1,1);stem(r,.);title(FFT相关计算);subplot(2,1,2);stem(r1,.);title(直接计算);3.功率谱计算信号设为： f1=42Hz, f2=45Hz, f3=47Hz,采样频率 fs=300Hz,采样点数 N=256。计算信号的功率谱，将功率谱图显示出来。功率谱定义为 P(k)= |X(k)|2clear;f1=42;f2=45;f3=47;fs=300;N=64;t=(0:1:N-1)/fs;x=sin(

26、2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t); X=abs(fft(x);P=X.2;stem(P,.);title(功率谱)1给定差分方程y(n)-y(n-1)+0.9y(n-2)=x(n)（1）计算并画出单位脉冲响应h(n)，n=0,1,127；（2）计算并画出单位阶跃响应s(n)，n=0,1,127；（3）判断系统是否稳定。clear;num=1 0 0;den=1 -1 0.9;subplot(2,1,1);dimpulse(num,den,128);title(单位冲激响应); %绘subplot(2,1,2);dstep(num,den,128)

27、;title(单位阶跃响应);y1=dimpulse(num,den,128) %求y2=dstep(num,den,128)以下为定义法% x1=zeros(1,130); %单位冲激函数% x1(3)=1;% y1(1)=0;y1(2)=0;% for k=1:128% y1(k+2)=x1(k+2)+y1(k+1)-0.9*y1(k);% end% % x2=ones(1,130); %单位阶跃函数% y2(1)=0;y2(2)=0;% for k=1:128% y2(k+2)=x2(k+2)+y2(k+1)-0.9*y2(k);% end% % figure% subplot(2,1,

28、1);plot(y1);title(单位脉冲响应) % subplot(2,1,2);plot(y2);title(单位阶跃响应)% y1 % y22给定系统函数（1）画出H(z)零-极点示意图；（2）求出并绘出H(z)的幅频响应、相频响应；（3）求出并绘出该系统的单位抽样响应 h(n)；（4）求系统的单位阶跃响应s(n)。clear;num=-0.2 0;den=1 0 0.7;sys=tf(num,den);figure;pzmap(sys);title(零极点分布图)figure;freqz(num,den);title(频率响应)figure;subplot(2,1,1);dimpulse(num,den);title(单位冲激响应);subplot(2,1,2);dstep(num,den);title(单位阶跃响应);y1=dimpulse(num,den,128) %求y2=d