傅里叶变换
1 #include <stdio.h> 2 #include <math.h> 3 4 typedef char int8_t; 5 typedef unsigned char uint8_t; 6 typedef unsigned short uint16_t; 7 typedef short int16_t; 8 9 #define N 64 //傅里叶变换的点数 10 #define M 6 //蝶形运算的级数,N = 2^M 11 #define N2 32 //N/2 12 #define N4 16 //N/4 13 14 #define PI 3.14159 //圆周率 15 #define PI2 6.28318 16 17 //定义复数结构体 18 typedef struct 19 { 20 float real; //实部 21 float imag; //虚部 22 }complex; 23 24 //傅里叶变换序列,一维 25 float pr[64] = 26 { 27 1234,1285,1151,1086,896,671,541,392,368,565,762,905,987,1103,1352,1601,1593,1565,1576,1565,1379,1152,1208,1150,1086,1124,1092,945,791,561,393,291,352,596,950,1307,1490,1707,1760,1494,1351,1510,1427,1191,1156,1090,953,917,944,847,655,580,457,495,679,877,1015,1084,1256,1519,1469,1411,1423,1509 28 }; 29 30 //正弦函数表 31 float sin_table[N4 + 1]; 32 33 34 float find_sin(float x) 35 { 36 int i = ((int)(N * x)) >> 1; 37 38 if (i > N4) //不会超过N/2 39 i = N2 - i; 40 41 return sin_table[i]; 42 } 43 44 float find_cos(float x) 45 { 46 int i = ((int)(N * x)) >> 1; 47 48 if (i < N4) 49 return sin_table[N4 - i]; 50 else //i>Npart4 && i<N/2 51 return -sin_table[i - N4]; 52 } 53 54 //dir:1 - 傅里叶变换, -1 - 傅里叶逆变换 55 void fft(complex data[], uint8_t num, uint8_t n, int8_t dir) 56 { 57 uint8_t i = 0, j = 0, k = 0, m = 0, t = num - 1; 58 float angle; 59 complex w, temp; 60 61 //变址,把自然顺序变为倒位序 62 for (i = 0; i < t; i++) 63 { 64 if (i < j) 65 { 66 temp = data[j]; 67 data[j] = data[i]; 68 data[i] = temp; 69 } 70 k = N2; 71 while (k <= j) 72 { 73 j -= k; 74 k >>= 1; 75 } 76 77 j += k; 78 } 79 80 for (i = 1; i <= n; i++) 81 { 82 uint8_t km, step = 1 << i; 83 m = step >> 1; 84 85 for (j = 0; j < m; j++) 86 { 87 angle = ((double)j) / m; 88 89 if (dir == 1) 90 w.imag = -find_sin(angle); 91 else if (dir == -1) 92 w.imag = find_sin(angle); 93 w.real = find_cos(angle); 94 95 for (k = j; k < num; k += step) 96 { 97 km = k + m; 98 //用下面两行直接计算复数乘法,省去函数调用开销 99 temp.real = data[km].real * w.real - data[km].imag * w.imag; 100 temp.imag = w.imag * data[km].real + w.real * data[km].imag; 101 102 data[km].real = data[k].real - temp.real; 103 data[km].imag = data[k].imag - temp.imag; 104 105 data[k].real = data[k].real + temp.real; 106 data[k].imag = data[k].imag + temp.imag; 107 } 108 } 109 } 110 111 if (dir == -1) 112 { 113 for (i = 0; i < num; i++) 114 { 115 data[i].real /= num; 116 } 117 } 118 } 119 120 121 int main() 122 { 123 int i; 124 complex data[N]; 125 126 //初始化输入数据 127 for (i = 0; i < N; i++) 128 { 129 data[i].real = pr[i]; 130 data[i].imag = 0; 131 } 132 133 //创建正弦表 134 for (i = 0; i <= N4; i++) 135 { 136 sin_table[i] = (float)sin(PI * i /N2); 137 } 138 139 //正变换 140 fft(data, N, M, 1); 141 142 for (i = 0; i < N; i++) 143 { 144 printf("%.f %.f\n", data[i].real, data[i].imag); 145 } 146 printf("\n"); 147 148 //波形处理 149 for (i = 5; i < N; i++) //10:3.9Hz,11:4.3Hz 150 { 151 data[i].imag= 0; 152 data[i].real = 0; 153 } 154 155 //逆变换 156 fft(data, N, M, -1); 157 158 for (i = 0; i < N; i++) 159 { 160 printf("%.f %.f\n", data[i].real, data[i].imag); 161 } 162 printf("\n"); 163 164 return 0; 165 }