FFT并行算法与应用-基于MPI(一)

首先，介绍下FFT算法。
参考：
http://www.gatevin.moe/acm/fft%E7%AE%97%E6%B3%95%E5%AD%A6%E4%B9%A0%E7%AC%94%E8%AE%B0/
博客写的很好，也很基础。

简单介绍以下FFT：

为什么需要FFT

　　第一个问题是为什么要创造FFT，简单的说，为了速度。我们承认DFT很有用，但是我们发现他的速度不是很快，1D的DFT原始算法的时间复杂度是O(n^2)，这个可以通过公式观察出来，对于2D的DFT其时间复杂度是O（n^4），这个速度真的很难接受，也就是说，你计算一幅1024x768的图像时，你将等大概。。。大概。。。我也没试过，反正很久。不信的自己去试试。所以找到一种快速方法的方法计算FFT势在必行。
　　以下为DFT公式：
这里写图片描述

　　计算一个4点DFT。计算量如下：
这里写图片描述

如何得到FFT

    　　通过观察DFT公式，我们发现DFT计算每一项X[u]的时候都遍历了完整的x[n]所以，我们的想法就是能不能通过其他的X[u+m]（m为一个整数，可正可负）得到X[u]，也就是充分利用之间的计算结构来构建现在的结果，这种方法就很容易表现成迭代的形式。本文介绍基2的FFT，及离散信号x[n]的个数为2的k次方，即如果完整的离散信号中有N个信分量{x1，x2，x3....xN}其中`(N=1<<k)`。

　　数学基础
　　FFT的数学基础其实就是：
这里写图片描述

　　旋转因子具有以下性质：

　　这些性质使得我们可以利用前面的计算结果来完成出后续的计算。

2-FFT

　　观察：一个只有两个值的离散信号，假设N=2，利用性质2-对称性可以得到

4-FFT

　　与上面2点的一样，推广到4点，将会是这样，其中方框内的操作为上述2-FFT：

　　最终算法：
这里写图片描述

8-FFT

　　废话不多说，和上面一样：
　　计算过程为：
这里写图片描述
　　结果为：

2^n-FFT

　　同理，推广到2^n，可以得到类似的结果，不过我们发现，为了使输出结果为顺序结果，输入的顺序经过了一系列的调整，而且每一步WN的幂次参数也是变化，我们必须得到其变化规律才能更好的编写程序：
　　观察WN的变化规律为：
这里写图片描述

　　节点距离（设为z）就是从WN的0次幂开始连续加到WN的z次幂。然后间隔为z个式子，再次从0次幂加到z次。重复以上过程：
这里写图片描述

　　例如第二级，间隔z=2，节点为实心黑色点，节点0,1，不做操作，节点2*W8^0,节点3*W8^2，节点4，5不做操作。。
　　其内在规律就是节点i是否乘以WN：
　　if（i%z==奇数）
　　节点i*WN^(step)；
　　step每次增加的数量由当前的所在的计算级决定
　　step=1<<（k（总级数）-i（当前级数））
　　输入参数重新排序
　　i行表示数组下标，蓝色数字表示实际数据：其内在规律就是，下标为偶数的将被映射到自己本身下标的1/2处，下表为奇数时被映射到数组后半段（size_n/2）的（下标-1）1/2处，将排列后的数据分为前后两个部分，递归次过程，直到只有两个元素。则停止。过程如下：

　　观察算法
　　观察至此，我们已经基本把FFT蝶形算法的所有特征都搞定了，接下来就是使用代码来实现了。
这里写图片描述
fft.c

#include "math.h"  
#include "fft.h"  
//精度0.0001弧度  

void conjugate_complex(int n,complex in[],complex out[])  
{  
  int i = 0;  
  for(i=0;i<n;i++)  
  {  
    out[i].imag = -in[i].imag;  
    out[i].real = in[i].real;  
  }   
}  

void c_abs(complex f[],float out[],int n)  
{  
  int i = 0;  
  float t;  
  for(i=0;i<n;i++)  
  {  
    t = f[i].real * f[i].real + f[i].imag * f[i].imag;  
    out[i] = sqrt(t);  
  }   
}  


