1.1 --- a/src/Doc/isac/jrocnik/FFT.thy Mon Sep 16 12:27:20 2013 +0200
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,532 +0,0 @@
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