Fft算法.
对于变换长度为N的序列x(n)其傅立叶变换可以表示如下：| | N | nk |
| X(k)=DFT[x(n)] = Σx(n)W | ||
| | n="0" | |
式（１）
其中，W="exp"(-2π/N)。
在下面的源码里,对应函数为 DFT
算法导论 里介绍了一种递归计算fft地方法, 函数为FFT_recursive
主要利用了 DFT(x) = DFT[0](x) + w*DFT[1](x)
继而导出迭代的fft计算方法,就是现在最常用的,函数为FFT_Iter
步骤为:
将输入数组元素反置变换
for 树的每一层:
for 这一层上需要计算的每一对数据:
根据这一对数据,及蝶形公式,计算上一层的数据
另外,代码里写了一个简单的复数类和向量计算函数
1
#include <vector>
2#include <math.h>
3#include <assert.h>
4//#include <windows.h>
5
6using namespace std;
7
8//************************************/
9// Fushu
10
/**************************************/
11class FuShu
12{
13
public:
14 FuShu(double r=0, double i=0): x(r), y(i)
15
{
16
}
17
18 FuShu & operator+= (FuShu const &rhs)
19 {
20 x += rhs.x;
21 y += rhs.y;
22 return *this;
23 }
24
25 FuShu & operator-= (FuShu const &rhs)
26 {
27 x -= rhs.x;
28 y -= rhs.y;
29 return *this;
30 }
31
32 FuShu & operator*= (FuShu const &rhs)
33 {
34 double x2 = x * rhs.x - y * rhs.y;
35 double y2 = x * rhs.y + y * rhs.x;
36 x=x2;
37 y=y2;
38 return *this;
39 }
40
41 //private:
42 double x,y;
43
};
44
45FuShu operator+ (FuShu const & lhs, FuShu const & rhs)
46{
47 FuShu res = lhs;
48 return res+=rhs;
49}
50
51FuShu operator- (FuShu const & lhs, FuShu const & rhs)
52{
53 FuShu res = lhs;
54 return res-=rhs;
55}
56
57FuShu operator* (FuShu const & lhs, FuShu const & rhs)
58{
59 FuShu res = lhs;
60 return res*=rhs;
61}
62
63double fabs(FuShu const & lhs)
64{
65 return fabs(lhs.x) + fabs(lhs.y);
66}
67// bool operator== (FuShu const & lhs, FuShu const & rhs)
68// {
69 // return lhs.x == rhs.x && lhs.y == rhs.y;
70// }
71
72void Print(FuShu const & lhs)
73{
74 if(lhs.y>=0)
75 printf("%g+%gj ", lhs.x, lhs.y);
76 else
77 printf("%g%gj ", lhs.x, lhs.y);
78}
79
80void Print(vector <FuShu> const & lhs)
81{
82 size_t n=lhs.size();
83 for(size_t i=0; i<n; i++)
84 Print(lhs[i]);
85 printf("\n");
86}
87
88
89//************************************/
90// vector
91/**************************************/
92
93template <typename T>
94vector <T> operator+ (vector <T> const & lhs, vector <T> const & rhs)
95{
96 assert(lhs.size() == rhs.size());
97 vector <T> res;
98 size_t n = lhs.size();
99 res.resize(n);
100 for(size_t i=0; i<n; i++)
101 res[i] = lhs[i] + rhs[i];
102 return res;
103}
104
105template <typename T>
106vector <T> operator- (vector <T> const & lhs, vector <T> const & rhs)
107{
108 assert(lhs.size() == rhs.size());
109 vector <T> res;
110 size_t n = lhs.size();
111 res.resize(n);
112 for(size_t i=0; i<n; i++)
113 res[i] = lhs[i] - rhs[i];
114 return res;
115}
116
117
118template <typename T>
119vector <T> operator* (FuShu const & lhs, vector <T> const & rhs)
120{
121 vector <T> res;
122 size_t n = rhs.size();
123 res.resize(n);
124 for(size_t i=0; i<n; i++)
125 res[i] = lhs[i] * rhs[i];
126 return res;
127}
128
129template <typename T>
130double fabs(vector <T> const & lhs)
131{
132 double res(0);
133 size_t n=lhs.size();
134
135 for(size_t i=0; i<n; i++)
136 res += fabs(lhs[i]);
137 return res;
138}
139// template <typename T>
140// bool operator== (FuShu const & lhs, vector <T> const & rhs)
141// {
142 // size_t n = lhs.size();
143 // if(n != rhs.size())
144 // return false;
145
146 // for(size_t i=0; i<n; i++)
147 // if(lhs[i] != rhs[i])
148 // return false;
149
150 // return true;
151// }
152
153template <typename T>
154vector <T> & operator<<(vector <T> & lhs, T const &rhs)
155{
156 lhs.push_back(rhs);
157 return lhs;
158}
159
160/***************************************************/
161// FFT
162/***************************************************/
163vector<FuShu> DFT(vector<FuShu> const &a)
164{
165 size_t n=a.size();
166 if(n==1) return a;
167
168 vector<FuShu> res(n);
169 for(size_t k=0; k<n; k++)
170 {
171 FuShu wnk(cos(6.28*k/n), sin(6.28*k/n));
172 FuShu w = 1;
173 for(size_t j=0; j<n; ++j, w*=wnk)
174 res[k] += a[j] * w;
175 }
176
177 return res;
178}
179
180vector<FuShu> FFT_recursive(vector<FuShu> const &a)
181{
182 size_t n=a.size();
183 if(n==1) return a;
184
185 FuShu wn(cos(6.28/n), sin(6.28/n));
186 FuShu w(1);
187 vector<FuShu> a0, a1;
188 a0.reserve(n/2);
189 a1.reserve(n/2);
190 for(size_t i=0; i<n/2; i++)
191 {
192 a0.push_back(a[i*2]);
193 a1.push_back(a[i*2+1]);
194 }
195
196 vector<FuShu> y0 = FFT_recursive(a0);
197 vector<FuShu> y1 = FFT_recursive(a1);
198
199 vector<FuShu> vy;
200 vy.resize(n);
201 for(size_t k=0; k<n/2; k++)
202 {
203 vy[k] = y0[k] + w * y1[k];
204 vy[k + n/2] = y0[k] - w * y1[k];
205 w = w*wn;
206 }
207
208 return vy;
209}
210
211unsigned int rev(unsigned int num, unsigned int nBits)
212{
213 unsigned int r = 0;
214 for(unsigned int i=0; i<nBits; i++)
215 {
216 if(num&(1<<i))
217 r |= 1<<(nBits-i-1);
218 }
219
220 return r;
221}
222
223vector<FuShu> FFT_Iter(vector<FuShu> const &a, unsigned int nBits)
224{
225 size_t n=a.size();
226 if(n==1) return a;
227
228 vector<FuShu> A;
229 A.reserve(n);
230 for(size_t i=0; i<n; i++)
231 A<<a[rev(i,nBits)];
232
233 size_t m=2;
234 for(size_t s=0; s<nBits; s++, m*=2)
235 {
236 FuShu wm(cos(6.28/m), sin(6.28/m));
237 for(size_t k=0; k<n; k+=m)
238 {
239 FuShu w=1;
240 for(size_t j=0; j<m/2; j++, w*=wm)
241 {
242 FuShu t = w * A[k+j+m/2];
243 FuShu u = A[k+j];
244 A[k+j] = u+t;
245 A[k+j+m/2] = u-t;
246 }
247 }
248 }
249
250 return A;
251}
252
253/***************************************************/
254// Main
255/***************************************************/
256
257void main()
258{
259
260 srand(clock());
261 for(int i=0; i<10; i++)
262 {
263 FuShu a(rand(),rand()), b(rand(),rand());
264
265 vector<FuShu> vA;
266 vA<<a <<b<<a <<b;
267
268 printf("input:\n");
269 Print(vA);
270 printf("FFT_recursive result:\n");
271 vector<FuShu> vB = FFT_recursive(vA);
272 Print(vB);
273 printf("DFT result:\n");
274 vector<FuShu> vBDft = DFT(vA);
275 Print(vBDft);
276
277 printf("FFT_Iter result:\n");
278 vector<FuShu> vBItr = FFT_Iter(vA,2);
279 Print(vBItr);
280
281 printf("diff: %g\n", fabs(vB - vBItr));
282 assert( fabs(vB - vBItr)<1e-1);
283 }
284#if 0
285 for(unsigned int i=0; i<8; i++)
286 printf("%d ", rev(i,3));
287#endif
288}
#include <vector>2
#include <math.h>3
#include <assert.h>4
//#include <windows.h>5
6
using namespace std;7
8
//************************************/9
// Fushu10
/**************************************/
11
class FuShu12
{
13
public:14
FuShu(double r=0, double i=0): x(r), y(i)15
{
16
}17
18
FuShu & operator+= (FuShu const &rhs)19
{20
x += rhs.x;21
y += rhs.y;22
return *this;23
}24
25
FuShu & operator-= (FuShu const &rhs)26
{27
x -= rhs.x;28
y -= rhs.y;29
return *this;30
}31
32
FuShu & operator*= (FuShu const &rhs)33
{34
double x2 = x * rhs.x - y * rhs.y;35
double y2 = x * rhs.y + y * rhs.x;36
x=x2;37
y=y2;38
return *this;39
}40
41
//private:42
double x,y;43
};44
45
FuShu operator+ (FuShu const & lhs, FuShu const & rhs)46
{47
FuShu res = lhs;48
return res+=rhs;49
} 50
51
FuShu operator- (FuShu const & lhs, FuShu const & rhs)52
{53
FuShu res = lhs;54
return res-=rhs;55
} 56
57
FuShu operator* (FuShu const & lhs, FuShu const & rhs)58
{59
FuShu res = lhs;60
return res*=rhs;61
} 62
63
double fabs(FuShu const & lhs)64
{65
return fabs(lhs.x) + fabs(lhs.y);66
}67
// bool operator== (FuShu const & lhs, FuShu const & rhs)68
// {69
// return lhs.x == rhs.x && lhs.y == rhs.y;70
// } 71
72
void Print(FuShu const & lhs)73
{74
if(lhs.y>=0)75
printf("%g+%gj ", lhs.x, lhs.y);76
else77
printf("%g%gj ", lhs.x, lhs.y);78
}79
80
void Print(vector <FuShu> const & lhs)81
{82
size_t n=lhs.size();83
for(size_t i=0; i<n; i++)84
Print(lhs[i]);85
printf("\n");86
}87
88
89
//************************************/90
// vector91
/**************************************/92
93
template <typename T>94
vector <T> operator+ (vector <T> const & lhs, vector <T> const & rhs)95
{96
assert(lhs.size() == rhs.size());97
vector <T> res;98
size_t n = lhs.size();99
res.resize(n);100
for(size_t i=0; i<n; i++)101
res[i] = lhs[i] + rhs[i];102
return res;103
}104
105
template <typename T>106
vector <T> operator- (vector <T> const & lhs, vector <T> const & rhs)107
{108
assert(lhs.size() == rhs.size());109
vector <T> res;110
size_t n = lhs.size();111
res.resize(n);112
for(size_t i=0; i<n; i++)113
res[i] = lhs[i] - rhs[i];114
return res;115
}116
117
118
template <typename T>119
vector <T> operator* (FuShu const & lhs, vector <T> const & rhs)120
{121
vector <T> res;122
size_t n = rhs.size();123
res.resize(n);124
for(size_t i=0; i<n; i++)125
res[i] = lhs[i] * rhs[i];126
return res;127
}128
129
template <typename T>130
double fabs(vector <T> const & lhs)131
{132
double res(0);133
size_t n=lhs.size();134
135
for(size_t i=0; i<n; i++)136
res += fabs(lhs[i]);137
return res;138
}139
// template <typename T>140
// bool operator== (FuShu const & lhs, vector <T> const & rhs)141
// {142
// size_t n = lhs.size();143
// if(n != rhs.size())144
// return false;145
146
// for(size_t i=0; i<n; i++)147
// if(lhs[i] != rhs[i])148
// return false;149
150
// return true;151
// }152
153
template <typename T>154
vector <T> & operator<<(vector <T> & lhs, T const &rhs)155
{156
lhs.push_back(rhs);157
return lhs;158
}159
160
/***************************************************/161
// FFT162
/***************************************************/163
vector<FuShu> DFT(vector<FuShu> const &a)164
{165
size_t n=a.size();166
if(n==1) return a;167
168
vector<FuShu> res(n);169
for(size_t k=0; k<n; k++)170
{171
FuShu wnk(cos(6.28*k/n), sin(6.28*k/n));172
FuShu w = 1;173
for(size_t j=0; j<n; ++j, w*=wnk)174
res[k] += a[j] * w;175
}176
177
return res;178
}179
180
vector<FuShu> FFT_recursive(vector<FuShu> const &a)181
{182
size_t n=a.size();183
if(n==1) return a;184
185
FuShu wn(cos(6.28/n), sin(6.28/n));186
FuShu w(1);187
vector<FuShu> a0, a1;188
a0.reserve(n/2);189
a1.reserve(n/2);190
for(size_t i=0; i<n/2; i++)191
{192
a0.push_back(a[i*2]);193
a1.push_back(a[i*2+1]);194
}195
196
vector<FuShu> y0 = FFT_recursive(a0);197
vector<FuShu> y1 = FFT_recursive(a1);198
199
vector<FuShu> vy;200
vy.resize(n);201
for(size_t k=0; k<n/2; k++)202
{203
vy[k] = y0[k] + w * y1[k];204
vy[k + n/2] = y0[k] - w * y1[k];205
w = w*wn;206
}207
208
return vy;209
}210
211
unsigned int rev(unsigned int num, unsigned int nBits)212
{213
unsigned int r = 0;214
for(unsigned int i=0; i<nBits; i++)215
{216
if(num&(1<<i))217
r |= 1<<(nBits-i-1);218
}219
220
return r;221
}222
223
vector<FuShu> FFT_Iter(vector<FuShu> const &a, unsigned int nBits)224
{225
size_t n=a.size();226
if(n==1) return a;227
228
vector<FuShu> A;229
A.reserve(n);230
for(size_t i=0; i<n; i++)231
A<<a[rev(i,nBits)];232
233
size_t m=2;234
for(size_t s=0; s<nBits; s++, m*=2)235
{236
FuShu wm(cos(6.28/m), sin(6.28/m));237
for(size_t k=0; k<n; k+=m)238
{239
FuShu w=1;240
for(size_t j=0; j<m/2; j++, w*=wm)241
{242
FuShu t = w * A[k+j+m/2];243
FuShu u = A[k+j];244
A[k+j] = u+t;245
A[k+j+m/2] = u-t;246
}247
}248
}249
250
return A;251
}252
253
/***************************************************/254
// Main255
/***************************************************/256
257
void main()258
{259
260
srand(clock());261
for(int i=0; i<10; i++)262
{263
FuShu a(rand(),rand()), b(rand(),rand());264
265
vector<FuShu> vA;266
vA<<a <<b<<a <<b;267
268
printf("input:\n");269
Print(vA);270
printf("FFT_recursive result:\n");271
vector<FuShu> vB = FFT_recursive(vA);272
Print(vB);273
printf("DFT result:\n");274
vector<FuShu> vBDft = DFT(vA);275
Print(vBDft);276
277
printf("FFT_Iter result:\n");278
vector<FuShu> vBItr = FFT_Iter(vA,2);279
Print(vBItr);280
281
printf("diff: %g\n", fabs(vB - vBItr));282
assert( fabs(vB - vBItr)<1e-1);283
}284
#if 0 285
for(unsigned int i=0; i<8; i++)286
printf("%d ", rev(i,3));287
#endif 288
}