C语言数据结构算法之实现快速傅立叶变换
本实例将实现二维快速傅立叶变换，同时也将借此实例学习用c语言实现矩阵的基本操作、复数的基本掾作，复习所学过的动态内存分配、文件操作、结构指针的函数调用等内容。
很久以来，傅立叶变换一直是许多领域，如线性系统、光学、概率论、量子物理、天线、数字图像处理和信号分析等的一个基本分析工具，但是即便使用计算速度惊人的计算机计算离散傅立叶变换所花费的时间也常常是难以接受的，因此导致了快速傅立叶变换（FFT)的产生。
本实例将对一个二维数组进行正、反快速傅立叶变换。正傅立叶变换时dfft()函数先调用fft()按行对数组进行变换，再对结果调用fft()按列进行变换，此时完成了快速傅立叶变换，再调用rdfft()函数进行傅立叶逆变换。如果程序设计正确的话，变换的结果应与原数组相同。
实例代码：
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
|
#include <stdio.h>#include <stdlib.h>#include <math.h>#define PI 3.14159265358979323846 struct COMPLEX{ float re; float im;} cplx , * Hfield , * S , * R , * w; int n , m;int ln , lm; void initiate ();void dfft ();void rdfft ();void showresult (); void fft (int l , int k);int reverse (int t , int k);void W (int l);int loop (int l);void conjugate (); void add (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z);void sub (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z);void mul (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z);struct COMPLEX * Hread(int i , int j);void Hwrite (int i , int j , struct COMPLEX x); void main (){ initiate (); printf("\n原始数据:\n"); showresult(); getchar (); dfft (); printf("\n快速复利叶变换后的结果:\n"); showresult (); getchar (); rdfft (); printf("\n快速复利叶逆变换后的结果:\n"); showresult (); getchar (); free (Hfield);} void initiate (){//程序初始化操作，包括分配内存、读入要处理的数据、进行显示等 FILE * df; df = fopen ("data.txt" , "r"); fscanf (df , "%5d" , &n); fscanf (df , "%5d" , &m); if ((ln = loop (n)) == -1) { printf (" 列数不是2的整数次幂 "); exit (1); } if ((lm = loop (m)) == -1) { printf (" 行数不是2的整数次幂 "); exit (1); } Hfield = (struct COMPLEX *) malloc (n * m * sizeof (cplx)); if (fread (Hfield , sizeof (cplx) , m * n , df) != (unsigned) (m * n)) { if (feof (df)) printf (" Premature end of file "); else printf (" File read error "); } fclose (df);} void dfft (){//进行二维快速复利叶变换 int i , j; int l , k; l = n; k = ln; w = (struct COMPLEX *) calloc (l , sizeof (cplx)); R = (struct COMPLEX *) calloc (l , sizeof (cplx)); S = (struct COMPLEX *) calloc (l , sizeof(cplx)); W (l); for ( i = 0 ; i < m ; i++ ) {//按行进行快速复利叶变换 for (j = 0 ; j < n ; j++) { S[j].re = Hread (i , j)->re; S[j].im = Hread (i , j)->im; } fft(l , k); for (j = 0 ; j < n ; j++) Hwrite (i , j , R[j]); } free (R); free (S); free (w); l = m; k = lm; w = (struct COMPLEX *) calloc (l , sizeof (cplx)); R = (struct COMPLEX *) calloc (l , sizeof (cplx)); S = (struct COMPLEX *) calloc (l , sizeof (cplx)); W (l); for (i = 0 ; i < n ; i++) {//按列进行快速复利叶变换 for(j = 0 ; j < m ; j++) { S[j].re = Hread(j , i)->re; S[j].im = Hread(j , i)->im; } fft(l , k); for (j = 0 ; j < m ; j++) Hwrite (j , i , R[j]); } free (R); free (S); free (w);} void rdfft (){ conjugate (); dfft (); conjugate ();} void showresult (){ int i , j; for (i = 0 ; i < m ; i++) { printf ( " \n第%d行\n " , i); for (j = 0 ; j < n ; j++) { if (j % 4 == 0) printf (" \n "); printf(" (%5.2f,%5.2fi) " , Hread (i , j)->re , Hread (i , j)->im); } }} void fft (int l , int k){ int i , j , s , nv , t; float c; struct COMPLEX mp , r; nv = l; c = (float) l; c = pow (c , 0.5); for (i = 0 ; i < k ; i++) { for (t = 0 ; t < l ; t += nv) { for (j = 0 ; j < nv / 2 ; j++) { s = (t + j) >> (k - i -1); s = reverse(s , k); r.re = S[t + j].re; r.im = S[t + j].im; mul (&w[s] , &S[t + j + nv / 2] , &mp);/////////讲解传递结构指针和结构本身的区别 add (&r , &mp , &S[t + j]); sub (&r , &mp , &S[t + j + nv / 2]); } } nv = nv >> 1; } for (i = 0 ; i < l ; i++) { j = reverse(i , k); R[j].re = S[i].re / c; R[j].im = S[i].im / c; }} int reverse (int t , int k){ int i , x , y; y = 0; for (i = 0 ; i < k ; i++) { x = t & 1; t = t >> 1; y = (y << 1) + x; } return y;} void W (int l){ int i; float c , a; c = (float) l; c = 2 * PI / c; for (i = 0 ; i < l ; i++) { a = (float) i; w[i].re = (float) cos(a * c); w[i].im = -(float) sin(a * c); }} int loop (int l){//检验输入数据是否为2的整数次幂，如果是返回用2进制表示时的位数 int i , m; if (l != 0) { for (i = 1 ; i < 32 ; i++) { m = l >> i; if (m == 0) break; } if (l == (1 << (i - 1))) return i - 1; } return -1;} void conjugate (){//求复数矩阵的共轭矩阵 int i , j; for (i = 0 ; i < m ; i++) { for (j = 0 ; j < n ; j++) { Hread (i , j)->im *= -1; } }} struct COMPLEX * Hread (int i , int j){//按读矩阵方式返回Hfield中指定位置的指针 return (Hfield + i * n + j);} void Hwrite (int i , int j , struct COMPLEX x){//按写矩阵方式将复数结构x写到指定的Hfield位置上 (Hfield + i * n + j)->re = x.re; (Hfield + i * n + j)->im = x.im;} void add (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z){//定义复数加法 z->re = x->re + y->re; z->im = x->im + y->im; } void sub (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z){//定义复数减法 z->re = x->re - y->re; z->im = x->im - y->im;} void mul (struct COMPLEX * x , struct COMPLEX * y , struct COMPLEX * z){//定义复数乘法 z->re = (x->re) * (y->re) - (x->im) * (y->im); z->im = (x->im) * (y->re) + (x->re) * (y->im);} |
感谢阅读，希望能帮助到大家，谢谢大家对本站的支持！
原文链接：http://www.okbase.net/doc/details/484617
