快速傅里叶变换c语言版

  1 /*
  2 用法：
  3      通过调用FFT(double *data)函数，在此函数中的临时复数数组_data[]就是
  4      傅里叶变换后的结果。想要得到幅值显示可以用sqrt()函数。幅值就是实部
  5      的平方加上虚部的平方再开方。
  6      注：宏定义中lengh必须要和传入的数组data的长度一样长，此函数支持的
  7      傅里叶变换的位数为2,4,8,16,32,64。理论上可以应用于单片机上（这本来
  8      也是出发目的）
  9 */
 10 #include "stdafx.h"
 11 #include <iostream>
 12 #include <stdio.h>
 13 #include <stdlib.h>
 14 #define lengh 8
 15 #define pi 3.1415
 16 
 17 typedef struct complex
 18 {
 19     double real;
 20     double imag;
 21 }complex;
 22 
 23 void daoxu(int a[])
 24 {
 25     int i,j,N = 0;
 26     unsigned int _daoxu,k;
 27     for (i = 0; i < 8; i++) {         //算出是多少次方  ，最多频谱是64位
 28         if ((1 << i) == lengh)
 29             N = i;
 30     }
 31 
 32     for (i = 0; i < lengh; i++) {    //初始化查表空间
 33         k = i;
 34         _daoxu = 0;
 35         for (j = 0; j < N; j++) {    //倒序
 36             _daoxu = _daoxu << 1;
 37             _daoxu += (1 & k);
 38             k = k >> 1;
 39         }
 40         a[i] = _daoxu;
 41     }
 42 }
 43 /*
 44    参数：data是输入数据数组实数
 45    返回值：无
 46 */
 47 void FFT(double *data)  
 48 {
 49     int i, j, N = 0,x1,x2,step,_lengh;
 50     int daoxu_data[lengh] = { 0 };
 51     complex Wn, tem_variable,_data[lengh];
 52 
 53 
 54     daoxu(daoxu_data);      //生成倒序表
 55     
 56 
 57     for (i = 0; i < 8; i++) {         //算出是多少次方  ，最多频谱是64位
 58         if ((1 << i) == lengh)
 59             N = i;
 60     }
 61     for (i = 0; i < lengh; i++) {         //把实数转换成复数,并且把倒序后的数据传过去
 62         _data[i].real = data[daoxu_data[i]];
 63         _data[i].imag = 0;
 64         printf("_data[%d].real:%f\n", i, _data[i].real);
 65     }
 66 
 67     for (i = 0; i < N; i++) {  //分级
 68         step = 1 << (i + 1);    //确定运算几大步，
 69         Wn.real = 1.0;
 70         Wn.imag = 0.0;
 71         printf("第几次蝶形运算：%d\n", i);
 72         for (j = 0; j < (step / 2); j++) {
 73             for (x1 = j; x1 < lengh; x1 += step) {
 74                 x2 = x1 + (step / 2);
 75 
 76 
 77                 tem_variable.real = 0.0;
 78                 tem_variable.imag = 0.0;
 79 
 80                 tem_variable.real = (_data[x1].real + _data[x2].real*Wn.real - _data[x2].imag*Wn.imag);
 81                 tem_variable.imag = (_data[x1].imag + _data[x2].real*Wn.imag + _data[x2].imag*Wn.real);
 82 
 83                 _data[x2].real = (_data[x1].real - _data[x2].real*Wn.real + _data[x2].imag*Wn.imag);
 84                 _data[x2].imag = (_data[x1].imag - _data[x2].real*Wn.imag - _data[x2].imag*Wn.real);
 85                 _data[x1].real = tem_variable.real;
 86                 _data[x1].imag = tem_variable.imag;
 87                 printf("x1:%d,x2:%d\n", x1, x2);
 88             }
 89             _lengh = lengh;
 90             printf("Wn实部：%f,Wn的虚部：%f\n", Wn.real, Wn.imag);
 91             Wn.real = (cos(2 * pi*(j + 1) / (_lengh >> (N - (i + 1)))));
 92             Wn.imag = -(sin(2 * pi*(j + 1) / (_lengh >> (N - (i + 1)))));
 93 
 94         }
 95     }
 96     for (i = 0; i < lengh; i++) {
 97         printf("实部：%f ,虚部：%f。\n", _data[i].real,_data[i].imag);
 98     }
 99 }
100 
101 int main(void)
102 {
103     double data[] = {0,1,2,3,4,5,6,7,8};
104     FFT(data);
105     return 0;
106 }

_data[0].real:0.000000
_data[1].real:4.000000
_data[2].real:2.000000
_data[3].real:6.000000
_data[4].real:1.000000
_data[5].real:5.000000
_data[6].real:3.000000
_data[7].real:7.000000
第几次蝶形运算：0
x1:0,x2:1
x1:2,x2:3
x1:4,x2:5
x1:6,x2:7
Wn实部：1.000000,Wn的虚部：0.000000
第几次蝶形运算：1
x1:0,x2:2
x1:4,x2:6
Wn实部：1.000000,Wn的虚部：0.000000
x1:1,x2:3
x1:5,x2:7
Wn实部：0.000046,Wn的虚部：-1.000000
第几次蝶形运算：2
x1:0,x2:4
Wn实部：1.000000,Wn的虚部：0.000000
x1:1,x2:5
Wn实部：0.707123,Wn的虚部：-0.707090
x1:2,x2:6
Wn实部：0.000046,Wn的虚部：-1.000000
x1:3,x2:7
Wn实部：-0.707058,Wn的虚部：-0.707156
实部：28.000000 ,虚部：0.000000。
实部：-4.000447 ,虚部：9.656985。
实部：-4.000185 ,虚部：4.000000。
实部：-4.000208 ,虚部：1.656777。
实部：-4.000000 ,虚部：0.000000。
实部：-3.999923 ,虚部：-1.656800。
实部：-3.999815 ,虚部：-3.999815。
实部：-3.999422 ,虚部：-9.656129。