FFT快速傅里叶变换与IFFT逆变换

1. /****************************************************************************/
   
2. /*                                                                          */
   
3. /*              快速傅里叶变换 / 快速傅里叶逆变换（实数）                   */
   
4. /*                                                                          */
   
5. /*              2015年07月21日                                              */
   
6. /*                                                                          */
   
7. /****************************************************************************/
   
8. #include                   // C 语言标准输入输出函数库
   
9. #include                    // C 数学函数库
   
10. 
    
11. #include "mathlib.h"                // DSP 数学函数库
    
12. #include "dsplib.h"                 // DSP 函数库
    
13. 
    
14. /****************************************************************************/
    
15. /*                                                                          */
    
16. /*              宏定义                                                      */
    
17. /*                                                                          */
    
18. /****************************************************************************/
    
19. // 软件断点
    
20. #define SW_BREAKPOINT     asm(" SWBP 0 ");
    
21. 
    
22. // 快速傅里叶变换
    
23. // π 及 浮点数极小值
    
24. #define PI                3.141592654
    
25. #define F_TOL             (1e-06)
    
26. 
    
27. /****************************************************************************/
    
28. /*                                                                          */
    
29. /*              全局变量                                                    */
    
30. /*                                                                          */
    
31. /****************************************************************************/
    
32. // 快速傅里叶变换测试
    
33. // 测试快速傅里叶变换点数
    
34. // 注意:TI DSP库 最大支持一次性计算 256K 个实数点的 FFT
    
35. #define N 1024
    
36. // 采样频率
    
37. #define Fs 1000.0
    
38. // 信号
    
39. #pragma DATA_ALIGN(Input, 8);
    
40. float Input[2 * N];
    
41. // 信号 副本
    
42. float InputOrig[2 * N];
    
43. // FFT 输出
    
44. #pragma DATA_ALIGN(FFT_Out, 8);
    
45. float FFT_Out[2 * N];
    
46. // IFFT 输出
    
47. #pragma DATA_ALIGN(IFFT_Out, 8);
    
48. float IFFT_Out[2 * N];
    
49. 
    
50. // 旋转因子
    
51. #pragma DATA_ALIGN(W, 8);
    
52. float W[2 * N];
    
53. 
    
54. // 二进制位翻转
    
55. #pragma DATA_ALIGN (brev, 8);
    
56. unsigned char brev[64] =
    
57. {
    
58.         0x00, 0x20, 0x10, 0x30, 0x08, 0x28, 0x18, 0x38,
    
59.         0x04, 0x24, 0x14, 0x34, 0x0c, 0x2c, 0x1c, 0x3c,
    
60.         0x02, 0x22, 0x12, 0x32, 0x0a, 0x2a, 0x1a, 0x3a,
    
61.         0x06, 0x26, 0x16, 0x36, 0x0e, 0x2e, 0x1e, 0x3e,
    
62.         0x01, 0x21, 0x11, 0x31, 0x09, 0x29, 0x19, 0x39,
    
63.         0x05, 0x25, 0x15, 0x35, 0x0d, 0x2d, 0x1d, 0x3d,
    
64.         0x03, 0x23, 0x13, 0x33, 0x0b, 0x2b, 0x1b, 0x3b,
    
65.         0x07, 0x27, 0x17, 0x37, 0x0f, 0x2f, 0x1f, 0x3f
    
66. };
    
67. 
    
68. /****************************************************************************/
    
69. /*                                                                          */
    
70. /*              函数声明                                                    */
    
71. /*                                                                          */
    
72. /****************************************************************************/
    
73. 
    
74. /****************************************************************************/
    
75. /*                                                                          */
    
76. /*              快速傅里叶变换测试                                          */
    
77. /*                                                                          */
    
78. /****************************************************************************/
    
79. // 产生旋转因子
    
80. void tw_gen(float *w, int n)
    
81. {
    
82.     int i, j, k;
    
83. 
    
84.     for(i = 1, k = 0; i > 2; i++)
    
85.     {
    
86.         w[k    ] = sin(2 * PI * i / n);
    
87.         w[k + 1] = cos(2 * PI * i / n);
    
88. 
    
89.         k += 2;
    
90.     }
    
91. 
    
92.     for(j = 1; j > 3; j = j > 3; i += j)
    
93.         {
    
94.             w[k]     = (float)sin( 4 * PI * i / n);
    
95.             w[k + 1] = (float)cos( 4 * PI * i / n);
    
96.             w[k + 2] = (float)sin( 8 * PI * i / n);
    
97.             w[k + 3] = (float)cos( 8 * PI * i / n);
    
98.             w[k + 4] = (float)sin(12 * PI * i / n);
    
99.             w[k + 5] = (float)cos(12 * PI * i / n);
    
100. 
     
101.             k += 6;
     
102.         }
     
103.     }
     
104. }
     
105. 
     
106. void tw_geni(float *w, int n)
     
107. {
     
108.     int i, j, k;
     
109. 
     
110.     for(i = 1, k = 0; i > 2; i++)
     
111.     {
     
112.         w[k    ] =  sin (2 * PI * i / n);
     
113.         w[k + 1] = -cos (2 * PI * i / n);
     
114. 
     
115.         k += 2;
     
116.     }
     
117. 
     
118.     for(j = 1; j > 3; j = j > 3; i += j)
     
119.         {
     
120.             w[k]     = (float) -sin ( 4 * PI * i / n);
     
121.             w[k + 1] = (float)  cos ( 4 * PI * i / n);
     
122.             w[k + 2] = (float) -sin ( 8 * PI * i / n);
     
123.             w[k + 3] = (float)  cos ( 8 * PI * i / n);
     
124.             w[k + 4] = (float) -sin (12 * PI * i / n);
     
125.             w[k + 5] = (float)  cos (12 * PI * i / n);
     
126. 
     
127.             k += 6;
     
128.         }
     
129.     }
     
130. }
     
131. 
     
132. /****************************************************************************/
     
133. /*                                                                          */
     
134. /*              主函数                                                      */
     
135. /*                                                                          */
     
136. /****************************************************************************/
     
137. int main(void)
     
138. {
     
139.         // 产生待测试信号
     
140.         unsigned int i;
     
141.         for (i = 0; i F_TOL)
     
142.                 {
     
143.                         printf("失败 %4d 原始数据 %10.5f 逆变换 %10.5f\r\n", i, InputOrig[i], IFFT_Out[i]);
     
144.                 }
     
145.         }
     
146. 
     
147.         // 断点
     
148.         SW_BREAKPOINT;
     
149. }
     

1. Matlab代码
   
2. fs=1000;                                     % 采样频率
   
3. N=1024;                                      % 数据点数
   
4. n=0:N-1;t=n/fs;                              % 时间序列
   
5. x=5*sin(2*pi*150*t)+15*sin(2*pi*350*t);      % 信号
   
6. subplot(2,1,1),plot(t,x);                    % 绘制时域信号图
   
7. y=fft(x,N);                                  % 对信号进行快速傅里叶变换
   
8. mag=abs(y);mag=mag*2/N;                      % 振幅
   
9. f=n*fs/N;                                    % 频率序列
   
10. subplot(2,1,2),plot(f,mag);                  % 绘制振幅图
    
11. xlabel('频率/Hz');
    
12. ylabel('振幅');title('N=1024');grid on;
    
13. grid on;                                     % 绘制振幅图

这是10年前玩的例程了

可以尝试，变旋转因子，尝试玩玩~~·

我也是刚玩这DSP，有时间一定要尝试尝试

尝试之后，可以玩玩旋转因子