1 #include 
     
        2 #include 
      
         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
       
        >= 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.3Hz150     {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 }