FFT算法实现实验报告

FFT算法实现实验报告
辛旸 PB10210006
实验目的
1、加深对快速傅里叶变换的理解。

2、掌握FFT 算法及其程序的编写。

3、掌握算法性能评测的方法。

实验内容
1.编写自己的FFT算法：
代码如下：
function [ X ] = Sampling( x,N )
%myFFT 实现FFT时域取样算法
% 输出：生成FFT序列X(k),输入：欲变换序列x（n），FFT变换长度N（可缺省）
（1） if ~exist('N','var'); %检查是否有变换长度N输入
（2） N=length(x); %若无，则令N等于序列长度
（3） end
（4） if N<length(x); %如果N小于序列长度，则对序列进行截短
（5） x=x(:,1:N);
（6） else
（7） x=[x,zeros(1,N-length(x))]; %如果N大于序列长度，对序列补零进行延长
（8） end;
（9） for i=1:1:length(x)/2+1; %判断N是2的多少次方
（10） if 2^i>=length(x); %若N不是2的整数幂
（11） N=2^i; %增大N为2的整数幂
（12） break;
（13） end
（14） end
（15） x=[x,zeros(1,N-length(x))]; %确保要变换的序列长度为2^i
（16） k1=zeros(1,N);
（17） X=zeros(1,N);
（18） w=zeros(1,N);
（19） for m=0:1:N-1; %确定反序序列k1和正序序列k的关系
（20） k=m;
（21） for n=i-1:-1:0; %从高位开始依次将各位移至反序位
（22） k1(m+1)=k1(m+1)+fix(k/(2^n))*(2^(i-1-n));
（23） k=rem(k,2^n);
（24） end;
（25） end
（26） for l=1:1:N;
（27） X(k1(l)+1)=x(l); %生成反序序列
（28） w(l)=exp(-1i*2*pi/N*(l-1)); %生成旋转因子
（29） end
（30） for l=0:1:i-1; %控制FFT运算级数
（31） for m=1:1:N; %每一级中有N/2个蝶形运算
（32） if rem((m-1),2^(l+1))<2^l; %找到蝶形运算的上半部分（33） b=X(m)+X(m+2^l)*w(2^(i-1-l)*rem((m-1),2^l)+1);
%将结果暂存至b
（34）
X(m+2^l)=a(m)-X(m+2^l)*w(2^(i-1-l)*rem((m-1),2^l)+1);
（35） X(m)=b; %实现原位运算
（36） end
（37） end
（38） end
2.选择实验1中的典型信号序列验证算法的有效性：
为方便比较两个算法，编写了myCompare函数计算两种算法的运行时间，并绘制频谱曲线
代码如下：
function [ t1,t2,e ] = myCompare( x,N )
%myCompare函数：比较自己编写的算法与系统自带算法的差异
% 输入：与变换信号序列x（n）和欲变换长度N
% 输出：自己编写的函数的执行时间t1，系统自带函数的执行时间t2，两者计算序列的差异平方和e tic;
X1=myFFT(x,N);
t1=toc;
tic;
X2=fft(x,N);
t2=toc;
subplot(1,2,1);plot(abs(X1));xlabel('k');ylabel('X(k');title('ÓÃ×Ô¼º±àдµÄº¯ÊýµÃµ½µÄ±ä»»ÐòÁÐƵÆ×');
subplot(1,2,2);plot(abs(X2));xlabel('k');ylabel('X(k');title('ϵͳ×Ô´øFFTº¯ÊýµÃµ½µÄ±ä»»ÐòÁÐƵÆ×');
e=sum((X1-X2).^2);
end
对理想采样信号A=444.128，α=50*2^(1/2)*π，Ω=50*2^(1/2)*π，T=1/1000，序列长度50，用自己编写的FFT算法和系统自带算法做64点FFT变换后绘制频域序列，如下：
对高斯序列，p=8，q=8，序列长度16，用自己编写的FFT算法和系统自带算法做16点FFT变换后绘制频域序列，如下：
对衰减正弦序列α=0.01，f=0.05，序列长度100，用自己编写的FFT算法和系统自带算法做128点FFT变换后绘制频域序列，如下：
由以上结果可知，自己编写的算法运行结果与系统自带算法一致，且可以对信号进行截断或补零后再做变换。

3.对所编制FFT算法进行性能评估：
与自己编写的DFT算法进行性能比较：
对N点序列进行DFT变换需要N²/2次复乘，而对N点序列做基2-FFT只需N/2*log2（N）次复乘，因此运算量减少了很多，且随着序列长度增加，运算量差异变大。

与系统自带FFT算法进行性能比较：
由于系统自带FFT函数用C语言实现，无法查看源代码，只知道效率更高，而且在计算任意点的DFT（不指定变换长度N）时，系统自带函数无需采取补零操作，而自己编写的函数会先补零再变换，改变了频域取样密度，会得到与系统自
带函数不同的结果。

比较自己编写的DFT算法、自己编写的FFT算法和系统自带算法三者运算时间，得到下表：
由此绘制曲线如图：
由比较结果可知，虽然运算时间曲线和理想曲线不完全相同，但三种算法的相对运算时间与理论预期一致。

实验总结及个人结论：
从对实验的比较结果中可知，自己编写的FFT算法有效且比自己编写的DFT算法快很多，但却始终比系统自带的算法慢，一方面是因为系统自带算法实现效率高，另一方面，在观察了自己编写的算法中各步执行时间后，发现函数用在生成反序序列的时间和实际进行FFT运算的时间相当，相当于多用了一倍的时间。

另外，对任意长序列，自己编写的FFT算法只能补零后计算，而不能像系统自带函数一样算出实际的N点DFT。