1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/doc-isac/jrocnik/FFT.thy Tue Sep 17 09:50:52 2013 +0200
1.3 @@ -0,0 +1,532 @@
1.4 +(* Title: Fast Fourier Transform
1.5 + Author: Clemens Ballarin <ballarin at in.tum.de>, started 12 April 2005
1.6 + Maintainer: Clemens Ballarin <ballarin at in.tum.de>
1.7 +*)
1.8 +
1.9 +theory FFT
1.10 +imports Complex_Main
1.11 +begin
1.12 +
1.13 +text {* We formalise a functional implementation of the FFT algorithm
1.14 + over the complex numbers, and its inverse. Both are shown
1.15 + equivalent to the usual definitions
1.16 + of these operations through Vandermonde matrices. They are also
1.17 + shown to be inverse to each other, more precisely, that composition
1.18 + of the inverse and the transformation yield the identity up to a
1.19 + scalar.
1.20 +
1.21 + The presentation closely follows Section 30.2 of Cormen \textit{et
1.22 + al.}, \emph{Introduction to Algorithms}, 2nd edition, MIT Press,
1.23 + 2003. *}
1.24 +
1.25 +
1.26 +section {* Preliminaries *}
1.27 +
1.28 +lemma of_nat_cplx:
1.29 + "of_nat n = Complex (of_nat n) 0"
1.30 + by (induct n) (simp_all add: complex_one_def)
1.31 +
1.32 +
1.33 +text {* The following two lemmas are useful for experimenting with the
1.34 + transformations, at a vector length of four. *}
1.35 +
1.36 +lemma Ivl4:
1.37 + "{0..<4::nat} = {0, 1, 2, 3}"
1.38 +proof -
1.39 + have "{0..<4::nat} = {0..<Suc (Suc (Suc (Suc 0)))}" by (simp add: eval_nat_numeral)
1.40 + also have "... = {0, 1, 2, 3}"
1.41 + by (simp add: atLeastLessThanSuc eval_nat_numeral insert_commute)
1.42 + finally show ?thesis .
1.43 +qed
1.44 +
1.45 +lemma Sum4:
1.46 + "(\<Sum>i=0..<4::nat. x i) = x 0 + x 1 + x 2 + x 3"
1.47 + by (simp add: Ivl4 eval_nat_numeral)
1.48 +
1.49 +
1.50 +text {* A number of specialised lemmas for the summation operator,
1.51 + where the index set is the natural numbers *}
1.52 +
1.53 +lemma setsum_add_nat_ivl_singleton:
1.54 + assumes less: "m < (n::nat)"
1.55 + shows "f m + setsum f {m<..<n} = setsum f {m..<n}"
1.56 +proof -
1.57 + have "f m + setsum f {m<..<n} = setsum f ({m} \<union> {m<..<n})"
1.58 + by (simp add: setsum_Un_disjoint ivl_disj_int)
1.59 + also from less have "... = setsum f {m..<n}"
1.60 + by (simp only: ivl_disj_un)
1.61 + finally show ?thesis .
1.62 +qed
1.63 +
1.64 +lemma setsum_add_split_nat_ivl_singleton:
1.65 + assumes less: "m < (n::nat)"
1.66 + and g: "!!i. [| m < i; i < n |] ==> g i = f i"
1.67 + shows "f m + setsum g {m<..<n} = setsum f {m..<n}"
1.68 + using less g
1.69 + by(simp add: setsum_add_nat_ivl_singleton cong: strong_setsum_cong)
1.70 +
1.71 +lemma setsum_add_split_nat_ivl:
1.72 + assumes le: "m <= (k::nat)" "k <= n"
1.73 + and g: "!!i. [| m <= i; i < k |] ==> g i = f i"
1.74 + and h: "!!i. [| k <= i; i < n |] ==> h i = f i"
1.75 + shows "setsum g {m..<k} + setsum h {k..<n} = setsum f {m..<n}"
1.76 + using le g h by (simp add: setsum_add_nat_ivl cong: strong_setsum_cong)
1.77 +
1.78 +lemma ivl_splice_Un:
1.79 + "{0..<2*n::nat} = (op * 2 ` {0..<n}) \<union> ((%i. Suc (2*i)) ` {0..<n})"
1.80 + apply (unfold image_def Bex_def)
1.81 + apply auto
1.82 + apply arith
1.83 + done
1.84 +
1.85 +lemma ivl_splice_Int:
1.86 + "(op * 2 ` {0..<n}) \<inter> ((%i. Suc (2*i)) ` {0..<n}) = {}"
1.87 + by auto arith
1.88 +
1.89 +lemma double_inj_on:
1.90 + "inj_on (%i. 2*i::nat) A"
1.91 + by (simp add: inj_onI)
1.92 +
1.93 +lemma Suc_double_inj_on:
1.94 + "inj_on (%i. Suc (2*i)) A"
1.95 + by (rule inj_onI) simp
1.96 +
1.97 +lemma setsum_splice:
1.98 + "(\<Sum>i::nat = 0..<2*n. f i) = (\<Sum>i = 0..<n. f (2*i)) + (\<Sum>i = 0..<n. f (2*i+1))"
1.99 +proof -
1.100 + have "(\<Sum>i::nat = 0..<2*n. f i) =
1.101 + setsum f (op * 2 ` {0..<n}) + setsum f ((%i. 2*i+1) ` {0..<n})"
1.102 + by (simp add: ivl_splice_Un ivl_splice_Int setsum_Un_disjoint)
1.103 + also have "... = (\<Sum>i = 0..<n. f (2*i)) + (\<Sum>i = 0..<n. f (2*i+1))"
1.104 + by (simp add: setsum_reindex [OF double_inj_on]
1.105 + setsum_reindex [OF Suc_double_inj_on])
1.106 + finally show ?thesis .
1.107 +qed
1.108 +
1.109 +
1.110 +section {* Complex Roots of Unity *}
1.111 +
1.112 +text {* The function @{term cis} from the complex library returns the
1.113 + point on the unity circle corresponding to the argument angle. It
1.114 + is the base for our definition of @{text root}. The main property,
1.115 + De Moirve's formula is already there in the library. *}
1.116 +
1.117 +definition root :: "nat => complex" where
1.118 + "root n == cis (2*pi/(real (n::nat)))"
1.119 +
1.120 +lemma sin_periodic_pi_diff [simp]: "sin (x - pi) = - sin x"
1.121 + by (simp add: sin_diff)
1.122 +
1.123 +lemma sin_cos_between_zero_two_pi:
1.124 + assumes 0: "0 < x" and pi: "x < 2 * pi"
1.125 + shows "sin x \<noteq> 0 \<or> cos x \<noteq> 1"
1.126 +proof -
1.127 + { assume "0 < x" and "x < pi"
1.128 + then have "sin x \<noteq> 0" by (auto dest: sin_gt_zero_pi) }
1.129 + moreover
1.130 + { assume "x = pi"
1.131 + then have "cos x \<noteq> 1" by simp }
1.132 + moreover
1.133 + { assume pi1: "pi < x" and pi2: "x < 2 * pi"
1.134 + then have "0 < x - pi" and "x - pi < pi" by arith+
1.135 + then have "sin (x - pi) \<noteq> 0" by (auto dest: sin_gt_zero_pi)
1.136 + with pi1 pi2 have "sin x \<noteq> 0" by simp }
1.137 + ultimately show ?thesis using 0 pi by arith
1.138 +qed
1.139 +
1.140 +
1.141 +subsection {* Basic Lemmas *}
1.142 +
1.143 +lemma root_nonzero:
1.144 + "root n ~= 0"
1.145 + apply (unfold root_def)
1.146 + apply (unfold cis_def)
1.147 + apply auto
1.148 + apply (drule sin_zero_abs_cos_one)
1.149 + apply arith
1.150 + done
1.151 +
1.152 +lemma root_unity:
1.153 + "root n ^ n = 1"
1.154 + apply (unfold root_def)
1.155 + apply (simp add: DeMoivre)
1.156 + apply (simp add: cis_def)
1.157 + done
1.158 +
1.159 +lemma root_cancel:
1.160 + "0 < d ==> root (d * n) ^ (d * k) = root n ^ k"
1.161 + apply (unfold root_def)
1.162 + apply (simp add: DeMoivre)
1.163 + done
1.164 +
1.165 +lemma root_summation:
1.166 + assumes k: "0 < k" "k < n"
1.167 + shows "(\<Sum>i=0..<n. (root n ^ k) ^ i) = 0"
1.168 +proof -
1.169 + from k have real0: "0 < real k * (2 * pi) / real n"
1.170 + by (simp add: zero_less_divide_iff
1.171 + mult_strict_right_mono [where a = 0, simplified])
1.172 + from k mult_strict_right_mono [where a = "real k" and
1.173 + b = "real n" and c = "2 * pi / real n", simplified]
1.174 + have realk: "real k * (2 * pi) / real n < 2 * pi"
1.175 + by (simp add: zero_less_divide_iff)
1.176 + txt {* Main part of the proof *}
1.177 + have "(\<Sum>i=0..<n. (root n ^ k) ^ i) =
1.178 + ((root n ^ k) ^ n - 1) / (root n ^ k - 1)"
1.179 + apply (rule geometric_sum)
1.180 + apply (unfold root_def)
1.181 + apply (simp add: DeMoivre)
1.182 + using real0 realk sin_cos_between_zero_two_pi
1.183 + apply (auto simp add: cis_def complex_one_def)
1.184 + done
1.185 + also have "... = ((root n ^ n) ^ k - 1) / (root n ^ k - 1)"
1.186 + by (simp add: power_mult [THEN sym] mult_ac)
1.187 + also have "... = 0"
1.188 + by (simp add: root_unity)
1.189 + finally show ?thesis .
1.190 +qed
1.191 +
1.192 +lemma root_summation_inv:
1.193 + assumes k: "0 < k" "k < n"
1.194 + shows "(\<Sum>i=0..<n. ((1 / root n) ^ k) ^ i) = 0"
1.195 +proof -
1.196 + from k have real0: "0 < real k * (2 * pi) / real n"
1.197 + by (simp add: zero_less_divide_iff
1.198 + mult_strict_right_mono [where a = 0, simplified])
1.199 + from k mult_strict_right_mono [where a = "real k" and
1.200 + b = "real n" and c = "2 * pi / real n", simplified]
1.201 + have realk: "real k * (2 * pi) / real n < 2 * pi"
1.202 + by (simp add: zero_less_divide_iff)
1.203 + txt {* Main part of the proof *}
1.204 + have "(\<Sum>i=0..<n. ((1 / root n) ^ k) ^ i) =
1.205 + (((1 / root n) ^ k) ^ n - 1) / ((1 / root n) ^ k - 1)"
1.206 + apply (rule geometric_sum)
1.207 + apply (simp add: nonzero_inverse_eq_divide [THEN sym] root_nonzero)
1.208 + apply (unfold root_def)
1.209 + apply (simp add: DeMoivre)
1.210 + using real0 realk sin_cos_between_zero_two_pi
1.211 + apply (auto simp add: cis_def complex_one_def)
1.212 + done
1.213 + also have "... = (((1 / root n) ^ n) ^ k - 1) / ((1 / root n) ^ k - 1)"
1.214 + by (simp add: power_mult [THEN sym] mult_ac)
1.215 + also have "... = 0"
1.216 + by (simp add: power_divide root_unity)
1.217 + finally show ?thesis .
1.218 +qed
1.219 +
1.220 +lemma root0 [simp]:
1.221 + "root 0 = 1"
1.222 + by (simp add: root_def cis_def)
1.223 +
1.224 +lemma root1 [simp]:
1.225 + "root 1 = 1"
1.226 + by (simp add: root_def cis_def)
1.227 +
1.228 +lemma root2 [simp]:
1.229 + "root 2 = Complex -1 0"
1.230 + by (simp add: root_def cis_def)
1.231 +
1.232 +lemma root4 [simp]:
1.233 + "root 4 = ii"
1.234 + by (simp add: root_def cis_def)
1.235 +
1.236 +
1.237 +subsection {* Derived Lemmas *}
1.238 +
1.239 +lemma root_cancel1:
1.240 + "root (2 * m) ^ (i * (2 * j)) = root m ^ (i * j)"
1.241 +proof -
1.242 + have "root (2 * m) ^ (i * (2 * j)) = root (2 * m) ^ (2 * (i * j))"
1.243 + by (simp add: mult_ac)
1.244 + also have "... = root m ^ (i * j)"
1.245 + by (simp add: root_cancel)
1.246 + finally show ?thesis .
1.247 +qed
1.248 +
1.249 +lemma root_cancel2:
1.250 + "0 < n ==> root (2 * n) ^ n = - 1"
1.251 + txt {* Note the space between @{text "-"} and @{text "1"}. *}
1.252 + using root_cancel [where n = 2 and k = 1]
1.253 + apply (simp only: mult_ac)
1.254 + apply (simp add: complex_one_def)
1.255 + done
1.256 +
1.257 +
1.258 +section {* Discrete Fourier Transformation *}
1.259 +
1.260 +text {*
1.261 + We define operations @{text DFT} and @{text IDFT} for the discrete
1.262 + Fourier Transform and its inverse. Vectors are simply functions of
1.263 + type @{text "nat => complex"}. *}
1.264 +
1.265 +text {*
1.266 + @{text "DFT n a"} is the transform of vector @{text a}
1.267 + of length @{text n}, @{text IDFT} its inverse. *}
1.268 +
1.269 +definition DFT :: "nat => (nat => complex) => (nat => complex)" where
1.270 + "DFT n a == (%i. \<Sum>j=0..<n. (root n) ^ (i * j) * (a j))"
1.271 +
1.272 +definition IDFT :: "nat => (nat => complex) => (nat => complex)" where
1.273 + "IDFT n a == (%i. (\<Sum>k=0..<n. (a k) / (root n) ^ (i * k)))"
1.274 +
1.275 +schematic_lemma "map (DFT 4 a) [0, 1, 2, 3] = ?x"
1.276 + by(simp add: DFT_def Sum4)
1.277 +
1.278 +text {* Lemmas for the correctness proof. *}
1.279 +
1.280 +lemma DFT_lower:
1.281 + "DFT (2 * m) a i =
1.282 + DFT m (%i. a (2 * i)) i +
1.283 + (root (2 * m)) ^ i * DFT m (%i. a (2 * i + 1)) i"
1.284 +proof (unfold DFT_def)
1.285 + have "(\<Sum>j = 0..<2 * m. root (2 * m) ^ (i * j) * a j) =
1.286 + (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j)) * a (2 * j)) +
1.287 + (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j + 1)) * a (2 * j + 1))"
1.288 + (is "?s = _")
1.289 + by (simp add: setsum_splice)
1.290 + also have "... = (\<Sum>j = 0..<m. root m ^ (i * j) * a (2 * j)) +
1.291 + root (2 * m) ^ i *
1.292 + (\<Sum>j = 0..<m. root m ^ (i * j) * a (2 * j + 1))"
1.293 + (is "_ = ?t")
1.294 + txt {* First pair of sums *}
1.295 + apply (simp add: root_cancel1)
1.296 + txt {* Second pair of sums *}
1.297 + apply (simp add: setsum_right_distrib)
1.298 + apply (simp add: power_add)
1.299 + apply (simp add: root_cancel1)
1.300 + apply (simp add: mult_ac)
1.301 + done
1.302 + finally show "?s = ?t" .
1.303 +qed
1.304 +
1.305 +lemma DFT_upper:
1.306 + assumes mbound: "0 < m" and ibound: "m <= i"
1.307 + shows "DFT (2 * m) a i =
1.308 + DFT m (%i. a (2 * i)) (i - m) -
1.309 + root (2 * m) ^ (i - m) * DFT m (%i. a (2 * i + 1)) (i - m)"
1.310 +proof (unfold DFT_def)
1.311 + have "(\<Sum>j = 0..<2 * m. root (2 * m) ^ (i * j) * a j) =
1.312 + (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j)) * a (2 * j)) +
1.313 + (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j + 1)) * a (2 * j + 1))"
1.314 + (is "?s = _")
1.315 + by (simp add: setsum_splice)
1.316 + also have "... =
1.317 + (\<Sum>j = 0..<m. root m ^ ((i - m) * j) * a (2 * j)) -
1.318 + root (2 * m) ^ (i - m) *
1.319 + (\<Sum>j = 0..<m. root m ^ ((i - m) * j) * a (2 * j + 1))"
1.320 + (is "_ = ?t")
1.321 + txt {* First pair of sums *}
1.322 + apply (simp add: root_cancel1)
1.323 + apply (simp add: root_unity ibound root_nonzero power_diff power_mult)
1.324 + txt {* Second pair of sums *}
1.325 + apply (simp add: mbound root_cancel2)
1.326 + apply (simp add: setsum_right_distrib)
1.327 + apply (simp add: power_add)
1.328 + apply (simp add: root_cancel1)
1.329 + apply (simp add: power_mult)
1.330 + apply (simp add: mult_ac)
1.331 + done
1.332 + finally show "?s = ?t" .
1.333 +qed
1.334 +
1.335 +lemma IDFT_lower:
1.336 + "IDFT (2 * m) a i =
1.337 + IDFT m (%i. a (2 * i)) i +
1.338 + (1 / root (2 * m)) ^ i * IDFT m (%i. a (2 * i + 1)) i"
1.339 +proof (unfold IDFT_def)
1.340 + have "(\<Sum>j = 0..<2 * m. a j / root (2 * m) ^ (i * j)) =
1.341 + (\<Sum>j = 0..<m. a (2 * j) / root (2 * m) ^ (i * (2 * j))) +
1.342 + (\<Sum>j = 0..<m. a (2 * j + 1) / root (2 * m) ^ (i * (2 * j + 1)))"
1.343 + (is "?s = _")
1.344 + by (simp add: setsum_splice)
1.345 + also have "... = (\<Sum>j = 0..<m. a (2 * j) / root m ^ (i * j)) +
1.346 + (1 / root (2 * m)) ^ i *
1.347 + (\<Sum>j = 0..<m. a (2 * j + 1) / root m ^ (i * j))"
1.348 + (is "_ = ?t")
1.349 + txt {* First pair of sums *}
1.350 + apply (simp add: root_cancel1)
1.351 + txt {* Second pair of sums *}
1.352 + apply (simp add: setsum_right_distrib)
1.353 + apply (simp add: power_add)
1.354 + apply (simp add: nonzero_power_divide root_nonzero)
1.355 + apply (simp add: root_cancel1)
1.356 + done
1.357 + finally show "?s = ?t" .
1.358 +qed
1.359 +
1.360 +lemma IDFT_upper:
1.361 + assumes mbound: "0 < m" and ibound: "m <= i"
1.362 + shows "IDFT (2 * m) a i =
1.363 + IDFT m (%i. a (2 * i)) (i - m) -
1.364 + (1 / root (2 * m)) ^ (i - m) *
1.365 + IDFT m (%i. a (2 * i + 1)) (i - m)"
1.366 +proof (unfold IDFT_def)
1.367 + have "(\<Sum>j = 0..<2 * m. a j / root (2 * m) ^ (i * j)) =
1.368 + (\<Sum>j = 0..<m. a (2 * j) / root (2 * m) ^ (i * (2 * j))) +
1.369 + (\<Sum>j = 0..<m. a (2 * j + 1) / root (2 * m) ^ (i * (2 * j + 1)))"
1.370 + (is "?s = _")
1.371 + by (simp add: setsum_splice)
1.372 + also have "... =
1.373 + (\<Sum>j = 0..<m. a (2 * j) / root m ^ ((i - m) * j)) -
1.374 + (1 / root (2 * m)) ^ (i - m) *
1.375 + (\<Sum>j = 0..<m. a (2 * j + 1) / root m ^ ((i - m) * j))"
1.376 + (is "_ = ?t")
1.377 + txt {* First pair of sums *}
1.378 + apply (simp add: root_cancel1)
1.379 + apply (simp add: root_unity ibound root_nonzero power_diff power_mult)
1.380 + txt {* Second pair of sums *}
1.381 + apply (simp add: nonzero_power_divide root_nonzero)
1.382 + apply (simp add: mbound root_cancel2)
1.383 + apply (simp add: setsum_divide_distrib)
1.384 + apply (simp add: power_add)
1.385 + apply (simp add: root_cancel1)
1.386 + apply (simp add: power_mult)
1.387 + apply (simp add: mult_ac)
1.388 + done
1.389 + finally show "?s = ?t" .
1.390 +qed
1.391 +
1.392 +text {* @{text DFT} und @{text IDFT} are inverses. *}
1.393 +
1.394 +declare divide_divide_eq_right [simp del]
1.395 + divide_divide_eq_left [simp del]
1.396 +
1.397 +lemma power_diff_inverse:
1.398 + assumes nz: "(a::'a::field) ~= 0"
1.399 + shows "m <= n ==> (inverse a) ^ (n-m) = (a^m) / (a^n)"
1.400 + apply (induct n m rule: diff_induct)
1.401 + apply (simp add: nonzero_power_inverse
1.402 + nonzero_inverse_eq_divide [THEN sym] nz)
1.403 + apply simp
1.404 + apply (simp add: nz)
1.405 + done
1.406 +
1.407 +lemma power_diff_rev_if:
1.408 + assumes nz: "(a::'a::field) ~= 0"
1.409 + shows "(a^m) / (a^n) = (if n <= m then a ^ (m-n) else (1/a) ^ (n-m))"
1.410 +proof (cases "n <= m")
1.411 + case True with nz show ?thesis
1.412 + by (simp add: power_diff)
1.413 +next
1.414 + case False with nz show ?thesis
1.415 + by (simp add: power_diff_inverse nonzero_inverse_eq_divide [THEN sym])
1.416 +qed
1.417 +
1.418 +lemma power_divides_special:
1.419 + "(a::'a::field) ~= 0 ==>
1.420 + a ^ (i * j) / a ^ (k * i) = (a ^ j / a ^ k) ^ i"
1.421 + by (simp add: nonzero_power_divide power_mult [THEN sym] mult_ac)
1.422 +
1.423 +theorem DFT_inverse:
1.424 + assumes i_less: "i < n"
1.425 + shows "IDFT n (DFT n a) i = of_nat n * a i"
1.426 + using [[linarith_split_limit = 0]]
1.427 + apply (unfold DFT_def IDFT_def)
1.428 + apply (simp add: setsum_divide_distrib)
1.429 + apply (subst setsum_commute)
1.430 + apply (simp only: times_divide_eq_left [THEN sym])
1.431 + apply (simp only: power_divides_special [OF root_nonzero])
1.432 + apply (simp add: power_diff_rev_if root_nonzero)
1.433 + apply (simp add: setsum_divide_distrib [THEN sym]
1.434 + setsum_left_distrib [THEN sym])
1.435 + proof -
1.436 + from i_less have i_diff: "!!k. i - k < n" by arith
1.437 + have diff_i: "!!k. k < n ==> k - i < n" by arith
1.438 +
1.439 + let ?sum = "%i j n. setsum (op ^ (if i <= j then root n ^ (j - i)
1.440 + else (1 / root n) ^ (i - j))) {0..<n} * a j"
1.441 + let ?sum1 = "%i j n. setsum (op ^ (root n ^ (j - i))) {0..<n} * a j"
1.442 + let ?sum2 = "%i j n. setsum (op ^ ((1 / root n) ^ (i - j))) {0..<n} * a j"
1.443 +
1.444 + from i_less have "(\<Sum>j = 0..<n. ?sum i j n) =
1.445 + (\<Sum>j = 0..<i. ?sum2 i j n) + (\<Sum>j = i..<n. ?sum1 i j n)"
1.446 + (is "?s = _")
1.447 + by (simp add: root_summation_inv nat_dvd_not_less
1.448 + setsum_add_split_nat_ivl [where f = "%j. ?sum i j n"])
1.449 + also from i_less i_diff
1.450 + have "... = (\<Sum>j = i..<n. ?sum1 i j n)"
1.451 + by (simp add: root_summation_inv nat_dvd_not_less)
1.452 + also from i_less have "... =
1.453 + (\<Sum>j\<in>{i} \<union> {i<..<n}. ?sum1 i j n)"
1.454 + by (simp only: ivl_disj_un)
1.455 + also have "... =
1.456 + (?sum1 i i n + (\<Sum>j\<in>{i<..<n}. ?sum1 i j n))"
1.457 + by (simp add: setsum_Un_disjoint ivl_disj_int)
1.458 + also from i_less diff_i have "... = ?sum1 i i n"
1.459 + by (simp add: root_summation nat_dvd_not_less)
1.460 + also from i_less have "... = of_nat n * a i" (is "_ = ?t")
1.461 + by (simp add: of_nat_cplx)
1.462 + finally show "?s = ?t" .
1.463 + qed
1.464 +
1.465 +
1.466 +section {* Discrete, Fast Fourier Transformation *}
1.467 +
1.468 +text {* @{text "FFT k a"} is the transform of vector @{text a}
1.469 + of length @{text "2 ^ k"}, @{text IFFT} its inverse. *}
1.470 +
1.471 +primrec FFT :: "nat => (nat => complex) => (nat => complex)" where
1.472 + "FFT 0 a = a"
1.473 +| "FFT (Suc k) a =
1.474 + (let (x, y) = (FFT k (%i. a (2*i)), FFT k (%i. a (2*i+1)))
1.475 + in (%i. if i < 2^k
1.476 + then x i + (root (2 ^ (Suc k))) ^ i * y i
1.477 + else x (i- 2^k) - (root (2 ^ (Suc k))) ^ (i- 2^k) * y (i- 2^k)))"
1.478 +
1.479 +primrec IFFT :: "nat => (nat => complex) => (nat => complex)" where
1.480 + "IFFT 0 a = a"
1.481 +| "IFFT (Suc k) a =
1.482 + (let (x, y) = (IFFT k (%i. a (2*i)), IFFT k (%i. a (2*i+1)))
1.483 + in (%i. if i < 2^k
1.484 + then x i + (1 / root (2 ^ (Suc k))) ^ i * y i
1.485 + else x (i - 2^k) -
1.486 + (1 / root (2 ^ (Suc k))) ^ (i - 2^k) * y (i - 2^k)))"
1.487 +
1.488 +text {* Finally, for vectors of length @{text "2 ^ k"},
1.489 + @{text DFT} and @{text FFT}, and @{text IDFT} and
1.490 + @{text IFFT} are equivalent. *}
1.491 +
1.492 +theorem DFT_FFT:
1.493 + "!!a i. i < 2 ^ k ==> DFT (2 ^ k) a i = FFT k a i"
1.494 +proof (induct k)
1.495 + case 0
1.496 + then show ?case by (simp add: DFT_def)
1.497 +next
1.498 + case (Suc k)
1.499 + assume i: "i < 2 ^ Suc k"
1.500 + show ?case proof (cases "i < 2 ^ k")
1.501 + case True
1.502 + then show ?thesis apply simp apply (simp add: DFT_lower)
1.503 + apply (simp add: Suc) done
1.504 + next
1.505 + case False
1.506 + from i have "i - 2 ^ k < 2 ^ k" by simp
1.507 + with False i show ?thesis apply simp apply (simp add: DFT_upper)
1.508 + apply (simp add: Suc) done
1.509 + qed
1.510 +qed
1.511 +
1.512 +theorem IDFT_IFFT:
1.513 + "!!a i. i < 2 ^ k ==> IDFT (2 ^ k) a i = IFFT k a i"
1.514 +proof (induct k)
1.515 + case 0
1.516 + then show ?case by (simp add: IDFT_def)
1.517 +next
1.518 + case (Suc k)
1.519 + assume i: "i < 2 ^ Suc k"
1.520 + show ?case proof (cases "i < 2 ^ k")
1.521 + case True
1.522 + then show ?thesis apply simp apply (simp add: IDFT_lower)
1.523 + apply (simp add: Suc) done
1.524 + next
1.525 + case False
1.526 + from i have "i - 2 ^ k < 2 ^ k" by simp
1.527 + with False i show ?thesis apply simp apply (simp add: IDFT_upper)
1.528 + apply (simp add: Suc) done
1.529 + qed
1.530 +qed
1.531 +
1.532 +schematic_lemma "map (FFT (Suc (Suc 0)) a) [0, 1, 2, 3] = ?x"
1.533 + by simp
1.534 +
1.535 +end