|
1 (* Title: Fast Fourier Transform |
|
2 Author: Clemens Ballarin <ballarin at in.tum.de>, started 12 April 2005 |
|
3 Maintainer: Clemens Ballarin <ballarin at in.tum.de> |
|
4 *) |
|
5 |
|
6 theory FFT |
|
7 imports Complex_Main |
|
8 begin |
|
9 |
|
10 text {* We formalise a functional implementation of the FFT algorithm |
|
11 over the complex numbers, and its inverse. Both are shown |
|
12 equivalent to the usual definitions |
|
13 of these operations through Vandermonde matrices. They are also |
|
14 shown to be inverse to each other, more precisely, that composition |
|
15 of the inverse and the transformation yield the identity up to a |
|
16 scalar. |
|
17 |
|
18 The presentation closely follows Section 30.2 of Cormen \textit{et |
|
19 al.}, \emph{Introduction to Algorithms}, 2nd edition, MIT Press, |
|
20 2003. *} |
|
21 |
|
22 |
|
23 section {* Preliminaries *} |
|
24 |
|
25 lemma of_nat_cplx: |
|
26 "of_nat n = Complex (of_nat n) 0" |
|
27 by (induct n) (simp_all add: complex_one_def) |
|
28 |
|
29 |
|
30 text {* The following two lemmas are useful for experimenting with the |
|
31 transformations, at a vector length of four. *} |
|
32 |
|
33 lemma Ivl4: |
|
34 "{0..<4::nat} = {0, 1, 2, 3}" |
|
35 proof - |
|
36 have "{0..<4::nat} = {0..<Suc (Suc (Suc (Suc 0)))}" by (simp add: eval_nat_numeral) |
|
37 also have "... = {0, 1, 2, 3}" |
|
38 by (simp add: atLeastLessThanSuc eval_nat_numeral insert_commute) |
|
39 finally show ?thesis . |
|
40 qed |
|
41 |
|
42 lemma Sum4: |
|
43 "(\<Sum>i=0..<4::nat. x i) = x 0 + x 1 + x 2 + x 3" |
|
44 by (simp add: Ivl4 eval_nat_numeral) |
|
45 |
|
46 |
|
47 text {* A number of specialised lemmas for the summation operator, |
|
48 where the index set is the natural numbers *} |
|
49 |
|
50 lemma setsum_add_nat_ivl_singleton: |
|
51 assumes less: "m < (n::nat)" |
|
52 shows "f m + setsum f {m<..<n} = setsum f {m..<n}" |
|
53 proof - |
|
54 have "f m + setsum f {m<..<n} = setsum f ({m} \<union> {m<..<n})" |
|
55 by (simp add: setsum_Un_disjoint ivl_disj_int) |
|
56 also from less have "... = setsum f {m..<n}" |
|
57 by (simp only: ivl_disj_un) |
|
58 finally show ?thesis . |
|
59 qed |
|
60 |
|
61 lemma setsum_add_split_nat_ivl_singleton: |
|
62 assumes less: "m < (n::nat)" |
|
63 and g: "!!i. [| m < i; i < n |] ==> g i = f i" |
|
64 shows "f m + setsum g {m<..<n} = setsum f {m..<n}" |
|
65 using less g |
|
66 by(simp add: setsum_add_nat_ivl_singleton cong: strong_setsum_cong) |
|
67 |
|
68 lemma setsum_add_split_nat_ivl: |
|
69 assumes le: "m <= (k::nat)" "k <= n" |
|
70 and g: "!!i. [| m <= i; i < k |] ==> g i = f i" |
|
71 and h: "!!i. [| k <= i; i < n |] ==> h i = f i" |
|
72 shows "setsum g {m..<k} + setsum h {k..<n} = setsum f {m..<n}" |
|
73 using le g h by (simp add: setsum_add_nat_ivl cong: strong_setsum_cong) |
|
74 |
|
75 lemma ivl_splice_Un: |
|
76 "{0..<2*n::nat} = (op * 2 ` {0..<n}) \<union> ((%i. Suc (2*i)) ` {0..<n})" |
|
77 apply (unfold image_def Bex_def) |
|
78 apply auto |
|
79 apply arith |
|
80 done |
|
81 |
|
82 lemma ivl_splice_Int: |
|
83 "(op * 2 ` {0..<n}) \<inter> ((%i. Suc (2*i)) ` {0..<n}) = {}" |
|
84 by auto arith |
|
85 |
|
86 lemma double_inj_on: |
|
87 "inj_on (%i. 2*i::nat) A" |
|
88 by (simp add: inj_onI) |
|
89 |
|
90 lemma Suc_double_inj_on: |
|
91 "inj_on (%i. Suc (2*i)) A" |
|
92 by (rule inj_onI) simp |
|
93 |
|
94 lemma setsum_splice: |
|
95 "(\<Sum>i::nat = 0..<2*n. f i) = (\<Sum>i = 0..<n. f (2*i)) + (\<Sum>i = 0..<n. f (2*i+1))" |
|
96 proof - |
|
97 have "(\<Sum>i::nat = 0..<2*n. f i) = |
|
98 setsum f (op * 2 ` {0..<n}) + setsum f ((%i. 2*i+1) ` {0..<n})" |
|
99 by (simp add: ivl_splice_Un ivl_splice_Int setsum_Un_disjoint) |
|
100 also have "... = (\<Sum>i = 0..<n. f (2*i)) + (\<Sum>i = 0..<n. f (2*i+1))" |
|
101 by (simp add: setsum_reindex [OF double_inj_on] |
|
102 setsum_reindex [OF Suc_double_inj_on]) |
|
103 finally show ?thesis . |
|
104 qed |
|
105 |
|
106 |
|
107 section {* Complex Roots of Unity *} |
|
108 |
|
109 text {* The function @{term cis} from the complex library returns the |
|
110 point on the unity circle corresponding to the argument angle. It |
|
111 is the base for our definition of @{text root}. The main property, |
|
112 De Moirve's formula is already there in the library. *} |
|
113 |
|
114 definition root :: "nat => complex" where |
|
115 "root n == cis (2*pi/(real (n::nat)))" |
|
116 |
|
117 lemma sin_periodic_pi_diff [simp]: "sin (x - pi) = - sin x" |
|
118 by (simp add: sin_diff) |
|
119 |
|
120 lemma sin_cos_between_zero_two_pi: |
|
121 assumes 0: "0 < x" and pi: "x < 2 * pi" |
|
122 shows "sin x \<noteq> 0 \<or> cos x \<noteq> 1" |
|
123 proof - |
|
124 { assume "0 < x" and "x < pi" |
|
125 then have "sin x \<noteq> 0" by (auto dest: sin_gt_zero_pi) } |
|
126 moreover |
|
127 { assume "x = pi" |
|
128 then have "cos x \<noteq> 1" by simp } |
|
129 moreover |
|
130 { assume pi1: "pi < x" and pi2: "x < 2 * pi" |
|
131 then have "0 < x - pi" and "x - pi < pi" by arith+ |
|
132 then have "sin (x - pi) \<noteq> 0" by (auto dest: sin_gt_zero_pi) |
|
133 with pi1 pi2 have "sin x \<noteq> 0" by simp } |
|
134 ultimately show ?thesis using 0 pi by arith |
|
135 qed |
|
136 |
|
137 |
|
138 subsection {* Basic Lemmas *} |
|
139 |
|
140 lemma root_nonzero: |
|
141 "root n ~= 0" |
|
142 apply (unfold root_def) |
|
143 apply (unfold cis_def) |
|
144 apply auto |
|
145 apply (drule sin_zero_abs_cos_one) |
|
146 apply arith |
|
147 done |
|
148 |
|
149 lemma root_unity: |
|
150 "root n ^ n = 1" |
|
151 apply (unfold root_def) |
|
152 apply (simp add: DeMoivre) |
|
153 apply (simp add: cis_def) |
|
154 done |
|
155 |
|
156 lemma root_cancel: |
|
157 "0 < d ==> root (d * n) ^ (d * k) = root n ^ k" |
|
158 apply (unfold root_def) |
|
159 apply (simp add: DeMoivre) |
|
160 done |
|
161 |
|
162 lemma root_summation: |
|
163 assumes k: "0 < k" "k < n" |
|
164 shows "(\<Sum>i=0..<n. (root n ^ k) ^ i) = 0" |
|
165 proof - |
|
166 from k have real0: "0 < real k * (2 * pi) / real n" |
|
167 by (simp add: zero_less_divide_iff |
|
168 mult_strict_right_mono [where a = 0, simplified]) |
|
169 from k mult_strict_right_mono [where a = "real k" and |
|
170 b = "real n" and c = "2 * pi / real n", simplified] |
|
171 have realk: "real k * (2 * pi) / real n < 2 * pi" |
|
172 by (simp add: zero_less_divide_iff) |
|
173 txt {* Main part of the proof *} |
|
174 have "(\<Sum>i=0..<n. (root n ^ k) ^ i) = |
|
175 ((root n ^ k) ^ n - 1) / (root n ^ k - 1)" |
|
176 apply (rule geometric_sum) |
|
177 apply (unfold root_def) |
|
178 apply (simp add: DeMoivre) |
|
179 using real0 realk sin_cos_between_zero_two_pi |
|
180 apply (auto simp add: cis_def complex_one_def) |
|
181 done |
|
182 also have "... = ((root n ^ n) ^ k - 1) / (root n ^ k - 1)" |
|
183 by (simp add: power_mult [THEN sym] mult_ac) |
|
184 also have "... = 0" |
|
185 by (simp add: root_unity) |
|
186 finally show ?thesis . |
|
187 qed |
|
188 |
|
189 lemma root_summation_inv: |
|
190 assumes k: "0 < k" "k < n" |
|
191 shows "(\<Sum>i=0..<n. ((1 / root n) ^ k) ^ i) = 0" |
|
192 proof - |
|
193 from k have real0: "0 < real k * (2 * pi) / real n" |
|
194 by (simp add: zero_less_divide_iff |
|
195 mult_strict_right_mono [where a = 0, simplified]) |
|
196 from k mult_strict_right_mono [where a = "real k" and |
|
197 b = "real n" and c = "2 * pi / real n", simplified] |
|
198 have realk: "real k * (2 * pi) / real n < 2 * pi" |
|
199 by (simp add: zero_less_divide_iff) |
|
200 txt {* Main part of the proof *} |
|
201 have "(\<Sum>i=0..<n. ((1 / root n) ^ k) ^ i) = |
|
202 (((1 / root n) ^ k) ^ n - 1) / ((1 / root n) ^ k - 1)" |
|
203 apply (rule geometric_sum) |
|
204 apply (simp add: nonzero_inverse_eq_divide [THEN sym] root_nonzero) |
|
205 apply (unfold root_def) |
|
206 apply (simp add: DeMoivre) |
|
207 using real0 realk sin_cos_between_zero_two_pi |
|
208 apply (auto simp add: cis_def complex_one_def) |
|
209 done |
|
210 also have "... = (((1 / root n) ^ n) ^ k - 1) / ((1 / root n) ^ k - 1)" |
|
211 by (simp add: power_mult [THEN sym] mult_ac) |
|
212 also have "... = 0" |
|
213 by (simp add: power_divide root_unity) |
|
214 finally show ?thesis . |
|
215 qed |
|
216 |
|
217 lemma root0 [simp]: |
|
218 "root 0 = 1" |
|
219 by (simp add: root_def cis_def) |
|
220 |
|
221 lemma root1 [simp]: |
|
222 "root 1 = 1" |
|
223 by (simp add: root_def cis_def) |
|
224 |
|
225 lemma root2 [simp]: |
|
226 "root 2 = Complex -1 0" |
|
227 by (simp add: root_def cis_def) |
|
228 |
|
229 lemma root4 [simp]: |
|
230 "root 4 = ii" |
|
231 by (simp add: root_def cis_def) |
|
232 |
|
233 |
|
234 subsection {* Derived Lemmas *} |
|
235 |
|
236 lemma root_cancel1: |
|
237 "root (2 * m) ^ (i * (2 * j)) = root m ^ (i * j)" |
|
238 proof - |
|
239 have "root (2 * m) ^ (i * (2 * j)) = root (2 * m) ^ (2 * (i * j))" |
|
240 by (simp add: mult_ac) |
|
241 also have "... = root m ^ (i * j)" |
|
242 by (simp add: root_cancel) |
|
243 finally show ?thesis . |
|
244 qed |
|
245 |
|
246 lemma root_cancel2: |
|
247 "0 < n ==> root (2 * n) ^ n = - 1" |
|
248 txt {* Note the space between @{text "-"} and @{text "1"}. *} |
|
249 using root_cancel [where n = 2 and k = 1] |
|
250 apply (simp only: mult_ac) |
|
251 apply (simp add: complex_one_def) |
|
252 done |
|
253 |
|
254 |
|
255 section {* Discrete Fourier Transformation *} |
|
256 |
|
257 text {* |
|
258 We define operations @{text DFT} and @{text IDFT} for the discrete |
|
259 Fourier Transform and its inverse. Vectors are simply functions of |
|
260 type @{text "nat => complex"}. *} |
|
261 |
|
262 text {* |
|
263 @{text "DFT n a"} is the transform of vector @{text a} |
|
264 of length @{text n}, @{text IDFT} its inverse. *} |
|
265 |
|
266 definition DFT :: "nat => (nat => complex) => (nat => complex)" where |
|
267 "DFT n a == (%i. \<Sum>j=0..<n. (root n) ^ (i * j) * (a j))" |
|
268 |
|
269 definition IDFT :: "nat => (nat => complex) => (nat => complex)" where |
|
270 "IDFT n a == (%i. (\<Sum>k=0..<n. (a k) / (root n) ^ (i * k)))" |
|
271 |
|
272 schematic_lemma "map (DFT 4 a) [0, 1, 2, 3] = ?x" |
|
273 by(simp add: DFT_def Sum4) |
|
274 |
|
275 text {* Lemmas for the correctness proof. *} |
|
276 |
|
277 lemma DFT_lower: |
|
278 "DFT (2 * m) a i = |
|
279 DFT m (%i. a (2 * i)) i + |
|
280 (root (2 * m)) ^ i * DFT m (%i. a (2 * i + 1)) i" |
|
281 proof (unfold DFT_def) |
|
282 have "(\<Sum>j = 0..<2 * m. root (2 * m) ^ (i * j) * a j) = |
|
283 (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j)) * a (2 * j)) + |
|
284 (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j + 1)) * a (2 * j + 1))" |
|
285 (is "?s = _") |
|
286 by (simp add: setsum_splice) |
|
287 also have "... = (\<Sum>j = 0..<m. root m ^ (i * j) * a (2 * j)) + |
|
288 root (2 * m) ^ i * |
|
289 (\<Sum>j = 0..<m. root m ^ (i * j) * a (2 * j + 1))" |
|
290 (is "_ = ?t") |
|
291 txt {* First pair of sums *} |
|
292 apply (simp add: root_cancel1) |
|
293 txt {* Second pair of sums *} |
|
294 apply (simp add: setsum_right_distrib) |
|
295 apply (simp add: power_add) |
|
296 apply (simp add: root_cancel1) |
|
297 apply (simp add: mult_ac) |
|
298 done |
|
299 finally show "?s = ?t" . |
|
300 qed |
|
301 |
|
302 lemma DFT_upper: |
|
303 assumes mbound: "0 < m" and ibound: "m <= i" |
|
304 shows "DFT (2 * m) a i = |
|
305 DFT m (%i. a (2 * i)) (i - m) - |
|
306 root (2 * m) ^ (i - m) * DFT m (%i. a (2 * i + 1)) (i - m)" |
|
307 proof (unfold DFT_def) |
|
308 have "(\<Sum>j = 0..<2 * m. root (2 * m) ^ (i * j) * a j) = |
|
309 (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j)) * a (2 * j)) + |
|
310 (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j + 1)) * a (2 * j + 1))" |
|
311 (is "?s = _") |
|
312 by (simp add: setsum_splice) |
|
313 also have "... = |
|
314 (\<Sum>j = 0..<m. root m ^ ((i - m) * j) * a (2 * j)) - |
|
315 root (2 * m) ^ (i - m) * |
|
316 (\<Sum>j = 0..<m. root m ^ ((i - m) * j) * a (2 * j + 1))" |
|
317 (is "_ = ?t") |
|
318 txt {* First pair of sums *} |
|
319 apply (simp add: root_cancel1) |
|
320 apply (simp add: root_unity ibound root_nonzero power_diff power_mult) |
|
321 txt {* Second pair of sums *} |
|
322 apply (simp add: mbound root_cancel2) |
|
323 apply (simp add: setsum_right_distrib) |
|
324 apply (simp add: power_add) |
|
325 apply (simp add: root_cancel1) |
|
326 apply (simp add: power_mult) |
|
327 apply (simp add: mult_ac) |
|
328 done |
|
329 finally show "?s = ?t" . |
|
330 qed |
|
331 |
|
332 lemma IDFT_lower: |
|
333 "IDFT (2 * m) a i = |
|
334 IDFT m (%i. a (2 * i)) i + |
|
335 (1 / root (2 * m)) ^ i * IDFT m (%i. a (2 * i + 1)) i" |
|
336 proof (unfold IDFT_def) |
|
337 have "(\<Sum>j = 0..<2 * m. a j / root (2 * m) ^ (i * j)) = |
|
338 (\<Sum>j = 0..<m. a (2 * j) / root (2 * m) ^ (i * (2 * j))) + |
|
339 (\<Sum>j = 0..<m. a (2 * j + 1) / root (2 * m) ^ (i * (2 * j + 1)))" |
|
340 (is "?s = _") |
|
341 by (simp add: setsum_splice) |
|
342 also have "... = (\<Sum>j = 0..<m. a (2 * j) / root m ^ (i * j)) + |
|
343 (1 / root (2 * m)) ^ i * |
|
344 (\<Sum>j = 0..<m. a (2 * j + 1) / root m ^ (i * j))" |
|
345 (is "_ = ?t") |
|
346 txt {* First pair of sums *} |
|
347 apply (simp add: root_cancel1) |
|
348 txt {* Second pair of sums *} |
|
349 apply (simp add: setsum_right_distrib) |
|
350 apply (simp add: power_add) |
|
351 apply (simp add: nonzero_power_divide root_nonzero) |
|
352 apply (simp add: root_cancel1) |
|
353 done |
|
354 finally show "?s = ?t" . |
|
355 qed |
|
356 |
|
357 lemma IDFT_upper: |
|
358 assumes mbound: "0 < m" and ibound: "m <= i" |
|
359 shows "IDFT (2 * m) a i = |
|
360 IDFT m (%i. a (2 * i)) (i - m) - |
|
361 (1 / root (2 * m)) ^ (i - m) * |
|
362 IDFT m (%i. a (2 * i + 1)) (i - m)" |
|
363 proof (unfold IDFT_def) |
|
364 have "(\<Sum>j = 0..<2 * m. a j / root (2 * m) ^ (i * j)) = |
|
365 (\<Sum>j = 0..<m. a (2 * j) / root (2 * m) ^ (i * (2 * j))) + |
|
366 (\<Sum>j = 0..<m. a (2 * j + 1) / root (2 * m) ^ (i * (2 * j + 1)))" |
|
367 (is "?s = _") |
|
368 by (simp add: setsum_splice) |
|
369 also have "... = |
|
370 (\<Sum>j = 0..<m. a (2 * j) / root m ^ ((i - m) * j)) - |
|
371 (1 / root (2 * m)) ^ (i - m) * |
|
372 (\<Sum>j = 0..<m. a (2 * j + 1) / root m ^ ((i - m) * j))" |
|
373 (is "_ = ?t") |
|
374 txt {* First pair of sums *} |
|
375 apply (simp add: root_cancel1) |
|
376 apply (simp add: root_unity ibound root_nonzero power_diff power_mult) |
|
377 txt {* Second pair of sums *} |
|
378 apply (simp add: nonzero_power_divide root_nonzero) |
|
379 apply (simp add: mbound root_cancel2) |
|
380 apply (simp add: setsum_divide_distrib) |
|
381 apply (simp add: power_add) |
|
382 apply (simp add: root_cancel1) |
|
383 apply (simp add: power_mult) |
|
384 apply (simp add: mult_ac) |
|
385 done |
|
386 finally show "?s = ?t" . |
|
387 qed |
|
388 |
|
389 text {* @{text DFT} und @{text IDFT} are inverses. *} |
|
390 |
|
391 declare divide_divide_eq_right [simp del] |
|
392 divide_divide_eq_left [simp del] |
|
393 |
|
394 lemma power_diff_inverse: |
|
395 assumes nz: "(a::'a::field) ~= 0" |
|
396 shows "m <= n ==> (inverse a) ^ (n-m) = (a^m) / (a^n)" |
|
397 apply (induct n m rule: diff_induct) |
|
398 apply (simp add: nonzero_power_inverse |
|
399 nonzero_inverse_eq_divide [THEN sym] nz) |
|
400 apply simp |
|
401 apply (simp add: nz) |
|
402 done |
|
403 |
|
404 lemma power_diff_rev_if: |
|
405 assumes nz: "(a::'a::field) ~= 0" |
|
406 shows "(a^m) / (a^n) = (if n <= m then a ^ (m-n) else (1/a) ^ (n-m))" |
|
407 proof (cases "n <= m") |
|
408 case True with nz show ?thesis |
|
409 by (simp add: power_diff) |
|
410 next |
|
411 case False with nz show ?thesis |
|
412 by (simp add: power_diff_inverse nonzero_inverse_eq_divide [THEN sym]) |
|
413 qed |
|
414 |
|
415 lemma power_divides_special: |
|
416 "(a::'a::field) ~= 0 ==> |
|
417 a ^ (i * j) / a ^ (k * i) = (a ^ j / a ^ k) ^ i" |
|
418 by (simp add: nonzero_power_divide power_mult [THEN sym] mult_ac) |
|
419 |
|
420 theorem DFT_inverse: |
|
421 assumes i_less: "i < n" |
|
422 shows "IDFT n (DFT n a) i = of_nat n * a i" |
|
423 using [[linarith_split_limit = 0]] |
|
424 apply (unfold DFT_def IDFT_def) |
|
425 apply (simp add: setsum_divide_distrib) |
|
426 apply (subst setsum_commute) |
|
427 apply (simp only: times_divide_eq_left [THEN sym]) |
|
428 apply (simp only: power_divides_special [OF root_nonzero]) |
|
429 apply (simp add: power_diff_rev_if root_nonzero) |
|
430 apply (simp add: setsum_divide_distrib [THEN sym] |
|
431 setsum_left_distrib [THEN sym]) |
|
432 proof - |
|
433 from i_less have i_diff: "!!k. i - k < n" by arith |
|
434 have diff_i: "!!k. k < n ==> k - i < n" by arith |
|
435 |
|
436 let ?sum = "%i j n. setsum (op ^ (if i <= j then root n ^ (j - i) |
|
437 else (1 / root n) ^ (i - j))) {0..<n} * a j" |
|
438 let ?sum1 = "%i j n. setsum (op ^ (root n ^ (j - i))) {0..<n} * a j" |
|
439 let ?sum2 = "%i j n. setsum (op ^ ((1 / root n) ^ (i - j))) {0..<n} * a j" |
|
440 |
|
441 from i_less have "(\<Sum>j = 0..<n. ?sum i j n) = |
|
442 (\<Sum>j = 0..<i. ?sum2 i j n) + (\<Sum>j = i..<n. ?sum1 i j n)" |
|
443 (is "?s = _") |
|
444 by (simp add: root_summation_inv nat_dvd_not_less |
|
445 setsum_add_split_nat_ivl [where f = "%j. ?sum i j n"]) |
|
446 also from i_less i_diff |
|
447 have "... = (\<Sum>j = i..<n. ?sum1 i j n)" |
|
448 by (simp add: root_summation_inv nat_dvd_not_less) |
|
449 also from i_less have "... = |
|
450 (\<Sum>j\<in>{i} \<union> {i<..<n}. ?sum1 i j n)" |
|
451 by (simp only: ivl_disj_un) |
|
452 also have "... = |
|
453 (?sum1 i i n + (\<Sum>j\<in>{i<..<n}. ?sum1 i j n))" |
|
454 by (simp add: setsum_Un_disjoint ivl_disj_int) |
|
455 also from i_less diff_i have "... = ?sum1 i i n" |
|
456 by (simp add: root_summation nat_dvd_not_less) |
|
457 also from i_less have "... = of_nat n * a i" (is "_ = ?t") |
|
458 by (simp add: of_nat_cplx) |
|
459 finally show "?s = ?t" . |
|
460 qed |
|
461 |
|
462 |
|
463 section {* Discrete, Fast Fourier Transformation *} |
|
464 |
|
465 text {* @{text "FFT k a"} is the transform of vector @{text a} |
|
466 of length @{text "2 ^ k"}, @{text IFFT} its inverse. *} |
|
467 |
|
468 primrec FFT :: "nat => (nat => complex) => (nat => complex)" where |
|
469 "FFT 0 a = a" |
|
470 | "FFT (Suc k) a = |
|
471 (let (x, y) = (FFT k (%i. a (2*i)), FFT k (%i. a (2*i+1))) |
|
472 in (%i. if i < 2^k |
|
473 then x i + (root (2 ^ (Suc k))) ^ i * y i |
|
474 else x (i- 2^k) - (root (2 ^ (Suc k))) ^ (i- 2^k) * y (i- 2^k)))" |
|
475 |
|
476 primrec IFFT :: "nat => (nat => complex) => (nat => complex)" where |
|
477 "IFFT 0 a = a" |
|
478 | "IFFT (Suc k) a = |
|
479 (let (x, y) = (IFFT k (%i. a (2*i)), IFFT k (%i. a (2*i+1))) |
|
480 in (%i. if i < 2^k |
|
481 then x i + (1 / root (2 ^ (Suc k))) ^ i * y i |
|
482 else x (i - 2^k) - |
|
483 (1 / root (2 ^ (Suc k))) ^ (i - 2^k) * y (i - 2^k)))" |
|
484 |
|
485 text {* Finally, for vectors of length @{text "2 ^ k"}, |
|
486 @{text DFT} and @{text FFT}, and @{text IDFT} and |
|
487 @{text IFFT} are equivalent. *} |
|
488 |
|
489 theorem DFT_FFT: |
|
490 "!!a i. i < 2 ^ k ==> DFT (2 ^ k) a i = FFT k a i" |
|
491 proof (induct k) |
|
492 case 0 |
|
493 then show ?case by (simp add: DFT_def) |
|
494 next |
|
495 case (Suc k) |
|
496 assume i: "i < 2 ^ Suc k" |
|
497 show ?case proof (cases "i < 2 ^ k") |
|
498 case True |
|
499 then show ?thesis apply simp apply (simp add: DFT_lower) |
|
500 apply (simp add: Suc) done |
|
501 next |
|
502 case False |
|
503 from i have "i - 2 ^ k < 2 ^ k" by simp |
|
504 with False i show ?thesis apply simp apply (simp add: DFT_upper) |
|
505 apply (simp add: Suc) done |
|
506 qed |
|
507 qed |
|
508 |
|
509 theorem IDFT_IFFT: |
|
510 "!!a i. i < 2 ^ k ==> IDFT (2 ^ k) a i = IFFT k a i" |
|
511 proof (induct k) |
|
512 case 0 |
|
513 then show ?case by (simp add: IDFT_def) |
|
514 next |
|
515 case (Suc k) |
|
516 assume i: "i < 2 ^ Suc k" |
|
517 show ?case proof (cases "i < 2 ^ k") |
|
518 case True |
|
519 then show ?thesis apply simp apply (simp add: IDFT_lower) |
|
520 apply (simp add: Suc) done |
|
521 next |
|
522 case False |
|
523 from i have "i - 2 ^ k < 2 ^ k" by simp |
|
524 with False i show ?thesis apply simp apply (simp add: IDFT_upper) |
|
525 apply (simp add: Suc) done |
|
526 qed |
|
527 qed |
|
528 |
|
529 schematic_lemma "map (FFT (Suc (Suc 0)) a) [0, 1, 2, 3] = ?x" |
|
530 by simp |
|
531 |
|
532 end |