实验2-matlab中基2-DIT-FFT的实现

电子科技大学实验报告学生姓名：学号：2023013080指导教师一、实验室名称：数字信号处理实验室二、实验工程名称：FFT的实现三、实验原理：FFT算法思想：DFT的定义：对于有限长离散数字信号{x[n]}，0nN-1，其离散谱{x[k]}可以由离散付氏变换〔DFT〕求得。DFT的定义为：，k=0,1,…N-1通常令，称为旋转因子。直接计算DFT的问题及FFT的根本思想：由DFT的定义可以看出，在x[n]为复数序列的情况下，完全直接运算N点DFT需要〔N-1〕2次复数乘法和N〔N-1〕次加法。因此，对于一些相当大的N值〔如1024〕来说，直接计算它的DFT所作的计算量是很大的。FFT的根本思想在于，将原有的N点序列分成两个较短的序列，这些序列的DFT可以很简单的组合起来得到原序列的DFT。例如，假设N为偶数，将原有的N点序列分成两个〔N/2〕点序列，那么计算N点DFT将只需要约[(N/2)2·2]=N2/2次复数乘法。即比直接计算少作一半乘法。因子〔N/2〕2表示直接计算〔N/2〕点DFT所需要的乘法次数，而乘数2代表必须完成两个DFT。上述处理方法可以反复使用，即〔N/2〕点的DFT计算也可以化成两个〔N/4〕点的DFT〔假定N/2为偶数〕，从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT运算的情况。基2按时间抽取〔DIT〕的FFT算法思想：设序列长度为，L为整数〔如果序列长度不满足此条件，通过在后面补零让其满足〕。将长度为的序列，先按n的奇偶分成两组：，r=0,1,…,N/2-1DFT化为：上式中利用了旋转因子的可约性，即：。又令，那么上式可以写成：(k=0,1,…,N/2-1)可以看出，分别为从中取出的N/2点偶数点和奇数点序列的N/2点DFT值，所以，一个N点序列的DFT可以用两个N/2点序列的DFT组合而成。但是，从上式可以看出，这样的组合仅表示出了前N/2点的DFT值，还需要继续利用表示的后半段本算法推导才完整。利用旋转因子的周期性，有：，那么后半段的DFT值表达式：，同样，〔k=0,1,…,N/2-1〕，所以后半段〔k=N/2,…,N-1〕的DFT值可以用前半段k值表达式获得，中间还利用到，得到后半段的值表达式为：〔k=0,1,…,N/2-1〕。这样，通过计算两个N/2点序列的N/2点DFT，可以组合得到N点序列的DFT值，其组合过程如下列图所示：-1比方，一个N=8点的FFT运算按照这种方法来计算FFT可以用下面的流程图来表示：基2按频率抽取〔DIF〕的FFT算法思想：设序列长度为，L为整数〔如果序列长度不满足此条件，通过在后面补零让其满足〕。在把按k的奇偶分组之前，把输入按n的顺序分成前后两半：因为，那么有，所以：按k的奇偶来讨论，k为偶数时：k为奇数时：前面已经推导过，所以上面的两个等式可以写为：通过上面的推导，的偶数点值和奇数点值分别可以由组合而成的N/2点的序列来求得，其中偶数点值为输入x[n]的前半段和后半段之和序列的N/2点DFT值，奇数点值为输入x[n]的前半段和后半段之差再与相乘序列的N/2点DFT值。令，，那么有：这样，也可以用两个N/2点DFT来组合成一个N点DFT，组合过程如下列图所示：-1在FFT计算中使用到的MATLAB命令：函数fft(x)可以计算R点序列的R点DFT值；而fft(x,N)那么计算R点序列的N点DFT，假设R>N，那么直接截取R点DFT的前N点，假设R<N，那么x先进行补零扩展为N点序列再求N点DFT。函数ifft(X)可以计算R点的谱序列的R点IDFT值；而ifft(X,N)同fft(x,N)的情况。四、实验目的：离散傅氏变换〔DFT〕的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆DFT变换到时域。FFT是DFT的一种快速算法。在数字信号处理系统中，FFT作为一个非常重要的工具经常使用，甚至成为DSP运算能力的一个考核因素。本实验通过直接计算DFT，利用FFT算法思想计算DFT，以及使用MATLAB函数中的FFT命令计算离散时间信号的频谱，以加深对离散信号的DFT变换及FFT算法的理解。五、实验内容：计算实数序列的256点DFT。计算周期为1kHz的方波序列〔占空比为50％，幅度取为＋/-512，采样频率为25kHz，取256点长度〕256点DFT。六、实验器材〔设备、元器件〕：安装MATLAB软件的PC机一台，DSP实验演示系统一套。七、实验步骤：先利用DFT定义式，编程直接计算2个要求序列的DFT值。利用MATLAB中提供的FFT函数，计算2个要求序列的DFT值。(拓展要求)不改变序列的点数，仅改变DFT计算点数〔如变为计算1024点DFT值〕，观察画出来的频谱与前面频谱的差异，并解释这种差异。通过这一步骤的分析，理解频谱分辨力的概念，解释如何提高频谱分辨力。利用FFT的根本思想〔基2－DIT或基2－DIF〕，自己编写FFT计算函数，并用该函数计算要求序列的DFT值。并对前面3个结果进行比照。(拓展要求)尝试对其他快速傅立叶变换算法〔如Goertzel算法〕进行MATLAB编程实现，并用它来计算要求的序列的DFT值。并与前面的结果进行比照。〔拓展要求〕在提供的DSP实验板上演示要求的2种序列的FFT算法〔基2－DIT〕，用示波器观察实际计算出来的频谱结果，并与理论结果比照。八、实验数据及结果分析：程序：对要求的2种序列直接进行DFT计算的程序%第一种序列的计算N=0:255;X=cos(5*pi*N/16);fora=1:256Y(a)=0;forb=1:256Y(a)=Y(a)+X(b)*exp(-j*2*pi*(b-1)*(a-1)/256);endendsubplot(2,1,1)stem(N,abs(Y))title('DFTµÄ½á¹û')subplot(2,1,2)Y2=fft(X);stem(N,Y2)title('FFTµÄ½á¹û')%第二种序列的计算N=0:1/(1000*25):255/(1000*25);X=512*square(2*pi*N*1000);fora=1:256Y(a)=0;forb=1:256Y(a)=Y(a)+X(b)*exp(-j*2*pi*(b-1)*(a-1)/256);endend[Y]%»­fftºÍÉÏÃæµÄ½á¹ûµÄ¶Ô±Èf=0:255;Y1=fft(X);subplot(2,1,1)stem(f,Y)title('DFTµÄ½á¹û')subplot(2,1,2)stem(f,Y1)title('FFTµÄ½á¹û')对要求的2种序列进行基2－DIT和基2－DIFFFT算法程序%基-2DIT-FFT的算法%»ù-2-DIT-FFTclearclctic%½«¶ÔÊäÈëÐòÁÐx²¹ÁãÖÁ2^M¸öN=input('ÇëÊäÈëµãÊýN=');x=input('ÇëÊäÈëÐòÁÐx(ÆäÖÐÒѾ­¶¨ÒåÁËn=0:N-1)=');M=nextpow2(length(x));N=2^M;n=0:N-1;x=[x,zeros(1,N-length(x))];%¶ÔÊäÈëµÄÐòÁÐx½øÐÐÖØÐÂÅÅÐòLH=N/2;j1=LH;N1=N-2;fori=1:N1if(i<j1)T=x(i+1);x(i+1)=x(j1+1);x(j1+1)=T;endk=LH;while(j1>=k)j1=j1-k;k=k/2;endj1=j1+k;end%½øÐеûÐÎÔËËãforL=1:M%¼¶ÊýµÄÑ­»·B=2^(L-1);fori=0:B-1%ͬһ¼¶ÐýתÒò×ÓµÄÑ­»·p=i*2^(M-L);fork=i:2^L:(2^(M-L)-1)*(2^L)+i%ͬһ¼¶Í¬Ò»¸öÐýתÒò×ÓµÄÑ­»·.£¨ÔÚL¼¶ÖÐͬһ¸öÐýתÒò×Ó³öÏÖ2^(M-L)´Î£©temp=x(k+1);x(k+1)=temp+x(k+1+B)*exp(-j*2*pi*p/N);x(k+1+B)=temp-x(k+1+B)*exp(-j*2*pi*p/N);endendendstem(n,x)title('»ù-2-DIT-FFTµÄ½á¹û')time=toc%基-2-DIF-FFT的算法%»ù-2-DIT-FFTclearclctic%½«¶ÔÊäÈëÐòÁÐx²¹ÁãÖÁ2^M¸öN=input('ÇëÊäÈëµãÊýN=');x=input('ÇëÊäÈëÐòÁÐx(ÆäÖÐÒѾ­¶¨ÒåÁËn=0:N-1)=');M=nextpow2(length(x));N=2^M;n=0:N-1;x=[x,zeros(1,N-length(x))];%¶ÔÊäÈëµÄÐòÁÐx½øÐÐÖØÐÂÅÅÐòLH=N/2;j1=LH;N1=N-2;fori=1:N1if(i<j1)T=x(i+1);x(i+1)=x(j1+1);x(j1+1)=T;endk=LH;while(j1>=k)j1=j1-k;k=k/2;endj1=j1+k;end%½øÐеûÐÎÔËËãforL=1:M%¼¶ÊýµÄÑ­»·B=2^(L-1);fori=0:B-1%ͬһ¼¶ÐýתÒò×ÓµÄÑ­»·p=i*2^(M-L);fork=i:2^L:(2^(M-L)-1)*(2^L)+i%ͬһ¼¶Í¬Ò»¸öÐýתÒò×ÓµÄÑ­»·.£¨ÔÚL¼¶ÖÐͬһ¸öÐýתÒò×Ó³öÏÖ2^(M-L)´Î£©temp=x(k+1);x(k+1)=temp+x(k+1+B)*exp(-j*2*pi*p/N);x(k+1+B)=temp-x(k+1+B)*exp(-j*2*pi*p/N);endendendstem(n,x)title('»ù-2-DIT-FFTµÄ½á¹û')time=toc对要求的2种序列用MATLAB中提供的FFT函数进行计算的程序N=0:255;w=0:4095;X=512*square(2*pi*1000*N/25000);y1=fft(X);y2=fft(X,4096);plot(N*2/256,abs(y1),'o',w*2/4096,abs(y2))结果：〔1〕对2种要求的序列直接进行DFT计算的频域波形〔2〕对2种要求的序列进行基2－DIT和基2－DIFFFT算法频域波形〔3〕对2种要求的序列用MATLAB中提供的FFT函数计算的频域波形。〔4〕〔拓展要求〕分析利用上面的方法画出的信号频谱与理论计算出来的频谱之间的差异，并解释这种差异。答：对于第一个函数，理论上的频谱应该有很多的褶皱的。而用DFT画出来的频谱就只有两点才有值，其他的值都为0，这就是因为采样没有采到其它的点。〔5〕〔拓展要求〕保持序列点数不变，改变DFT计算点数〔变为1024点〕，观察频谱的变化，并分析这种变化，由此讨论如何提高频谱分辨力的问题。答：如下列图所示，我们将第一个函数的采样点提高了。这个时候，我们可以看到其实提高点数，就是增加在频域的采样点数，在频域0到2*pi范围内的采样点数。我们增加FFT变换的点数，实际上就是增加了采样率。九、实验结