离散傅里叶变换和快速傅里叶变换

第二次实验离散傅里叶变换和快速傅里叶变换1.2掌握快速傅里叶变换（FFT）的原理和实现，二、实验原理序列x(n)的离散事件傅里叶变换（DTFT）表示为n=-伪x(n)的离散傅里叶变换（DFT）表达式为可以完成由一组有限个信号采样值x(n)直接计算得到一组有限个频谱采样值X(k)。X(k)的幅度谱为X(k)=X2(k)+X2(k)，其中下标R和I分别表示取实部、虚部的运算。X(k)的相位谱为RII(k)。X(k)R离散傅里叶反变换（IDFT）定义为N和对称性以及一些特殊值来减少DFT的运算量，可使DFT的运算量下降几个数量级，从而使数字信号处且抗混叠滤波器的截止频率不得高于与采样频率的一半。三、实验内容与相关分析（共6道）说明：为了便于老师查看，现将各题的内容写在这里——3.1求有限长离散时间信号x(n)的离散时间傅里叶变换（DTFT）X(ejΩ)并绘图。l0其他【解答】思路：这是DTFT的变换，按照定义编写DTFT的M文件即可。算机计算，计算时只取三个周期[-2π,4π]中均匀的1000个点用于绘图。理论计算的各序列DTFT表达式，请见本题的分析。M文件源代码（我的Matlab源文件不支持中文注释，抱歉functionDTFT(n1,n2,x)%ThisisaDTFTfunctionformyexperimentofSignalProcessing&Analysis.w=0:2*pi/1000:2*pi;%Definethebracketofomegaforplotting.X=zeros(size(w));%DefinetheinitialvaluesofX.fori=n1:n2X=X+x(i-n1+1)*exp((-1)*j*w*i);%ItisthedefinitionofDTFT.endAmp=abs(X);%Acquiretheamplification.Phs=angle(X);%Acquirethephaseangle(radian).subplot(1,2,1);plot(w,Amp,'r');xlabel('\Omega');ylabel('Amplification');holdon;%Plotamplificationontheleft.subplot(1,2,2);plot(w,Phs,'b');xlabel('\Omega');ylabel('PhaseAngle(radian)');holdoff;%Plotphaseangleontheright.end命令窗口中的运行及其结果（理论计算的各序列DTFT表达式，请见本题的分析第（1）小题>>n=(-2:2);>>x=1.^n;>>DTFT(-2,2,x);A543210nAnes43210>>n=(0:10);>>x=2.^n;>>DTFT(0,10,x);A【分析】 Ω43210Ω它的幅度频谱相应地应该是周期连续信号。事实上，我们可计算出它的表达式：X(Ω)=x(n)e-jΩn=Σ2e-jΩn-5jΩ)=常1-e-jΩ11-e-5jΩ1-e-jΩ和次极大，两个主极大间有|5-1|=4个极小，即有3个次级大。而对于它的相位频谱，则是周期性地在-π、对于第（2）小题，由于是离散非周期的信号。它的幅度频谱相应地应该是周期连续信号。而它的表Σx(n)e-jΩn1-2e-jΩ1-2e-jΩ1-2e-jΩ,因此主极大之间只有而由DTFT的定义可知，频谱都是以2π为周期向两边无限延伸的。由于DTFT是连续谱，对于计算机处理来说特别困难，因此我们才需要离散信号的频谱也离散，由此构造出DFT（以及为加速计算DFT【解答】M文件源代码：functionDFT(N,x)%ThisisaDFTfunctionformyexperimentofSignalProcessing&Analysis.k=(0:N-1);%DefinevariablekforDFT.X=zeros(size(k));%DefinetheinitialvalvesofX.fori=0:N-1X=X+x(i+1)*exp((-1)*j*2*k*pi/N*i);%ItisthedefinitionofDFT.endAmp=abs(X);%Acquiretheamplification.Phs=angle(X);%Acquirethephaseangle(radian).subplot(1,2,1);stem(k,Amp,'.',’MarkerSize’,18);xlabel('k');ylabel('Amplification');holdon;%Plotamplificationontheleft.subplot(1,2,2);stem(k,Phs,'*');xlabel('k');ylabel('PhaseAngle(radian)');holdoff;%Plotphaseangleontheright.endfunctionmyFFT(N,x)%Thisisabase-2FFTfunction.lov=(0:N-1);j1=0;fori=1:N%indexedaddressingifi<j1+1temp=x(j1+1);x(j1+1)=x(i);x(i)=temp;endk=N/2;whilek<=j1j1=j1-k;k=k/2;endj1=j1+k;enddigit=0;k=N;whilek>1digit=digit+1;k=k/2;endn=N/2;%Nowwestartthe"butterfly-shaped"process.formu=1:digitdif=2^(mu-1);%Differncebetweentheindexesofthetargetvariables.idx=1;fori=1:nidx1=idx;idx2=1;forj1=1:N/(2*n)r=(idx2-1)*2^(digit-mu);wn=exp(j*(-2)*pi*r/N);%Itisthe"circulatingcoefficients".temp=x(idx);x(idx)=temp+x(idx+dif)*wn;x(idx+dif)=temp-x(idx+dif)*wn;idx=idx+1;idx2=idx2+1;endidx=idx1+2*dif;endn=n/2;endAmp=abs(x);%Acquiretheamplification.Phs=angle(x);%Acquirethephaseangle(radian).subplot(1,2,1);stem(lov,Amp,'.',’MarkerSize’,18);xlabel('FFTk');ylabel('FFTAmplification');holdon;%Plottheamplification.subplot(1,2,2);stem(lov,Phs,'*');xlabel('FFTk');ylabel('FFTPhaseAngle(radian)');holdoff;end>>x=[8,7,9,5,1,7,9,5];>>DFT(8,x);nAaaAes3210kkk>>x=[8,7,9,5,1,7,9,5];>>myFFT(8,x);nAFaaAesTFF3210【分析】DFT是离散信号、离散频谱之间的映射。在这里我们可以看到序列的频谱也被离散化。事实上，我们首先，将序列做N=8周期延拓，成为离散周期信号。然后利用DFS计算得到延拓后的频谱：--=Σx(n)ejkn--的主值区间得到DFT，与图一致。因此计算正确。的幅度、相位图，并将结果与X(k)的幅度、相位图，并将结果与X(ω)相比较。【解答】M文件源代码：i)采样函数：functionxs=sampJune3(N,Ts,g)%ThisisafunctionappliedtoProblem3&4.%N:numberofsamplingpoints;Ts:samplingperiod;g=0:Nogaussiannoise;g=1:gussiannoiseexists.n=1:N;fori=1:N%Notethatimuststartat0inanalysis.ThusImadeaadaptation.sample(i)=3*cos(8*pi*Ts*(i-1))+g*randn;%InMatlab,indexstartsat1,whichisnotconsistentwithourhabit.ThusImadeaadaptationincodes.%Itisasamplingprocesswith(g=1)/without(g=0)noise.endxs=sample;plot(n-1,sample,'.',’MarkerSize’,18);xlabel('n');ylabel('value');holdoff;%Mustuse(n-1),becauseinMatlab,indexstartsat1.endii)DFT和基2-FFT函数的代码，请见第3.2节。不需再新编一个。命令窗口中的运行及其结果：>>xs=sampJune3(12,0.1,0);%末尾的0表示无噪声。>>DFT(12,xs);>>myFFT(12,xs);321 v246n8246n82.522.520nAnAa88Aes6420kkknF864200aaAesTFF2100>>xs=sampJune3(20,0.1,0);%末尾的0表示无噪声。>>DFT(20,xs);>>myFFT(20,xs); v32100n443322nAanAa0A0es505kk4332naa FATFF0505321 v5n5n44332nAnAaAeAes50kkknF500aaAesTFF43210【分析】可见原信号被0.1s采样周期的采样信号离散化之后，谱线以20π为周期重复，并且只在（20k±8）π倍，儿12点采样频谱的主谱线高度是次谱线高度的3倍。可见，在无法保证采样时间是信号周期整数倍的情况下，增加采样时间有助于减轻频谱泄漏的程度。3.4对第3步中所述连续时间信号叠加高斯白噪声信号【解答】不需要再新编程序。可以直接引用上面的函数：sampJune3(N,Ts,g)，取g=1，以体现存在白噪声DFT(N,x)myFFT(N,x)>>xs=sampJune3(12,0.1,1);%末尾的1表示有噪声。>>DFT(12,xs);>>myFFT(12,xs);euv543210nnA86420kaaAes33210knF86420aaAesTFF33210>>xs=sampJune3(20,0.1,1);%末尾的1表示有噪声。>>DFT(20,xs);>>myFFT(20,xs);nnA4321 vn50kaaAes33210knF50aAesTFF33210>>xs=sampJune3(28,0.1,0);%末尾的1表示有噪声。>>DFT(28,xs);>>myFFT(28,xs); v4321005nnAAaaAes3210nnATFFaaAesTFF43210由于叠加了噪声，所以频谱都受到了一定的干扰。由于白噪声在各个频率的功率相等，因此频谱上各以明显区分出原信号的谱线。第二好的是28点采样，因为采样时间较长，即使存在频谱泄漏也能较好地区分原信号的谱线。而最差的是12点采样，由于噪声的存在和严重的频谱泄漏，它的次谱线与主谱线的高度相差不大，使原信号不明显。jNk。N6【解答】调整既可用于计算。首先生成x(n)的6点DFT，再按照各小题分别转换，最后求相应的IDFT。M文件源代码：functionX=ExportDFT(N,x)%ThisisaDFTfunctionthatexportsthesolutiontovectorX.k=(0:N-1);X=zeros(size(k));fori=0:N-1X=X+x(i+1)*exp((-1)*j*2*k*pi/N*i);%DefinevariablekforDFT.%DefinetheinitialvalvesofX.%ItisthedefinitionofDFT.endendii)算第（1）小题的Y(k)的函数：functionY=Convertor1(X)%Thisisamathematicalconvertorforthesubproblem(1).fork=1:6Y(k)=exp((-j)*2*pi*(-4*(k-1))/6)*X(k);%Herewemustuse(k-1),becauseinMatlabindexstartsat1.endendiii)算第（2）小题的W(k)的函数：functionW=Convertor2(X)%Thisisamathematicalconvertorforthesubproblem(2).W=real(X);%AcquiretherealpartofX.endfunctionQ=Convertor3(X)%Thisisamathematicalconvertorforthesubproblem(3).Q=zeros(3);%Detailedexplanationgoesherefortmp=1:3Q(tmp)=X(2*tmp);endendfunctionx=ExportIDFT(N,X)%ThisistheIDFTfunctionformyexperiment.n=(0:N-1);%DefinevariablenforIDFT.x=zeros(size(n));%Definetheinitialvalvesofx.fork=0:N-1x=x+X(k+1)*exp(j*2*k*pi/N*n);endx=x/N;a=real(x);real(x),thoughwemayALREADYknowxisreal.%It'scausedbynumericcalculation(notanalyticcalculation)inMatlab.stem(n,a,'.','MarkerSize',18);xlabel('n');ylabel('Solution');end>>x=[4,3,2,1,0,0];>>X=ExportDFT(6,x);第（1）小题>>Y=Convertor1(X);>>y=ExportIDFT(6,Y)y=Columns1through40.0000+0.0000i0.0000+0.0000i4.0000+0.0000i3.0000+0.0000iColumns5through62.0000+0.0000i1.0000-0.0000i%虚部都是0，说明是实数%事实上，在Matlab中，由于数值计算的截断误差，对原复数做乘法后，答案的虚部可能有一极小的量。noiS43210012n345>>W=Convertor2(X);>>w=ExportIDFT(6,W)Columns1through44.0000+0.0000i1.5000+0.0000i1.0000+0.0000i1.0000+0.0000iColumns5through61.0000+0.0000i1.5000+0.0000i%虚部都是0，说明是实数；%事实上，在Matlab中，由于数值计算的截断误差，对原复数做乘法后，答案的虚部可能有一极小的量。noiS43.532.5210.5n>>Q=Convertor3(X);>>q=ExportIDFT(6,Q)q=Columns1through41.5000-0.0000i-0.1667-0.2887i0.7500-1.2990i0.8333-0.0000iColumns5through6-0.5000-0.8660i1.0833-1.8764i【分析】第（1）小题，根据DFT的性质可以判断，对原频谱乘上旋转因子W-4k之后进行IDFT得到的y(n)，6N第（2）小题，由于对原频谱取了实部，那么根据DFT的奇偶虚实性知，得到的w(n)也是实数的。第（3）小题，对原信号进行了尺度变换（“抽取”导致丢失了一些谱线，使得无法通过IDFT得到原来的序列x(n)。说明频谱记录了原有信号的信息，若频谱发生变化，则对应的时域信号也随之改变。2t)3t)，其中f1=4Hz、f2=4.02Hz、f3=5Hz，采用采样频【解答】M文件源代码：i)采样函数（其中Plus表示采样后补零的个数）functionxs=sampJune6(N,Plus)%ThisisafunctionappliedtoProblem6%N:samlingpoints;Plus:quantityofadditinalzeropoints.Ts=1/20;w1=2*pi*4;w2=2*pi*4.02;w3=2*pi*5;sample=zeros(1,N+Plus);n=(1:N);sample=sin(w1*Ts*(n-1))+sin(w2*Ts*(n-1))+sin(w3*Ts*(n-1));forp=(N+1):(N+Plus)sample(p)=0;%Addzeropoints.endxs=sample;%ReturnendfunctionmyFFT(N,x)%Thisisabase-2FFTfunctionwrotebymyself....Amp=abs(x);%Acquiretheamplification.Phs=angle(x);%Acquirethephaseangle(radian).subplot(1,2,1);stem(lov,Amp,'.');xlabel('FFTk');ylabel('FFTAmplification');holdon;%Plottheamplification.subplot(1,2,2);stem(lov,Phs,'*');xlabel('FFTk');ylabel('FFTPhaseAngle(radian)');holdoff;endfunctionmyFFT(N,x)%Thisisabase-2FFTfunctionwrotebymyself....Amp=abs(x);%Acquiretheamplification.stem(lov,Amp,'.');xlabel('FFTk');ylabel('FFTAmplification');%Plottheamplification.end>>x512=sampJune6(512,0);>>x2048=sampJune6(2048,0);>>myFFT(512,x512);>>myFFT(2048,x2048);512Points300250200nATATFF500050100150200250300350400450500nnATFFnF0>>x32p1N=sampJune6(32,32*1);>>x32p4N=sampJune6(32,32*4);>>x32p8N=sampJune6(32,32*8);>>x32p16N=sampJune6(32,32*16);>>myFFT(32+32*1,x32p1N);>>myFFT(32+32*4,x32p4N);>>myFFT(32+32*8,x32p8N);>>myFFT(32+32*16,x32p16N);%32点采样，补零N=32个，共64个数据点%32点采样，补零4N=128个，共160个数据点%32点采样，补零8N=256个，共288个数据点%32点采样，补零16N=128个，共544个数据点50nATFF5000nATFF5035302520nAT35302520nATFF500200300400500600FFTk0【分析】请注意，题目只要求绘制幅度频谱。首先，由于采样时间都不是原有信号周期的整数倍，两个采样方式对应的频谱均发生了泄漏。不过由首先，由于采样时间都不是原有信号周期的整数倍，频谱均发生了泄漏。而且由于采样时间较短，频谱泄漏比第（1）小题的两个频谱更加严重。最后，虽然都是32点采样，但由于补0个数的不同，各频谱谱线间距各不相同。例如，补零最多的序列一共有544个数据点，因此谱线间距小。由此还可以得出结论：数据点个数越多，则频谱越倾向于连可见，当采样时间不是原信号周期整数倍而且采样时间较短时，频谱泄漏都有了幅值即能量，可见当取样信号的样点数取得不够时，原信号所携带的信息就不能被完全取得。而若4.1关于各个实验的分析，请见第3部分每道题的末尾。分注意到了这个问题。4.3由于Matlab进行数值计算的过程中存在截断误差，所以最后算得的值并不是准确值。例如，对一个复数z，即使f(z)的虚部为零，但由于截断误差的存在（特别是z的虚部为无穷4.4通过利用软件对离散信号的各种变换DTFT、DFT以及其快速算法FFT进行计算，使得在实验中比较难以实现的信号分析过程（离散信号的采集和显示都是比较困难的）在计算机计算中实现，证明了理论的正确性，说明仿真计算是一种十分有效的辅助手段。（2）采样频率要高于原始信号的奈奎斯特频率。对于频谱不受限的信号，为了避免频谱