数字信号处理

第一次讨论课内容（1）时域信号（a）如何由模拟信号产生时域离散信号；模拟信号：信号的自变量和函数值都为连续值。一段连续的时间间隔内，其代表信息的特征量可以在任意瞬间呈现为任意数值的信号。时域离散信号：自变量取离散值，函数值取连续值。离散信号是一个序列，即其自变量是“离散”的。这个序列的每一个值都可以被看作是连续信号的一个采样。模拟信号，以采样间隔T对它进行等间隔采样，得到时域离散信号。即： (取整数， )采样是将时间上、幅值上都连续的模拟信号，在采样脉冲的作用，转换成时间上离散（时间上有固定间隔）、但幅值上仍连续的离散模拟信号。所以采样又称为波形的离散化过程。 对模拟信号进行采样可以看作一个模拟信号通过一个电子开关S。设电子开关的作用等效成一宽度为，周期为T的矩形脉冲串，采样信号就是与相乘的结果。理想采样：0，脉冲串变为单位冲激串。中每个单位冲激处在采样点上，强度为1，理想采样是与相乘的结果。即： 式中是单位冲激信号，上式只有在时，才可能是非零值，所以可写成:由模拟信号经采样产生时域离散信号的MATLAB程序： 一个连续的周期性三角波信号频率为50HZ，信号幅度在0+2V之间，在窗口显示2个周期信号波形，对信号的一个周期进行16点采样来获取离散信号。代码：f=50;Um=1;ts=2; %输入信号的频率、振幅、显示周期N=16; %信号一个采样周期的采样点数为16T=1/f; %信号周期Tdt=T/N; %采样时间间隔n=0:ts*N-1; %建立离散时间的时间序列列tn=n*dt; %确定时间序列样点在时间轴上的位置y=Um*sawtooth(2*f*pi*tn,1/2)+1; %三角波subplot(2,1,1),stem(tn,y); %显示采样后的信号title(离散信号);subplot(2,1,2),plot(tn,y); %显示原连续信号title(连续信号);运行结果： 一个连续的周期性正弦信号频率为200Hz，信号幅度在-1+1V之间，在窗口上显示2个周期信号波形，用Fs=4KHZ的频率对连续信号进行采样，试显示连续信号和采样获得的离散信号波形。代码：f=200;Um=1;nt=2; %输入信号的频率、振幅、显示周期Fs=4000; %采样频率N=Fs/f; %采样点数T=1/f; %信号周期 dt=T/N; %采样时间间隔n=0:nt*N-1; %建立离散时间的时间序列tn=n*dt; %确定时间序列样点在时间轴上的位置y=Um*sin(2*f*pi*tn); %正弦波subplot(2,1,1),stem(tn,y); %显示采样后的信号title(离散信号);subplot(2,1,2),plot(tn,y); %显示原连续信号title(连续信号);运行结果：抽样定理是连续时间信号和离散时间信号之间的桥梁，在时域该系统实现了输入信号与抽样序列的相乘，完成了时间轴的离散，在频域实现了原信号频谱的周期延拓。在奈奎斯特抽样定理的条件下(抽样频率不小于被抽样带限信号最高频率的两倍)，一个连续时间信号完全可以用该信号在等时间间隔点上的样本来表示，在频率轴上实现了原信号频谱无混叠的周期化。采样定理：设连续信号的的最高频率为，如果采样频率，那么采样信号可以唯一的恢复出原连续信号，否则会造成采样信号中的频谱混叠现象，不可能无失真地恢复原连续信号。验证采样定理的MATLAB程序：画出连续时间信号的时域波形及其幅频特性曲线，信号为 对信号进行采样，得到采样序列，绘制采样频率分别为80Hz，120 Hz，150 Hz时的采样序列波形；、对不同采样频率下的采样序列进行频谱分析，绘制其幅频曲线，对比各频率下采样序列和的幅频曲线有无差别。由采样序列恢复出连续时间信号，画出其时域波形，对比与原连续时间信号的时域波形。代码：%实现采样频谱分析绘图函数function fz=caiyang(fy,fs) %第一个输入变量是原信号函数，信号函数fy以字符串的格式输入%第二个输入变量是采样频率fs0=10000; tp=0.1; t=-tp:1/fs0:tp; k1=0:999; k2=-999:-1; m1=length(k1); m2=length(k2); f=fs0*k2/m2,fs0*k1/m1; %设置原信号的频率数组w=-2*pi*k2/m2,2*pi*k1/m1; fx1=eval(fy); FX1=fx1*exp(-j*1:length(fx1)*w); %求原信号的离散时间傅里叶变换figure % 画原信号波形subplot(2,1,1),plot(t,fx1,r) title(原信号), xlabel(时间t (s) axis(min(t),max(t),min(fx1),max(fx1)% 画原信号幅度频谱subplot(2,1,2),plot(f,abs(FX1),r) title(原信号幅度频谱) , xlabel(频率f (Hz) axis(-100,100,0,max(abs(FX1)+5) % 对信号进行采样Ts=1/fs; %采样周期t1=-tp:Ts:tp; %采样时间序列f1=fs*k2/m2,fs*k1/m1; %设置采样信号的频率数组t=t1; %变量替换fz=eval(fy);%获取采样序列FZ=fz*exp(-j*1:length(fz)*w); %采样信号的离散时间傅里叶变换figure % 画采样序列波形subplot(2,1,1),stem(t,fz,.), title(取样信号) , xlabel(时间t (s) line(min(t),max(t),0,0) % 画采样信号幅度频谱subplot(2,1,2),plot(f1,abs(FZ),m) title(取样信号幅度频谱) , xlabel(频率f (Hz)function fh=huifu(fz,fs)%第一个输入变量是采样序列%第二个输入变量是得到采样序列所用的采样频率T=1/fs; dt=T/10; tp=0.1; t=-tp:dt:tp; n=-tp/T:tp/T; TMN=ones(length(n),1)*t-n*T*ones(1,length(t); fh=fz*sinc(fs*TMN); % 由采样信号恢复原信号k1=0:999; k2=-999:-1; m1=length(k1); m2=length(k2); w=-2*pi*k2/m2,2*pi*k1/m1; FH=fh*exp(-j*1:length(fh)*w); % 恢复后的信号的离散时间傅里叶变换figure % 画恢复后的信号的波形subplot(2,1,1),plot(t,fh,g), st1=sprintf(由取样频率fs=%d,fs); st2=恢复后的信号; st=st1,st2; title(st) , xlabel(时间t (s) axis(min(t),max(t),min(fh),max(fh) line(min(t),max(t),0,0) % 画重构信号的幅度频谱f=10*fs*k2/m2,10*fs*k1/m1; %设置频率数组subplot(2,1,2),plot(f,abs(FH),g) title(恢复后信号的频谱) , xlabel(频率f (Hz) axis(-100,100,0,max(abs(FH)+2);%主函数f1=sin(2*pi*60*t)+cos(2*pi*25*t)+cos(2*pi*30*t);%输入一个信号fs0=caiyang(f1,80);%欠采样fr0=huifu(fs0,80); fs1=caiyang(f1,120);%临界采样fr1=huifu(fs1,120); fs2=caiyang(f1,150);%过采样fr2=huifu(fs2,150); 运行结果：(1) 采样频率时，为原信号的欠采样信号和恢复，采样频率不满足时域采样定理，那么频移后的各相临频谱会发生相互重叠，这样就无法将他们分开，因而也不能再恢复原信号。频谱重叠的现象被称为混叠现象。欠采样信号的离散波形及频谱、恢复后信号见下图。(2) 采样频率 时，为原信号的临界采样信号和恢复，下图为其采样的离散波形和频谱，从下图恢复后信号和原信号先对比可知，只恢复了低频信号，高频信号未能恢复。（3）采样频率时，为原信号的过采样信号和恢复，由图采样信号离散波形和频谱，可以看出采样信号的频谱是原信号频谱进行周期延拓形成的，从图采样恢复后的波形和频谱，可看出与原信号误差很小了，说明恢复信号的精度已经很高。(b) 线性时不变系统输入和输出之间的关系设系统的输入用表示，表示成单位脉冲序列移位加权和为： 则系统输出为： 根据线性系统叠加性质：又根据时不变性质 即：线性时不变系统的输出等于输入序列和该系统的单位脉冲响应的卷积。验证例子：设模拟正弦信号：，现以T为周期，对其进行采样，得到正弦序列为： 其中为数字频率，为数字角频率。产生信号： N取大于500设离散时间线性时不变系统为：令通过此系统MATLAB程序如下：n=0:500;xa=sin(0.0364*pi*n);figure(1);stem(n,xa,.);xlabel(n);ylabel(xa(n);hn(1:50)=1;y1=conv(xa,hn);figure(2);stem(y1,.);xlabel(n);ylabel(y1(n);y2=filter(hn,1,xa);figure(3);stem(y2,.);axis(0,600,-20,20)xlabel(n);ylabel(y2(n);运行结果：的时域波性卷积conv方法：通过系统的输出波形滤波Filter函数方法：通过系统的输出波形两种函数运行结果相同，证明了线性时不变系统的输出等于输入序列和该系统的单位脉冲响应的卷积。（2）离散傅立叶变换(a) 离散傅立叶变换的性质；设是一个长度为M的有限长序列，则定义的N点离散傅里叶变换为： 1、 线性性质如果两个有限长序列分别为和x2(n)，长度分别为N1和N2，且y(n)ax1(n)bx2(n) (a、b均为常数)则该y(n)的N点DFT为Y(k)DFTy(n)aX1(k)bX2(k) 0kN1其中：NmaxN1，N2，X1(k)和X2(k)分别为x1(n)和x2(n)的N点DFT。验证线性性质例子：已知x1(n)0，1，2，4，x2(n)1，0，1，0，1，(1)y(n)2x1(n)3x2(n)，再由y(n)的N点DFT获得Y1(k)；(2)由x1(n)、x2(n)求X1(k)、X2(k)，再求Y2(k)2X1(k)3X2(k)。MATLAB程序：xn1=0,1,2,4; %建立xn1序列xn2=1,0,1,0,1;%建立xn2序列N1=length(xn1);N2=length(xn2);N=max(N1,N2); %确定Nif N1N2 xn2=xn2,zeros(1,N1-N2);%对长度短的序列补0elseif N2N1 xn1=xn1,zeros(1,N2-N1);endyn=2*xn1+3*xn2;%计算ynn=0:N-1;k=0:N-1;Yk1=yn*(exp(-j*2*pi/N).(n*k) %求yn的N点DFTXk1=xn1*(exp(-j*2*pi/N).(n*k); %求xn1的N点DFTXk2=xn2*(exp(-j*2*pi/N).(n*k); %求xn2的N点DFTYk2=2*Xk1+3*Xk2%由Xk1、Xk2求Yksubplot(3,2,1),stem(n,xn1);title(x1(n);subplot(3,2,2),stem(n,Xk1);title(X1(k);subplot(3,2,3),stem(n,xn2);title(x2(n);subplot(3,2,4),stem(n,Xk2);title(X1(k);%subplot(4,2,5),stem(n,yn);%title(yn);subplot(3,2,6),stem(n,Yk2);axis(-1,5,-20,30)title(2*Xk1+3*Xk2);subplot(3,2,5),stem(n,Yk1);axis(-1,5,-20,30)title(DFTy(n);运行结果：Yk1 = 23.0000 + 0.0000i -7.5902 + 1.5388i 3.5902 - 0.3633i 3.5902 + 0.3633i -7.5902 - 1.5388iYk2 = 23.0000 + 0.0000i -7.5902 + 1.5388i 3.5902 - 0.3633i 3.5902 + 0.3633i -7.5902 - 1.5388i由运算结果和波形图都可知，Y1(K)=Y2(K)，即y(n)2x1(n)3x2(n)，DFTy(n)= 2X1(k)3X2(k),可验证线性性质。2、 循环移位性质如果有限长序列为x(n)，长度为N，将x(n)的循环移位定义为：x(n)左移m位的过程可由以下步骤获得：(1)将x(n)以N为周期进行周期延拓，得到 ；(2)将左移m位，得到 ；(3)取的主值序列，得到x(n)循环移位序列y(n)。验证例子：已知有限长序列x(n)1，2，3，4，5，6，求x(n)左移2位成为新的向量y(n)，并画出循环移位的中间过程。MATLAB程序：xn=1,2,3,4,5,6; %建立xn序列Nx=length(xn);nx=0:Nx-1;nx1=-Nx:2*Nx-1;%设立周期延拓的范围x1=xn(mod(nx1,Nx)+1);%建立周期延拓序列ny1=nx1-2;y1=x1;%将x1左移2位,得到y1RN=(nx1=0)&(nx1=0)&(ny1Nx);%在y1的位置向量ny1上设置主值窗subplot(4,1,1),stem(nx1,RN.*x1);%画出x1的主值部分subplot(4,1,2),stem(nx1,x1);%画出x1subplot(4,1,3),stem(ny1,y1);%画出y1subplot(4,1,4),stem(ny1,RN1.*y1); %画出y1的主值部分运行结果：可见，循环移位的实质是将x(n)左移m位时，移出主值区的序列值又依次从右端进入主值区。有限长序列的移位也称为循环移位。3、时域和频域循环卷积特性离散傅里叶变换的循环卷积特性也称为圆周卷积，分为时域卷积和频域卷积两类。1)时域循环卷积假定有限长序列长度分别为N1、N2，。 2)频域循环卷积利用时域和频域的对称性，可以得到频域卷积特性。若 ，则 下面重点讨论时域循环卷积。时域循环卷积的方法有多种：方法1：直接使用时域循环卷积。由于有限长序列可以看成是周期序列的主值，因此，时域圆周卷积的结果可以由对应的周期序列卷积和取主值部分获得。方法2：用频域DFT相乘再求逆变换。即先分别求x1(n)、x2(n)的DFTX1(k)、X2(k)，再求Y(k)的IDFT获得y(n)。方法3：用FFT和IFFT进行循环卷积。基本思路同方法2，但直接使用了MATLAB提供的fft和ifft子函数来实现。举例：已知的两个时域周期序列分别取主值，得到x11，1，1，0，0，0，x20，1，2，3，0，0，求时域循环卷积y(n)并用图形表示。xn1=0,1,2,3,0,0; %建立x1(n)序列xn2=1,1,1,0,0,0; %建立x2(n)序列N=length(xn1);n=0:N-1;k=0:N-1;Xk1=xn1*(exp(-j*2*pi/N).(n*k); %由x1(n)的DFT求X1(k)Xk2=xn2*(exp(-j*2*pi/N).(n*k); %由x2(n)的DFT求X2(k)Yk=Xk1.*Xk2; %Y(k)=X1(k)X2(k)yn=Yk*(exp(j*2*pi/N).(n*k)/N; %由Y(k)的IDFT求y(n)yn=abs(yn) %取模值,消除DFT带来的微小复数影响subplot(2,3,1),stem(n,xn1);title(x1(n)subplot(2,3,2),stem(n,xn2);title(x2(n)subplot(2,3,3),stem(n,yn);title(y(n)subplot(2,3,4),stem(n,Xk1);title(X1(k)subplot(2,3,5),stem(n,Xk2);title(X2(k)subplot(2,3,6),stem(n,Yk);title(Y(k)运行结果：得到：yn = 0.0000 1.0000 3.0000 6.0000 5.0000 3.00004、共轭对称性由于序列x(n)及其离散傅里叶变换X(k)的定义在主值为0N1的区间，因此DFT的循环对称性对时间序列是指关于n0和nN/2的对称性，对频谱序列是关于数字频率为0和p的对称性。本实验重点分析实序列的循环对称性。实序列x(n)可以分解为循环偶序列xe(n)和循环奇序列xo(n)：x(n)xe(n)xo(n) 0nN1其中：设DFTx(n)X(k)ReX(k)j*ImX(k)，则有即实序列中的偶序列xe(n)对应于x(n)的离散傅里叶变换X(k)的实部，而实序列中的奇序列xo(n)对应于x(n)的离散傅里叶变换X(k)的虚部。举例：已知一个定义在主值区间的实序列xones(1，4)，zeros(1，4)，试将其分解成为偶对称序列和奇对称序列，并求它们的DFT，验证离散傅里叶变换的循环对称性。MATLAB代码：x=ones(1,5),zeros(1,5) %建立x(n)序列N=length(x);n=0:N-1;k=0:N-1;xr=x(mod(-n,N)+1); %求x(-n)xe=0.5*(x+xr) %求x(n)的偶序列xo=0.5*(x-xr) %求x(n)的奇序列X=x*(exp(-j*2*pi/N).(n*k); %由x(n)的DFT求X(k)Xe=xe*(exp(-j*2*pi/N).(n*k); %由xe(n)的DFT求Xe(k)Xo=xo*(exp(-j*2*pi/N).(n*k); %由xo(n)的DFT求Xo(k)error1=(max(abs(real(X)-Xe) %计算X(k)的实部与Xe(k)的差值error2=(max(abs(j*imag(X)-Xo) %计算X(k)的虚部与Xo(k)的差值subplot(2,4,1),stem(n,x);title(x(n)subplot(2,4,2),stem(n,xr);title(x(-n)subplot(2,4,3),stem(n,xe);title(xe(n)subplot(2,4,4),stem(n,xo);title(xo(n)subplot(2,4,5),stem(n,real(X);title(X(k)的实部)subplot(2,4,6),stem(n,imag(X);title(X(k)的虚部) subplot(2,4,7),stem(n,Xe);title(Xe(k)=DFT(xe(n)subplot(2,4,8),stem(n,Xo);title(Xo(k)=DFT(xo(n)x = 1 1 1 1 1 0 0 0 0 0xe = 1.0000 0.5000 0.5000 0.5000 0.5000 0 0.5000 0.5000 0.5000 0.5000xo = 0 0.5000 0.5000 0.5000 0.5000 0 -0.5000 -0.5000 -0.5000 -0.5000error1 = 3.8778e-15error2 = 3.9475e-15(b) 如何利用离散傅立叶变换对信号进行谱分析；谱分析的误差及改进方法；DFT的高分辨率频谱与高密度频谱之间的区别。利用离散傅立叶变换对信号进行谱分析：时域连续信号离散傅里叶分析的基本步骤如下图所示：其中低通滤波器LPF（预滤波器）的引入，是为了消除或减少时域连续信号转换成序列时，可能出现的频谱混叠的影响。实际工作中，时域离散信号x(n)的时宽是很长的甚至是无限长。由于DFT的需要，必须把x(n)限制在一定的时间间隔之内，即进行数据截断。数据的截断相当于加窗处理。因此，在计算x(n)的DFT之前，用一个时域有限的窗函数w(n)加到x(n)上。通过A/D变换器转换成采样序列x(n)时域离散序列。其频谱用表示，它是频率的周期函数，即 由于进行DFT的需要，对序列x(n)进行加窗处理，即v(n)=x(n).w(n)。加窗对频域的影响，用周期卷积表示。 最后进行DFT运算。加窗后的DFT是 例子：先用MATLAB产生数字信号：然后逐个用DFT进行谱分析，分别取DFT的长度N16，32，画出信号的幅谱图，分析实验结果。MATLAB程序：例2：先用MATLAB产生出如下模拟信号：设采样频率Hz，然后用DFT进行谱分析，分别取DFT的长度N16，32，64，画出信号的幅谱图，横轴打印模拟频率。讨论DFT长度与频率分辨率的关系，要在频谱图上分离出上述3个频率，DFT长度至少取多大？代码：N1=16;N2=32;N3=64;n=0:N1-1;Fs=64;t=n/Fs;x4=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);figure,subplot(231),stem(n,x4);gridX4=abs(fft(x4,N1);subplot(232),stem(n,X4);gridn=0:N2-1;Fs=64;t=n/Fs;x4=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);subplot(233),stem(n,x4);gridX4=abs(fft(x4,N2);subplot(234),stem(n,X4);gridn=0:N3-1;Fs=64;t=n/Fs;x4=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);subplot(235),stem(n,x4);gridaxis(0 31 0 4);X4=abs(fft(x4,N3);subplot(236),stem(n,X4);gridaxis(0 31 0 40);运行结果：可见，只有当N=64时，3个频率在频谱图上被分离出来。分析： 谱分析的误差及改进方法：（2）栅栏效应（3）截断效应DFT的高分辨率频谱与高密度频谱之间的区别:高分辨率谱通过增加信号的记录长度得到，频率分辨率F与信号的实际长度成反比。信号越长，F值越小，即分辨率（频域取样中两相邻点间的频率间隔）越高。高密度谱是指当信号的时间域长度不变时，在频域内对它的频谱进行提高采样频率而得到的高密度谱，它只可以更细化当前分辨率下的频谱，克服栅栏效应，但不能改变DFT的分辨率。另外，采用尾部补零的方法不能提高DFT的高分辨率。即增大M,在原序列后面加0，采样频率不变。区别1、高密度谱呈许多谱线型,补充0越多,谱线也越密集，但不改变频谱分辨率。2、通过增加采样点数N，可提高物理分辨率，得到高分辨率谱。实例分析：用两个不同频率的正弦信号开研究分辨率的大小，能分辨的的两个正弦信号频率越靠近，表明分辨率越高。设序列，序列有两个主要频率0.48和0.52。画出其10点DFT、64点DFT及10点后序列补零至64点的DFT。10点DFT程序如下：clearn=0:9;xn=cos(0.48*pi*n)+cos(0.52*pi*n);xn1=xn(1:1:10);N=10;n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn1*WNnk;subplot(121)stem(n,xn1);xlabel(n);ylabel(x(n);axis(0 9 -2.5 2.5);title(正弦序列);n1=0:1:9;subplot(122)plot(n1,abs(Xk(1:1:10)title(幅度谱);xlabel(k);ylabel(|X(k)|); 10点DFT正弦序列及频谱图64点DFT程序：clearn=0:1:63;xn=cos(0.48*pi*n)+cos(0.52*pi*n);xn1=xn(1:1:64);N=64;n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn1*WNnk;subplot(121)stem(n,xn1);xlabel(n);ylabel(x(n);axis(0 63 -2.5 2.5);title(正弦序列);n1=0:1:63;subplot(122)plot(n1,abs(Xk(1:1:64)title(幅度谱);xlabel(k);ylabel(|X(k)|);64点DFT正弦序列及频谱图10点序列补零值64点DFT的MATLAB程序：clearn=0:1:9;xn=cos(0.48*pi*n)+cos(0.52*pi*n);xn1=xn(1:1:10),zeros(1,54);N=64;n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn1*WNnk;subplot(121)stem(n,xn1);xlabel(n);ylabel(x(n);axis(0 63 -2.5 2.5);title(正弦序列);n1=0:1:63;subplot(122)plot(n1,abs(Xk(1:1:64)title(幅度谱);xlabel(k);ylabel(|X(k)|);10点序列补零值64点DFT正弦序列及频谱图 :DFT的高分辨率频谱与高密度频谱之间的区别:高分辨率谱通过增加信号的记录长度得到，频率分辨率F与信号的实际长度成反比。信号越长，F值越小，即分辨率（频域取样中两相邻点间的频率间隔）越高。高密度谱是指当信号的时间域长度不变时，在频域内对它的频谱进行提高采样频率而得到的高密度谱，它只可以更细化当前分辨率下的频谱，克服栅栏效应，但不能改变DFT的分辨率。另外，采用尾部补零的方法不能提高DFT的高分辨率。即增大M,在原序列后面加0，采样频率不变。区别1、高密度谱呈许多谱线型,补充0越多,谱线也越密集，但不改变频谱分辨率。2、通过增加采样点数N，可提高物理分辨率，得到高分辨率谱。已知序列： x(n)=cos(0.82n)+2sin(0.43n) ，试确定一合适样本数N，利用MATLAB的fft函数分析计算信号()xn的频谱。 解：序列x(n)是一个周期序列。为了说明高密度谱和高分辨率谱之间的区别，分以下几种情况进行讨论： 先取x(n)的前10个样本，10点DFT的MATLAB程序如下：n=0:1:9;x = cos(0.82*pi*n)+2*sin(0.43*pi*n); subplot(2,1,1); stem(n,x,.); title(x(n), 0 = n = 9); xlabel(n); ylabel(x(n); axis(0,10,-2.5,2.5); Xk = fft(x); magXk = abs(Xk(1:1:6); k1 = 0:1:5; w1 = 2*pi/10*k1; subplot(2,1,2); stem(w1/pi,magXk,.); title(Samples of DTFT Magnitude); xlabel(frequency in pi units); ylabel(|X(k)|); axis(0,1,0,10); 由于样本数不足，难以获得足够的信息而得到正确的结论。即从频谱图无法观测到原复合余弦信号x(n)的w=0.43和w=0.82两个频率分量. 在先前()xn的前10个样本后补90个零，以期得到一个更高密度的频谱。补零后程序：n = 0:1:9; x = cos(0.82*pi*n)+2*sin(0.43*pi*n); n1 = 0:1:99; x1 = x zeros(1,90); subplot(2,1,1); stem(n1,x1,.); title(x(n), 0 = n = 9 + 90 zeros); xlabel(n); ylabel(x(n); axis(0,100,-2.5,