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>
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
18 The presentation closely follows Section 30.2 of Cormen \textit{et
19 al.}, \emph{Introduction to Algorithms}, 2nd edition, MIT Press,
23 section {* Preliminaries *}
26 "of_nat n = Complex (of_nat n) 0"
27 by (induct n) (simp_all add: complex_one_def)
30 text {* The following two lemmas are useful for experimenting with the
31 transformations, at a vector length of four. *}
34 "{0..<4::nat} = {0, 1, 2, 3}"
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 .
43 "(\<Sum>i=0..<4::nat. x i) = x 0 + x 1 + x 2 + x 3"
44 by (simp add: Ivl4 eval_nat_numeral)
47 text {* A number of specialised lemmas for the summation operator,
48 where the index set is the natural numbers *}
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}"
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 .
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}"
66 by(simp add: setsum_add_nat_ivl_singleton cong: strong_setsum_cong)
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)
76 "{0..<2*n::nat} = (op * 2 ` {0..<n}) \<union> ((%i. Suc (2*i)) ` {0..<n})"
77 apply (unfold image_def Bex_def)
83 "(op * 2 ` {0..<n}) \<inter> ((%i. Suc (2*i)) ` {0..<n}) = {}"
87 "inj_on (%i. 2*i::nat) A"
88 by (simp add: inj_onI)
90 lemma Suc_double_inj_on:
91 "inj_on (%i. Suc (2*i)) A"
92 by (rule inj_onI) simp
95 "(\<Sum>i::nat = 0..<2*n. f i) = (\<Sum>i = 0..<n. f (2*i)) + (\<Sum>i = 0..<n. f (2*i+1))"
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 .
107 section {* Complex Roots of Unity *}
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. *}
114 definition root :: "nat => complex" where
115 "root n == cis (2*pi/(real (n::nat)))"
117 lemma sin_periodic_pi_diff [simp]: "sin (x - pi) = - sin x"
118 by (simp add: sin_diff)
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"
124 { assume "0 < x" and "x < pi"
125 then have "sin x \<noteq> 0" by (auto dest: sin_gt_zero_pi) }
128 then have "cos x \<noteq> 1" by simp }
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
138 subsection {* Basic Lemmas *}
142 apply (unfold root_def)
143 apply (unfold cis_def)
145 apply (drule sin_zero_abs_cos_one)
151 apply (unfold root_def)
152 apply (simp add: DeMoivre)
153 apply (simp add: cis_def)
157 "0 < d ==> root (d * n) ^ (d * k) = root n ^ k"
158 apply (unfold root_def)
159 apply (simp add: DeMoivre)
162 lemma root_summation:
163 assumes k: "0 < k" "k < n"
164 shows "(\<Sum>i=0..<n. (root n ^ k) ^ i) = 0"
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)
182 also have "... = ((root n ^ n) ^ k - 1) / (root n ^ k - 1)"
183 by (simp add: power_mult [THEN sym] mult_ac)
185 by (simp add: root_unity)
186 finally show ?thesis .
189 lemma root_summation_inv:
190 assumes k: "0 < k" "k < n"
191 shows "(\<Sum>i=0..<n. ((1 / root n) ^ k) ^ i) = 0"
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)
210 also have "... = (((1 / root n) ^ n) ^ k - 1) / ((1 / root n) ^ k - 1)"
211 by (simp add: power_mult [THEN sym] mult_ac)
213 by (simp add: power_divide root_unity)
214 finally show ?thesis .
219 by (simp add: root_def cis_def)
223 by (simp add: root_def cis_def)
226 "root 2 = Complex -1 0"
227 by (simp add: root_def cis_def)
231 by (simp add: root_def cis_def)
234 subsection {* Derived Lemmas *}
237 "root (2 * m) ^ (i * (2 * j)) = root m ^ (i * j)"
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 .
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)
255 section {* Discrete Fourier Transformation *}
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"}. *}
263 @{text "DFT n a"} is the transform of vector @{text a}
264 of length @{text n}, @{text IDFT} its inverse. *}
266 definition DFT :: "nat => (nat => complex) => (nat => complex)" where
267 "DFT n a == (%i. \<Sum>j=0..<n. (root n) ^ (i * j) * (a j))"
269 definition IDFT :: "nat => (nat => complex) => (nat => complex)" where
270 "IDFT n a == (%i. (\<Sum>k=0..<n. (a k) / (root n) ^ (i * k)))"
272 schematic_lemma "map (DFT 4 a) [0, 1, 2, 3] = ?x"
273 by(simp add: DFT_def Sum4)
275 text {* Lemmas for the correctness proof. *}
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))"
286 by (simp add: setsum_splice)
287 also have "... = (\<Sum>j = 0..<m. root m ^ (i * j) * a (2 * j)) +
289 (\<Sum>j = 0..<m. root m ^ (i * j) * a (2 * j + 1))"
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)
299 finally show "?s = ?t" .
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))"
312 by (simp add: setsum_splice)
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))"
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)
329 finally show "?s = ?t" .
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)))"
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))"
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)
354 finally show "?s = ?t" .
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)))"
368 by (simp add: setsum_splice)
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))"
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)
386 finally show "?s = ?t" .
389 text {* @{text DFT} und @{text IDFT} are inverses. *}
391 declare divide_divide_eq_right [simp del]
392 divide_divide_eq_left [simp del]
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)
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)
411 case False with nz show ?thesis
412 by (simp add: power_diff_inverse nonzero_inverse_eq_divide [THEN sym])
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)
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])
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
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"
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)"
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)
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" .
463 section {* Discrete, Fast Fourier Transformation *}
465 text {* @{text "FFT k a"} is the transform of vector @{text a}
466 of length @{text "2 ^ k"}, @{text IFFT} its inverse. *}
468 primrec FFT :: "nat => (nat => complex) => (nat => complex)" where
471 (let (x, y) = (FFT k (%i. a (2*i)), FFT k (%i. a (2*i+1)))
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)))"
476 primrec IFFT :: "nat => (nat => complex) => (nat => complex)" where
479 (let (x, y) = (IFFT k (%i. a (2*i)), IFFT k (%i. a (2*i+1)))
481 then x i + (1 / root (2 ^ (Suc k))) ^ i * y i
483 (1 / root (2 ^ (Suc k))) ^ (i - 2^k) * y (i - 2^k)))"
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. *}
490 "!!a i. i < 2 ^ k ==> DFT (2 ^ k) a i = FFT k a i"
493 then show ?case by (simp add: DFT_def)
496 assume i: "i < 2 ^ Suc k"
497 show ?case proof (cases "i < 2 ^ k")
499 then show ?thesis apply simp apply (simp add: DFT_lower)
500 apply (simp add: Suc) done
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
510 "!!a i. i < 2 ^ k ==> IDFT (2 ^ k) a i = IFFT k a i"
513 then show ?case by (simp add: IDFT_def)
516 assume i: "i < 2 ^ Suc k"
517 show ?case proof (cases "i < 2 ^ k")
519 then show ?thesis apply simp apply (simp add: IDFT_lower)
520 apply (simp add: Suc) done
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
529 schematic_lemma "map (FFT (Suc (Suc 0)) a) [0, 1, 2, 3] = ?x"