快速傅里叶变换FFT

快速傅里叶变换FFT

DFT是信号分析与处理中的一种重要变换。但直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。
1.直接计算DFT
长度为N的有限长序列x(n)的DFT为：

2.减少运算量的思路和方法
思路：N点DFT的复乘次数等于N2。把N点DFT分解为几个较短的DFT，可使乘法次数大大减少。另外，旋转因子WmN具有周期性和对称性。
(考虑x(n)为复数序列的一般情况，对某一个k值，直接按上式计算X(k)值需要N次复数乘法、(N-1)次复数加法.)
方法：
分解N为较小值：把序列分解为几个较短的序列，分别计算其DFT值，可使乘法次数大大减少；
利用旋转因子WNk的周期性、对称性进行合并、归类处理，以减少DFT的运算次数。
周期性：
对称性：
3.FFT算法思想
不断地把长序列的DFT分解成几个短序列的DFT，并利用旋转因子的周期性和对称性来减少DFT的运算次数。








再次分解，对N=8点，可分解三次。


c语言程序:

// FFT.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<windows.h>
#define N 1000
typedef struct{
    double real;
    double img;
 }complex;
void fft();/*快速傅里叶变换*/
void ifft();
void initW();
void change();
void add(complex, complex, complex *);/*复数加法*/
void mul(complex,complex,complex *);    /*复数乘法*/
void sub(complex,complex,complex *);/*复数减法*/
void divi(complex, complex,complex *);/*复数除法*/
void output();   /*输出结果*/
complex x[N], *W;/*输出序列的值*/
int size_x = 0;/*输入序列的长度，只限2的N次方*/
double PI;
int main(){
    int i, method;
    system("cls");
    PI = atan(1.0) * 4;
    printf("Please input the size of x: \n");/*输入序列的长度*/
    scanf("%d", &size_x);
    printf("Please input the data in x[N](such as:5 6):\n");
    /*输入序列对应的值*/
    for (i = 0; i<size_x; i++)
          scanf("%lf %lf", &x[i].real, &x[i].img);
    initW();
    /*选择FFT或逆FFT运算*/
    printf("Use FFT(0) or IFFT(1) ?\n");
    scanf("%d", &method);
      if (method == 0)
          fft();
      else
          ifft();
      output();
      system("pause");
      return 0;
}
    /*进行基-2 FFT运算*/
void fft(){
      int i = 0, j = 0, k = 0, l = 0;
      complex up, down, product;
      change();
      for (i = 0; i < log((float)size_x) / log(2.0); i++)
      {
          l = 1 << i;
          for (j = 0; j<size_x; j += 2 * l)
          {
              for (k = 0; k<l; k++)
              {
                  mul(x[j + k + l], W[size_x*k / 2 / l], &product);
                  add(x[j + k], product, &up);
                  sub(x[j + k], product, &down);
                  x[j + k] = up;
                  x[j + k + l] = down;
              }
          }
      }
  }
void ifft()
{
    int i = 0, j = 0, k = 0;

    complex up, down;
    for (i = 0; i < log((float)size_x) / log(2.0); i++)  /*蝶形运算*/
    {
        int l = size_x;
        l/= 2;
        for (j = 0; j < size_x; j += 2 * l)
        {
            for (k = 0; k < l; k++)
            {
                add(x[j + k], x[j + k + l], &up);
                up.real /= 2;
                up.img /= 2;
                sub(x[j + k], x[j + k + l], &down);
                down.real /= 2;
                down.img /= 2;
                divi(down, W[size_x*k / 2 / l], &down);
                x[j + k] = up;
                x[j + k + l] = down;
            }
        }
    }
    change();
}
/*初始化变化核*/
void   initW()
{
    int i;
    W = (complex   *)malloc(sizeof(complex)*   size_x);
    for (i = 0; i < size_x; i++)
    {
        W[i].real = cos(2 * PI / size_x*i);
        W[i].img = -1 * sin(2 * PI / size_x*i);
    }
}
    /*变址计算，将x(n)码位倒置*/
    void change()
    {
        complex temp;
        unsigned short i = 0, j = 0, k = 0;
        double t;
        for (i = 0; i<size_x; i++)
        {
            k = i;
            j = 0;
            t = (log((float)size_x) / log(2.0));
            while ((t--)>0)
            {
                j = j << 1;
                j |= (k & 1);
                k = k >> 1;
            }
            if (j > i)
            {
                temp = x[i];
                x[i] = x[j];
                x[j] = temp;
            }
        }
}
void output()/*输出结果*/
{
        int i;
        printf("The result are as follows: \n");
        for (i = 0; i<size_x; i++)
        {
            printf("%.4f", x[i].real);
            if (x[i].img >= 0.0001)
                printf("+%.4fi\n", x[i].img);
            else if(fabs(x[i].img)<0.0001)
                printf("\n");
            else
                printf("%.4fi\n", x[i].img);
        }
    }
void add(complex a, complex b, complex *c)
{
        c->real = a.real + b.real;
        c->img = a.img + b.img;
}
void mul(complex a, complex b, complex *c)
{
        c->real = a.real*b.real- a.img*b.img;
        c->img = a.real*b.img+ a.img*b.real;
}
void sub(complex a, complex b, complex *c)
{
        c->real = a.real - b.real;
        c->img = a.img - b.img;
}
void divi(complex a, complex b, complex *c)
{
    c->real = (a.real*b.real + a.img*b.img) / (b.real*b.real + b.img*b.img);
    c->img = (a.img*b.real - a.real*b.img) / (b.real*b.real + b.img*b.img);
}