快速傅里叶变换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。
运行结果与我用numpy自带的fft()方法结果相差无几
标签：虚部,变换,傅里叶,实部,Wn,x2,x1,data,语言版 来源： https://www.cnblogs.com/The-Shining/p/11613090.html