Fft算法.
对于变换长度为N的序列x(n)其傅立叶变换可以表示如下: | N | nk |
X(k)=DFT[x(n)] = Σx(n)W | ||
| n="0" | |
式(1)
其中,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
6
using namespace std;
7
8
//************************************/
9
// Fushu
10
/**************************************/
11
class 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
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
else
77
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
// vector
91
/**************************************/
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
// FFT
162
/***************************************************/
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
// Main
255
/***************************************************/
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
}

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

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288
