快速傅里叶变换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 1000typedef 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);}