数字信号处理实验讲义

实验一 序列的产生及绘图一、实验目的1．熟悉信号处理软件MATLAB 的使用。

2．离散信号的基本运算实现。

3．了解基本序列及复杂序列的产生方法。

4．运用卷积方法观察系统的时域特性。

5．掌握线性时不变系统的频域表示方法。

二、实验内容 1．熟悉扩展函数 2．运行例题程序 3．编程实现下列内容（1）利用扩展函数产生序列并画图(a) )4()2(*2)(--+=n n n x δδ -5<=n<=5(b) )04.0cos()(n n x π=和)(2.0)04.0cos()(n w n n y +=π 0<=n<=50w(n)为白噪声 函数为 w=randn(size(n)) （2）设线性移不变系统的抽样响应为 )()9.0()(n u n h n =输入序列为 )10()()(--=n u n u n x 求系统输出 y(n)并画图 提示: 输出为输入和抽样响应的卷积三、实验报告要求1．记录例题程序的实验结果、图形。

2．写出自己编写的程序并记录结果、图形。

注：以下程序中所有以 % 开头的行均为注释， 所有汉字均为注释，%后的内容不用写入程序 %如果要了解哪个函数的应用方法请用help 命令 如help zreos%本软件中 * 表示乘法, 卷积用函数 conv 或修改后的卷积 conv_m %以下是7个扩展函数%扩展函数1~7的用法和该软件自带函数用法一致,即在调用时要将实参代入 %例：应用扩展函数3需要输入x1(n),x2(n)的值。

%在Command Window （命令窗口）中输入 % n1=1:5; % n2=2:6;% x1=[1 3 5 7 9]; % x2=[2 4 6 8 10];% [y,n]=sigadd(x1,n1,x2,n2) 即可得两序列相加的结果%7个扩展函数要分别存到不同的文件中,并且文件名要和该扩展函数的函数名一致 %如产生单位取样序列的函数所存文件的文件名必须为 impseq %1.单位取样序列 x(n)=delta(n-n0) 要求n1<=n0<=n2 function[x,n]=impseq(n0,n1,n2) n=[n1:n2]; x=[(n-n0)==0];%2.单位阶跃序列 x(n)=u(n-n0) 要求n1<=n0<=n2 function[x,n]=stepseq(n0,n1,n2) n=[n1:n2]; x=[(n-n0)>=0];%3.信号加 y(n)=x1(n)+x2(n) %find 函数：找出非零元素的索引号%x1:第一个序列的值，n1:序列x1的索引号 %x2:第二个序列的值，n2:序列x2的索引号 function[y,n]=sigadd(x1,n1,x2,n2)n=min(min(n1),min(n2)):max(max(n1),max(n2));y1=zeros(1,length(n)); y2=y1;y1(find((n>=min(n1))&(n<=max(n1))==1))=x1; y2(find((n>=min(n2))&(n<=max(n2))==1))=x2; y=y1+y2;%4.信号乘 y(n)=x1(n)*x2(n) function[y,n]=sigmult(x1,n1,x2,n2)n=min(min(n1),min(n2)):max(max(n1),max(n2)); y1=zeros(1,length(n)); y2=y1;y1(find((n>=min(n1))&(n<=max(n1))==1))=x1; y2(find((n>=min(n2))&(n<=max(n2))==1))=x2; y=y1.*y2;%5.移位 y(n)=x(n-n0)function[y,n]=sigshift(x,m,n0) n=m+n0; y=x;%6.翻褶 y(n)=x(-n)function[y,n]=sigfold(x,n) y=fliplr(x); n=-fliplr(n);%7.修改后的卷积(系统自带卷积函数(conv)只能得到输出值,而无法表示取值范围(n),修改后的卷 %积既给出了卷积值,也给出了它的取值范围.) function[y,ny]=conv_m(x,nx,h,nh) nyb=nx(1)+nh(1);nye=nx(length(x))+nh(length(h)); ny=[nyb:nye]; y=conv(x,h);%例1 序列的基本运算 nx1=1:5; nx2=2:6;x1=[1 3 5 7 9]; x2=[2 4 6 8 10];[y1,n1]=sigadd(x1,nx1,x2,nx2) %序列相加 [y2,n2]=sigmult(x1,nx1,x2,nx2) %序列相乘 n0=3;[y3,n3]=sigshift(x1,nx1,n0) %序列移位，n0为移动的位数 [y4,n4]=sigfold(x2,nx2) %序列的翻褶 y5=sum(x1) %序列和 ∑==Nn n x y 1)( y6=prod(x1) %序列积 ∏==Nn n x y 1)(ex=sum(x1.*conj(x1)) %或 ex=sum(abs(x).^2);%信号能量 ∑∑====Nn Nn n x n xn x Ex 121*|)(|)()(%例2 产生序列并画图%)]20()10([*10)]10()([*)()10(3.0---+--=--n u n u e n u n u n n x n 0<=n<=20%subplot 函数: 将一个图形分为几个小区域,每次选择一个作为当前区域画图 % subplot(m,n,p) 将图形分为m 行, n 列的区域, p 为当前区域 %axis 函数：控制坐标系的刻度和形式 title 函数：图形标题 %xlabel 函数：X 轴标记 ylabel 函数：Y 轴标记 %stem 函数：画离散序列图 n=[0:20];x1=n.*(stepseq(0,0,20)-stepseq(10,0,20));x2=10*exp(-0.3*(n-10)).*(stepseq(10,0,20)-stepseq(20,0,20)); x=x1+x2;subplot(2,1,1);stem(n,x);xlabel('n');ylabel('x(n)');axis([0,20,-1,11]);%例3 产生复信号 n j e n x )3.01.0()(+-= -10<=n<=10 %并画出复序列的实部、虚部、幅值和相位离散图 figure(2);clfn=[-10:10]; alpha=-0.1+0.3j; x=exp(alpha*n);subplot(2,2,1);stem(n,real(x));title('real');xlabel('n') subplot(2,2,2);stem(n,imag(x));title('imag');xlabel('n') subplot(2,2,3);stem(n,abs(x));title('magtitude'); xlabel('n')subplot(2,2,4);stem(n,(180/pi)*angle(x));title('phase'); xlabel('n')%例4 线性时不变系统的频域表示 b=[0.2 0.1 ];a=[1 -0.4 -0.5]; %系统函数的系数 iN i i Mi ii z a zb z H -==-∑∑+=101)(h=impz(b,a); %系统的单位取样响应 figure(1)stem(h) %画出单位取样响应 title('h(n)') figure(2)fs=1000;[H,f]=freqz(b,a,256,fs); %求出系统的频率响应 mag=abs(H); %幅度响应 ph=angle(H); %相位响应 ph=ph*180/pi;subplot(2,1,1),plot(f,mag);grid %画出幅度响应 xlabel('frequency(Hz)'); ylabel('magnitude');subplot(2,1,2);plot(f,ph);grid %画出相位响应 xlabel('frequency(Hz)'); ylabel('phase'); figure(3)zr=roots(b) %求出系统的零点 pk=roots(a) %求出系统的极点zplane(b,a); %zplane 函数画出零极点图实验二 离散傅立叶变换及谱分析一、 实验目的1．掌握离散傅里叶变换的计算机实现方法。

2．检验实序列傅里叶变换的性质。

3．掌握计算序列的圆周卷积的方法。

4．熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。

5．学习用DFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差，以便在实际中正确应用DFT 。

二、 实验内容1．实现离散傅里叶变换。

2．计算序列圆周卷积。

3．计算实序列傅里叶变换并检验DFT 性质。

4．实现连续信号傅里叶变换以及由不同采样频率采样得到的离散信号的傅里叶变换。

5．实现补零序列的傅里叶变换。

6．实现高密度谱和高分辨率谱，并比较二者的不同。

三、 实验报告要求 见各程序要求%以下为4个扩展函数% （1）离散傅立叶变换 ∑-==1)()(N k nk N W n x k X 采用矩阵相乘的方法function [Xk]=dft(xn,N) n=[0:1:N-1]; k=[0:1:N-1];WN=exp(-j*2*pi/N); nk=n'*k;WNnk=WN.^(nk); Xk=xn*WNnk;%（2）逆离散傅立叶变换 ∑-=-=10)(1)(N k nk N W k X Nn xfunction [xn]=idft(Xk,N) n=[0:1:N-1]; k=[0:1:N-1];WN=exp(-j*2*pi/N); nk=n'*k;WNnk=WN.^(-nk); xn=(Xk*WNnk)/N;% (3) 实序列的分解% 实序列可分解为共扼对称分量 ]x((-n))[x(n)*(1/2)xec N += % 和共扼反对称分量 ]x((-n))-[x(n)*(1/2)xoc N = function [xec,xoc]=circevod(x) N=length(x); n=0:(N-1);xec=0.5*(x+x(mod(-n,N)+1)); %根据上面的公式计算，mod 函数为取余 xoc=0.5*(x-x(mod(-n,N)+1));% (4) 序列的循环移位 N m n x n y ))(()(-= function y=cirshftt(x,m,N)if length(x)>Nerror('N mustbe >= the length of x') %要求移位周期大于信号长度endx=[x zeros(1,N-length(x))];n=[0:1:N-1];n=mod(n-m,N);y=x(n+1);%例1 本例检验实序列的性质DFT[xec(n)]=Re[X(k)] DFT[xoc(n)]=Im[X(k)] % 设 x(n)=10*(0.8).^n 0<=n<=10 将x(n)分解为共扼对称及共扼反对称部分%实验报告要求：（1）将实验结果画出（2）实验结果说明什么n=0:10;x=10*(0.8).^n;[xec,xoc]=circevod(x);subplot(2,1,1);stem(n,xec); %画出序列的共扼对称分量title('Circular -even component')xlabel('n');ylabel('xec(n)');axis([-0.5,10.5,-1,11])subplot(2,1,2);stem(n,xoc); %画出序列的共扼反对称分量title('Circular -odd component')xlabel('n');ylabel('xoc(n)');axis([-0.5,10.5,-4,4])figure(2)X=dft(x,11); %求出序列的DFTXec=dft(xec,11); %求序列的共扼对称分量的DFTXoc=dft(xoc,11); %求序列的共扼反对称分量的DFTsubplot(2,2,1);stem(n,real(X));axis([-0.5,10.5,-5,50])title('Real{DFT[x(n)]}');xlabel('k'); %画出序列DFT的实部subplot(2,2,2);stem(n,imag(X));axis([-0.5,10.5,-20,20])title('Imag{DFT[x(n)]}');xlabel('k'); %画出序列DFT的虚部subplot(2,2,3);stem(n,Xec);axis([-0.5,10.5,-5,50])title('DFT[xec(n)]');xlabel('k');subplot(2,2,4);stem(n,imag(Xoc));axis([-0.5,10.5,-20,20])title('DFT[xoc(n)]');xlabel('k');% 例2 本例为计算序列的圆周卷积程序% 运行之前应在命令窗口输入 x1,x2,N 的值%实验报告要求：自己选择2个序列进行计算，将实验结果写出if length(x1)>Nerror('N must be >= the length of x1')endif length(x2)>Nerror('N must be >= the length of x2')endx1=[x1 zeros(1,N-length(x1))]; %将x1,x2补0成为N长序列x2=[x2 zeros(1,N-length(x2))];m=[0:1:N-1];x2=x2(mod(-m,N)+1); %该语句的功能是将序列翻褶，延拓，取主值序列H=zeros(N,N);for n=1:1:N %该for循环的功能是得到x2序列的循环移位矩阵H(n,:)=cirshftt(x2,n-1,N); %和我们手工计算圆周卷积得到的表是一致的endy=x1*H' %用矩阵相乘的方法得到结果% 例3 本例验证采样定理%令||1000)(t a e t x -=，绘制其傅立叶变换)(Ωj X a 。