快速傅里叶变换的递归实现

转载于：https://blog.csdn.net/a493823882/article/details/78731540

快速傅里叶变换算法原理简述 基于递归的fft实现

2017年12月06日 16:27:42 SethChai 阅读数 7620

一、FFT：Fast Fourier Transform

按定义式计算DFT的复杂度为O(n^2)，我们可以使用FFT将复杂度优化到O(nlog2N)。

长度为8的FFT的过程由如下图展示，称为蝶形运算：

对于长为N=2^n的序列x[n]，根据DFT：

令

那么DFT可写为：

所以FFT的思想是每次都把序列分为奇序列与偶序列。

在接下来的程序中，把取x[n]的部分序列的编号称为index，如index=(4,3)，即对应子序列x[4n+3]。

长度为8的序列分解如下图，每个都展示了index值：

这便是FFT算法

二、FFT的递归实现

1.核心代码：

采用C++实现，核心递归解释位于程序注释

 
 

   
   

    
    

   
   
   
   

    
    

     
     #include <iostream>
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     #include <complex>
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     #include <cmath>
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     using 
     
     namespace 
     
     std;
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     #define comp complex<double>
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     #define PI acos(-1.0)
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     #define I comp(0,1)
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     const 
     
     int N=
     
     8;
     
     //设长度8
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     int x[N];
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     comp X[N];
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     struct index{
     
     //表示取x[n]的index，eg:mul=4,add=3,即取x[4n+3]
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     int mul,add;
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         index(
     
     int m,
     
     int a){
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
             mul=m;
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
             add=a;
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         }
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         index 
     
     operator+(
     
     int a){
    
    
   
   

   
   

    
    

   
   
   
   

    
    
        
     
     return index(mul,add+a);
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         }
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         index 
     
     operator*(
     
     int m){
    
    
   
   

   
   

    
    

   
   
   
   

    
    
        
     
     return index(mul*m,add);
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         }
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     };
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     comp WN(double n,double N){
     
     //旋转因子
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     return 
     
     exp(
     
     -1.0*I*
     
     2.0*PI/N*n);
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     }
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     comp fft(int k,index i)
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     {
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     if(i.mul==N)
    
    
   
   

   
   

    
    

   
   
   
   

    
    
        
     
     return x[i.add];
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     return fft(k,i*
     
     2)+fft(k,i*
     
     2+i.mul)*WN(i.mul*k,N);
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     /* 核心递归：fft的思想是把序列分为奇与偶两部分，也就是相邻的分开，
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         * 因为对于当前index为mul的序列，在x[n]中每隔一个mul取一个数，那么对于当前序列x[mul*n+add]
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         * 每个数在x[n]中相邻mul，那么令mul=m, add=a，x[mn+a]=x[2*mn+a]+x[2*mn+a+m]
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
         */
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     }
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     int main()
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     {
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     for(
     
     int i=
     
     0;i<N;i++)
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
             x[i]=i;
     
     //设序列为0,1,2,...,7
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     for(
     
     int i=
     
     0;i<N;i++)
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     cout<<fft(i,index(
     
     1,
     
     0))<<
     
     endl;
    
    
   
   

   
   

    
    

   
   
   
   

    
    
    
     
     return 
     
     0;
    
    
   
   

   
   

    
    

   
   
   
   

    
    

     
     }
    
    
   
   

   
   

    
    

   
   
   
   

    
    
 
    
    
   
   
 
 

注：用complex类运行fft时间较长，只作为演示代码，不用于使用。

2.正确检验

长度为8的序列0~7的结果：

Matlab运行结果：

长度为16的序列0~15的结果：

Matlab运行结果：

可见程序可靠。

3.递归顺序展示：

展示了计算一次X[k]的值，每次程序输出的是当前递归的index值