傅立叶变换—FFT

FFT(快速傅立叶变换)使用“分而治之”的策略来计算一个n阶多项式的n阶DFT系数的值。定义n为2的整数幂数，为了计算一个n阶多项式f（x）,算法定义了连个新的n/2阶多项式，函数f^[0](x)包含了f(x)中的x偶次幂项，函数f^[1](x)f(x)中的x奇次幂项。

f^[0]=a₀+a₂x+a₄x²+ ...+a_n-2x^n/2-1

f^[1]=a₁+a₃x+a₅x²+ ...+a_n-1x^n/2-1

则f(x) = f^[0](x²)+ xf^[1](x²),因此w_n⁰,w_n¹,...w_n^n-1点计算f(x)的值得问题转化成计算f^[0]和f^[1]在(w_n⁰)²,(w_n¹)²,...(w_n^n-1)²点的问题，然后计算f(x) = f^[0](x²)+ xf^[1](x²)。

FFT Code：

#include "stdio.h"
#include "math.h"

#define  LENGTH 4
#define  PI 3.1415926

typedef struct Complex
{
  float real;
  float img;
}Complex;

void Recursive_FFT(float *a,Complex *y,int len);
Complex Mul(Complex w,Complex y1_var);
Complex Add(Complex y0_var,int op,Complex mul_result );

int main()
{
  float a[LENGTH] = {1,2,3,4};
  
  Complex f[LENGTH];
  Recursive_FFT(a,f,LENGTH);
  
  int i;
  for(i=0;i<LENGTH;i++)
  {
    if(f[i].real !=0)
    {
      printf("%3.1f",f[i].real);
    }
    if(f[i].img !=0)
    {
      printf("+%3.1fi",f[i].img);
    } 
    printf("
");
  }
}

//递归求解，a为输入的初始矩阵，y为计算出来的频率矩阵
void Recursive_FFT(float *a,Complex *y,int len)
{
  Complex w0,wn;
  Complex y0[len/2],y1[len/2];
  
  w0.real = 1.0;
  w0.img = 0.0;
  
  wn.real = cos(-2 * PI /(float) len);
  wn.img = sin(-2 * PI / (float) len);
  
  float a0[len/2];
  float a1[len/2];
  int count_a0 = 0;
  int count_a1 = 0;
  
  int i;
  if(len == 1)
  {
     y[0].real = a[0];
     y[0].img = 0;
  }
  else
  {
    for(i=0;i<len;i++)
    {
      if(i % 2)
      {
    a0[count_a0++] = a[i];
      }
      else
      {
    a1[count_a1++] = a[i];
      }
    }
    
    Recursive_FFT(a0,y0,len/2);
    Recursive_FFT(a1,y1,len/2);
    
    int k;
    Complex w = w0;;
    for(k=0;k<len/2;k++)
    {
      y[k] = Add(y0[k],1,Mul(w,y1[k])); 
      y[k+len/2] = Add(y0[k],-1,Mul(w,y1[k]));
      w = Mul(w,wn);
    }
  }
  
}

//乘法运算
Complex Mul(Complex w,Complex y1_var)
{
  Complex result;
  result.real = w.real * y1_var.real - w.img * y1_var.img;
  result.img = w.real * y1_var.img + w.img * y1_var.real;
  return result;
}

//op为1则为加法运算，-1为减法运算
Complex Add(Complex y0_var,int op,Complex mul_result )
{
  Complex result;
  if(op == 1)
  {
    result.real = y0_var.real + mul_result.real;
    result.img = y0_var.img + mul_result.img;
  }
  else
  {
    result.real = y0_var.real - mul_result.real;
    result.img = y0_var.img - mul_result.img;
  }
  
  return result;
}

时间复杂度：O（n*logn）。