void c_plus(complex a,complex b,complex *c)  
{  
  c->real = a.real + b.real;  
  c->imag = a.imag + b.imag;  
}  

void c_sub(complex a,complex b,complex *c)  
{  
  c->real = a.real - b.real;  
  c->imag = a.imag - b.imag;   
}  

void c_mul(complex a,complex b,complex *c)  
{  
  c->real = a.real * b.real - a.imag * b.imag;  
  c->imag = a.real * b.imag + a.imag * b.real;     
}  

void c_div(complex a,complex b,complex *c)  
{  
  c->real = (a.real * b.real + a.imag * b.imag)/(b.real * b.real +b.imag * b.imag);  
  c->imag = (a.imag * b.real - a.real * b.imag)/(b.real * b.real +b.imag * b.imag);  
}  

#define SWAP(a,b)  tempr=(a);(a)=(b);(b)=tempr  

void Wn_i(int n,int i,complex *Wn,char flag)  
{  
  Wn->real = cos(2*PI*i/n);  
  if(flag == 1)  
  Wn->imag = -sin(2*PI*i/n);  
  else if(flag == 0)  
  Wn->imag = -sin(2*PI*i/n);  
}  

//傅里叶变化  
void fft(int N,complex f[])  
{  
  complex t,wn;//中间变量  
  int i,j,k,m,n,l,r,M;  
  int la,lb,lc;  
  /*----计算分解的级数M=log2(N)----*/  
  for(i=N,M=1;(i=i/2)!=1;M++);   
  /*----按照倒位序重新排列原信号----*/  
  for(i=1,j=N/2;i<=N-2;i++)  
  {  
    if(i<j)  
    {  
      t=f[j];  
      f[j]=f[i];  
      f[i]=t;  
    }  
    k=N/2;  
    while(k<=j)  
    {  
      j=j-k;  
      k=k/2;  
    }  
    j=j+k;  
  }  

  /*----FFT算法----*/  
  for(m=1;m<=M;m++)  
  {  
    la=pow(2,m); //la=2^m代表第m级每个分组所含节点数       
    lb=la/2;    //lb代表第m级每个分组所含碟形单元数  
                 //同时它也表示每个碟形单元上下节点之间的距离  
    /*----碟形运算----*/  
    for(l=1;l<=lb;l++)  
    {  
      r=(l-1)*pow(2,M-m);     
      for(n=l-1;n<N-1;n=n+la) //遍历每个分组，分组总数为N/la  
      {  
        lc=n+lb;  //n,lc分别代表一个碟形单元的上、下节点编号       
        Wn_i(N,r,&wn,1);//wn=Wnr  
        c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算  
        c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr  
        c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr  
      }  
    }  
  }  
}  

//傅里叶逆变换  
void ifft(int N,complex f[])  
{  
  int i=0;  
  conjugate_complex(N,f,f);  
  fft(N,f);  
  conjugate_complex(N,f,f);  
  for(i=0;i<N;i++)  
  {  
    f[i].imag = (f[i].imag)/N;  
    f[i].real = (f[i].real)/N;  
  }  
}

fft.h

#ifndef __FFT_H__  
#define __FFT_H__  

typedef struct complex //复数类型  
{  
  float real;       //实部  
  float imag;       //虚部  
}complex;  

#define PI 3.1415926535897932384626433832795028841971  


///////////////////////////////////////////  
void conjugate_complex(int n,complex in[],complex out[]);  
void c_plus(complex a,complex b,complex *c);//复数加  
void c_mul(complex a,complex b,complex *c) ;//复数乘  
void c_sub(complex a,complex b,complex *c); //复数减法  
void c_div(complex a,complex b,complex *c); //复数除法  
void fft(int N,complex f[]);//傅立叶变换 输出也存在数组f中  
void ifft(int N,complex f[]); // 傅里叶逆变换  
void c_abs(complex f[],float out[],int n);//复数数组取模  
////////////////////////////////////////////  
#endif