DSP研究性学习报告频谱计算

《近代数字信号处理》课程研究性学习报告姓名舒晓萌学号08213047同组成员指导教师魏杰时间2010-4-10Matlab习题补充：M2-1利用DFT的性质，编写一MATLAB程序，计算下列序列的循环卷积。（1）g[k]={1,-3,4,2,0,-2},h[k]={3,0,1,-1,2,1};（2）x[k]=cos(pi*k/2),y=3k,k=0,1,2,3,4,5.（1）N=6;g=[1,-3,4,2,0,-2];k=[3,0,1,-1,2,1];x1=[gzeros(1,N-length(g))];x2=[kzeros(1,N-length(k))];m=0:N-1;x1=x1(mod(-m,N)+1);H=zeros(N,N);forn=1:N;H(1,:)=x1;H(n+1,:)=circshift(H(n,:),[0,1]);endy=H*x2';y=y’y=6-317-27-136自主学习方面：利用到mod求余函数，还利用到circshift求循环位移函数。b=circshift(a,[x,y])其中a为待移位的矩阵，x表示上下移位位数（正数表示向下移位），y表示左右移位位数（正数表示向右移位）（2）N=6;k=0:N-1;x=cos(pi*k/2);y=3.^k;x1=[gzeros(1,N-length(x))];x2=[kzeros(1,N-length(y))];x1=x1(mod(-k,N)+1);H=zeros(N,N);forn=1:N;H(1,:)=x1;H(n+1,:)=circshift(H(n,:),[0,1]);endy=H*x2';y=y'y=5253-7-595M2-2已知序列cos(/2),[]0,kNkNxk其他（1）计算该序列DTFT的表达式Xje，并画出N=10时，Xje的曲线。（2）编写一MATLAB程序，利用fft函数，计算N=10时，序列x[k]的DTFT在2mmN的抽样值。利用hold函数，将抽样点画在Xje曲线上。(1)N=10;k=0:N-1;L=512;x=cos(k*pi/2/N);X=fftshift(fft(x,L));plot(2*pi*(0:L-1)/L,abs(X))gridon;(2)N=10;k=0:N-1;L=512;x=cos(k*pi/2/N);X=fftshift(fft(x,L));w=2*pi*(-1/2+(0:L-1)/L);plot(w,abs(X),'r');gridon;holdon;m=N*(-1/2+(0:N-1)/N);stem(m,abs(X),'b');holdoff;怎么规定横坐标？DFT近似计算信号频谱专题研讨【目的】(1)掌握利用DFT近似计算不同类型信号频谱的原理和方法。(2)理解误差产生的原因及减小误差的方法。(3)培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】基本题M2-0已知一离散序列为31,,1,0),π2.0sin(][kkkx(1)用L=32点DFT计算该序列的频谱，求出频谱中谱峰的频率；(2)对序列进行补零，然后分别用L=64、128、256、512点DFT计算该序列的频谱，求出频谱中谱峰的频率；(3)讨论所获得的结果，给出你的结论。该结论对序列的频谱计算有何指导意义？【题目分析】本题讨论补零对离散序列频谱计算的影响。【温磬提示】在计算离散非周期序列频谱时常用/作为横坐标，称/为归一化频率normalizedfrequency)。在画频谱时需给出横坐标。每幅图下都需给出简要的文字说明。由于离散非周期序列频谱是周期的，所以在计算时不需要用fftshift函数对fft计算的结果进行重新排列。【序列频谱计算的基本方法】【仿真结果】(1)peakvalue14.5561frequencyofpeakvalue3(2)peakvalue14.5561frequencyofpeakvalue6peakvalue15.2671frequencyofpeakvalue13peakvalue15.2671frequencyofpeakvalue26peakvalue15.3563frequencyofpeakvalue51【结果分析】DFT点数越多，则L越大，即fsam/N越小，显示分辨率越高。【自主学习内容】已知幅度谱，用matlab求谱峰所对应的频率值。利用：[a1,b1]=max(X);其中X是一个向量，返回的b1值对应的是向量X中的最大值，返回的a1值为最大值b1对应的索引。[y,v]=max(a)returnstheindicesofthemaximumvaluesinvectorv.Ifthevaluesalongthefirstnonsingletondimensioncontainmorethanonemaximalelement,theindexofthefirstoneisreturned.【阅读文献】数字信号处理，陈后金。【发现问题】(专题研讨或相关知识点学习中发现的问题)：对应的题中离散非周期序列的DFT应该是离散的，应该利用stem画离散谱还是利用plot画连续谱？【问题探究】【仿真程序】(1)N=32;k=0:N-1;xk=sin(0.2*pi*k);Xm=abs(fft(xk,N));stem(2*k/N,Xm);xlabel('Normalizedfrequencyof32DFT');ylabel('Magnitude');gridon;[a1b1]=max(Xm);disp('peakvalue');disp(a1');disp('frequencyofpeakvalue');disp(k(b1)');（2）N=32k=0:N-1;xk=sin(0.2*pi*k);L=input('thelengthofDFT:L=');Xm=abs(fft(xk,L));m=0:L-1;stem(2*m/L,Xm);xlabel(['NormalizedfrequencyofL=',num2str(L)]);holdon;ylabel('Magnitude');gridon;[a1b1]=max(Xm);disp('peakvalue');disp(a1');disp('frequencyofpeakvalue');disp(m(b1)');M2-3已知一离散序列为x[k]=Acos0k+Bcos0+)k)。用长度N=64的哈明窗对信号截短后近似计算其频谱。试用不同的A和B的取值，确定用哈明窗能分辩的最小的谱峰间隔Ncπ2Δw中c的值。【题目分析】本题讨论用哈明窗计算序列频谱时的频率分辨率。【仿真结果】【结果分析】将实验结果与教材中定义的窗函数的有效宽度相比较，发表你的看法。=(-wsam*T/2+(0:L-1)*wsam*T/L),wsam*T=2m==(-1/2+(0:L-1)/L);Ncπ2Δwm=c/N对应画出的频谱中，两个谱峰横坐标的间隔设为m,m=c/N.m=0.0293,N=64c=0.0293*64=1.8752近似等于2【自主学习内容】始终令A=1，改变B的值（BA），控制变量来求得分辨情况。【阅读文献】【发现问题】(专题研讨或相关知识点学习中发现的问题)：fftshift的作用正是让正半轴部分和负半轴部分的图像分别关于各自的中心对称。因为直接用fft得出的数据与频率不是对应的，fftshift可以纠正过来【问题探究】在离散序列频谱计算中为何要用窗函数？用不同的窗函数对计算结果有何影响？与矩形窗相比哈明窗有何特点?如何选择窗函数？答：有一些离散序列长度无限长，计算机无法处理，所以要利用窗函数进行截短。用不同的窗函数得到的计算结果不同，与矩形窗相比哈明窗减小了旁瓣，却加宽了主瓣宽度。【仿真程序】A=input('A=');B=input('B=');x=A*cos(100*pi*k)+B*cos((100*pi+0.1*pi)*k);N=64;L=512;wn=(hamming(N))';xn=x.*wn;Xn=fftshift(fft(xn,L));m=(-1/2+(0:L-1)/L);plot(m,abs(Xn));xlabel('Normalizedfrequency');ylabel('Magnitude');gridon;title(['A=',num2str(A)'''B=',num2str(B)]);M2-4已知一离散序列为x[k]=cos(0k)+0.75cos(1k),0k63其中0=0.4,1=0+/64(1)对x[k]做64点FFT,画出此时信号的频谱。(2)如果(1)中显示的谱不能分辨两个谱峰，是否可对(1)中的64点信号补零而分辨出两个谱峰。通过编程进行证实，并解释其原因。(3)给出一种能分辨出信号中两个谱峰的计算方案，并进行仿真实验。【题目分析】分析影响谱峰分辨率的主要因数，进一步认识补零在在频谱计算中的作用。【仿真结果】（1）(2)（3）【结果分析】不能对(1)中的64点信号补零而分辨出两个谱峰，这样的方法只能改变屏幕分辨率。可以通过加hamming窗来实现对谱峰的分辨。【自主学习内容】【阅读文献】【发现问题】(专题研讨或相关知识点学习中发现的问题)：【问题探究】【仿真程序】（1）w1=0.4*pi;w2=w1+pi/64;N=64k=0:N-1;x=cos(w1*k)+0.75*cos(w2*k);X=fft(x);plot(k/N,abs(X));gridon;xlabel('Normalizedfrequencyof64DFT');ylabel('Magnitude');(2)w1=0.4*pi;w2=w1+pi/64;N=64;L=input('²¹ÁãÖ®ºóµÄÐòÁеij¤¶ÈL=');k=0:N-1;x=cos(w1*k)+0.75*cos(w2*k);X=fft(x,L);plot((0:L-1)/N,abs(X));gridon;xlabel(['Normalizedfrequencyof',num2str(L)]);ylabel('Magnitude');(3）w1=0.4*pi;w2=w1+pi/64;N=128;L=128;k=0:N-1;x=cos(w1*k)+0.75*cos(w2*k);wn=(boxcar(N))';xn=x.*wn;Xn=fft(xn,L);m=0:L-1;plot(m/L,abs(Xn));xlabel('Normalizedfrequency');ylabel('Magnitude');gridon;M2-6试用DFT近似计算高斯信号)exp()(2dttg的频谱抽样值。高斯信号频谱的理论值为)4exp(π)j(2ddG通过与理论值比较，讨论信号的时域截取长度和抽样频率对计算误差的影响。【题目分析】连续非周期信号频谱计算的基本方法。计算中出现误差的主要原因及减小误差的方法。【仿真结果】【结果分析】由于信号及频谱都有理论表达式，在进行误差分析时希望给出一些定量的结果。【自主学习内容】【阅读文献】【发现问题】(专题研讨或相关知识点学习中发现的问题)：【问题探究】【仿真程序】扩展题M2-7本题研究连续周期信号频谱的近似计算问题。周期为T0的连续时间周期信号x(t)可用Fourier级数表示为ntnnXtx0j0e)()(其中ttxTnXtnTde)(~1)(00j00X(n0)称为连续时间周期信号x(t)的频谱函数。000π2π2fT称为信号的基频（基波），0n称为信号的谐波。如果信号x(t)函数表达式已知，则可由积分得出信号的频谱。如果信号x(t)函数表达式未知，或者x(t)函数表达式非常复杂，则很难由积分得信号的频谱。本题的目的就是研究如何利用DFT近似计算连续时间周期信号的频谱。(1)若在信