1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/doc-src/isac/jrocnik/FFT.thy Tue Sep 06 15:57:25 2011 +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
2.1 --- a/doc-src/isac/jrocnik/Test_Z_Transform.thy Mon Sep 05 14:43:38 2011 +0200
2.2 +++ b/doc-src/isac/jrocnik/Test_Z_Transform.thy Tue Sep 06 15:57:25 2011 +0200
2.3 @@ -132,7 +132,8 @@
2.4 term2str s_1;
2.5 term2str s_2;
2.6 *}
2.7 -ML {*
2.8 +
2.9 +ML {* (*Solutions as Denominator --> Denominator1 = z - Zeropoint1, Denominator2 = z-Zeropoint2,...*)
2.10 val xx = HOLogic.dest_eq s_1;
2.11 val s_1' = HOLogic.mk_binop "Groups.minus_class.minus" xx;
2.12 val xx = HOLogic.dest_eq s_2;
2.13 @@ -144,17 +145,22 @@
2.14 subsubsection {*build expression*}
2.15 text {*in isac's CTP-based programming language: let s_1 = Take numerator / (s_1 * s_2)*}
2.16 ML {*
2.17 +(*The Main Denominator is the multiplikation of the partial fraction denominators*)
2.18 val denominator' = HOLogic.mk_binop "Groups.times_class.times" (s_1', s_2') ;
2.19 -val SOME n3 = parseNEW ctxt "3::real";
2.20 -val expression' = HOLogic.mk_binop "Rings.inverse_class.divide" (n3, denominator');
2.21 +val SOME numerator = parseNEW ctxt "3::real";
2.22 +
2.23 +val expression' = HOLogic.mk_binop "Rings.inverse_class.divide" (numerator, denominator');
2.24 term2str expression';
2.25 *}
2.26
2.27 -subsubsection {*Ansatz*}
2.28 +subsubsection {*Ansatz - create partial fractions out of our expression*}
2.29 +
2.30 axiomatization where
2.31 ansatz2: "n / (a*b) = A/a + B/(b::real)" and
2.32 multiply_eq2: "(n / (a*b) = A/a + B/b) = (a*b*(n / (a*b)) = a*b*(A/a + B/b))"
2.33 +
2.34 ML {*
2.35 +(*we use our ansatz2 to rewrite our expression and get an equilation with our expression on the left and the partial fractions of it on the right side*)
2.36 val SOME (t1,_) = rewrite_ @{theory Isac} e_rew_ord e_rls false @{thm ansatz2} expression';
2.37 term2str t1;
2.38 atomty t1;
2.39 @@ -162,10 +168,12 @@
2.40 term2str eq1;
2.41 *}
2.42 ML {*
2.43 +(*eliminate the demoninators by multiplying the left and the right side with the main denominator*)
2.44 val SOME (eq2,_) = rewrite_ @{theory Isac} e_rew_ord e_rls false @{thm multiply_eq2} eq1;
2.45 term2str eq2;
2.46 *}
2.47 ML {*
2.48 +(*simplificatoin*)
2.49 val SOME (eq3,_) = rewrite_set_ @{theory Isac} false norm_Rational eq2;
2.50 term2str eq3; (*?A ?B not simplified*)
2.51 *}
2.52 @@ -189,6 +197,7 @@
2.53 subsubsection {*get first koeffizient*}
2.54
2.55 ML {*
2.56 +(*substitude z with the first zeropoint to get A*)
2.57 val SOME (eq4_1,_) = rewrite_terms_ @{theory Isac} e_rew_ord e_rls [s_1] eq3'';
2.58 term2str eq4_1;
2.59 *}
2.60 @@ -202,6 +211,7 @@
2.61
2.62 *}
2.63 ML {*
2.64 +(*solve the simple linear equilation for A*)
2.65 val (p,_,f,nxt,_,pt) = CalcTreeTEST [(fmz, (dI',pI',mI'))];
2.66 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
2.67 val (p,_,f,nxt,_,pt) = me nxt p [] pt;
2.68 @@ -238,6 +248,7 @@
2.69 subsubsection {*get second koeffizient*}
2.70
2.71 ML {*
2.72 +(*substitude z with the second zeropoint to get B*)
2.73 val SOME (eq4b_1,_) = rewrite_terms_ @{theory Isac} e_rew_ord e_rls [s_2] eq3'';
2.74 term2str eq4_1;
2.75 *}
2.76 @@ -248,6 +259,8 @@
2.77 *}
2.78
2.79 ML {*
2.80 +(*solve the simple linear equilation for B*)
2.81 +
2.82 val fmz = ["equality (3 = -3 * B / (4::real))", "solveFor B","solutions L"];
2.83 val (dI',pI',mI') =("Isac", ["univariate","equation"], ["no_met"]);
2.84
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
3.2 +++ b/doc-src/isac/jrocnik/bakkarbeit_jrocnik.tex Tue Sep 06 15:57:25 2011 +0200
3.3 @@ -0,0 +1,186 @@
3.4 +\documentclass[a4paper,12pt]{article}
3.5 +
3.6 +\usepackage[german]{babel}
3.7 +\usepackage[T1]{fontenc}
3.8 +\usepackage[latin1]{inputenc}
3.9 +
3.10 +\usepackage{graphicx}
3.11 +
3.12 +\bibliographystyle{alpha}
3.13 +
3.14 +\def\isac{${\cal I}\mkern-2mu{\cal S}\mkern-5mu{\cal AC}$}
3.15 +\def\sisac{\footnotesize${\cal I}\mkern-2mu{\cal S}\mkern-5mu{\cal AC}$}
3.16 +
3.17 +\begin{document}
3.18 +
3.19 +\title{
3.20 + \Large{
3.21 + \bf Interactive Course Material for Signal Processing based on Isabelle/\isac\\~\\
3.22 + }
3.23 + \sisac-Projektteam des Instituts für Softwaretechnologie,\\Technische Universität Graz\\
3.24 + \vspace{0.7cm}
3.25 + \large{
3.26 + Betreuer: Dr. Walther Neuper
3.27 + }
3.28 +}
3.29 +\author{Jan Simon Rocnik\\{\tt jan.rocnik@student.tugraz.at}}
3.30 +
3.31 +\date{\today}
3.32 +\maketitle
3.33 +\clearpage
3.34 +\tableofcontents
3.35 +\clearpage
3.36 +
3.37 +
3.38 +\section{Introduction}
3.39 +
3.40 +todo
3.41 +
3.42 +motivation from \textbf{practice of mathematics learning} ... STEOP
3.43 +
3.44 +\textbf{mathematics applied} in signal processing (SP)
3.45 +
3.46 +mathematics mechanized in Computer Theorem Provers \textbf{(CTP)} ... (almost) traditional mathematical notation (predicate calculus) for axioms, definitions, lemmas, theorems. Recent developments provide also proofs in a humand readable format \cite{TODO}
3.47 +
3.48 +This thesis tries to \textbf{connect these two worlds} ... this trial is one of the first; others see related work
3.49 +
3.50 +\subsection{Mechanization of Mathematics}
3.51 +
3.52 +todo
3.53 +
3.54 +hughe theories of mathematics
3.55 +
3.56 +still a hugh gap between these theories and ``real applications'' e.g. in SP
3.57 +
3.58 +? academic engineering starts from physics (experimentation, measurement) and then proceeds to mathematical modelling --- mechanized math starts from mathematical models and (hopefully !) proceeds to match physics.
3.59 +
3.60 +CTP Isabelle ... survey of knowledge, links to knowledge
3.61 +
3.62 +\paragraph{\sisac}
3.63 +TODO
3.64 +
3.65 +adds an ``application'' axis (formal specifications of problems) and an ``algorithmic'' axis to the ``deductive'' axis of knowledge ... 3-dimensional universe of mathematics.
3.66 +
3.67 +\subsection{Goals of the Thesis}
3.68 +
3.69 +todo
3.70 +
3.71 +\subsection{Milestones for the Project}
3.72 +Die Planung des Projekts teilt sich in folgende Iterationen:
3.73 +\begin{enumerate}
3.74 +\item \textbf{Sammeln von Informationen über Themengebiete und deren Realisierbarkeit } (29.06. -- 27.07.)
3.75 +identify problems relevant for certain SP lectures
3.76 +
3.77 +estimate chances to realized them within the scope of this thesis
3.78 +
3.79 +order for implementing the problems negotiated with lecturers
3.80 +
3.81 +
3.82 +\item \textbf{1. Präsentation - Auswählen der realisierbaren Themengebiete} (27.07.)
3.83 +\item \textbf{Ausarbeiten der Aufgaben in \isac} (01.09. -- 11.11.)
3.84 +\item \textbf{Dokumentation der Aufgaben} (14.11. -- 02.12.)
3.85 +\item \textbf{Ausarbeitung in Latex, Bakkarbeit} (05.12. -- todo)
3.86 +\item \textbf{2. Präsentation - Abschluss der Arbeit} (todo)
3.87 +\end{enumerate}
3.88 +
3.89 +\subsection{Structure of the Thesis}
3.90 +
3.91 +todo
3.92 +
3.93 +\section{Mechanization of Mathematics for SP Problems}
3.94 +todo
3.95 +
3.96 +\subsection{Relevant Knowledge available in Isabelle}
3.97 +todo
3.98 +
3.99 +\paragraph{example FFT}, describe in detail !!!!
3.100 +
3.101 +? different meaning: FFT in Maple
3.102 +
3.103 +gap between what is available and what is required (@)!
3.104 +
3.105 +traditional notation ?
3.106 +
3.107 +\subsection{Relevant Knowledge available in \isac}
3.108 +todo
3.109 +
3.110 +specifications (``application axis'') and methods (``algorithmic axis'')
3.111 +
3.112 +partial fractions, cancellation of multivariate rational terms, ...
3.113 +
3.114 +\subsection{Survey: Available Knowledge and Selected Problems}
3.115 +todo
3.116 +
3.117 +estimate gap (@) for each problem (tables)
3.118 +
3.119 +conclusion: following order for implementing the problems ...
3.120 +
3.121 +\subsection{Formalization of missing knowledge in Isabelle}
3.122 +todo
3.123 +
3.124 +axiomatization ... where ... and
3.125 +
3.126 +\subsection{Notes on Problems with Traditional Notation}
3.127 +todo
3.128 +
3.129 +u[n] !!
3.130 +
3.131 +f x = why not f(x) ?!?!
3.132 +
3.133 +...
3.134 +
3.135 +\section{Implementation of Certain SP Problems}
3.136 +todo
3.137 +
3.138 +\subsection{Formal Specification of Problems}
3.139 +todo
3.140 +
3.141 +\subsection{Methods Solving the Problems}
3.142 +todo
3.143 +
3.144 +\subsection{Integration of Subproblems available in \isac}
3.145 +todo
3.146 +
3.147 +\subsection{Examples and Multimedia Content}
3.148 +todo
3.149 +
3.150 +
3.151 +\section{Related Work and Open Questions}
3.152 +todo
3.153 +
3.154 +See ``introduction'': This thesis tries to connect these two worlds ... this trial is one of the first; others see related work
3.155 +
3.156 +
3.157 +
3.158 +\section{Beschreibung der Meilensteine}\label{ms-desc}
3.159 +todo
3.160 +\section{Bericht zum Projektverlauf}
3.161 +todo
3.162 +\section{Abschliesende Bemerkungen}
3.163 +todo
3.164 +
3.165 +\clearpage
3.166 +
3.167 +\bibliography{bib}
3.168 +
3.169 +\clearpage
3.170 +
3.171 +\appendix
3.172 +%\section*{Anhang}
3.173 +\section{Demobeispiel}\label{demo-code}
3.174 +\begin{verbatim}
3.175 +
3.176 +bsp
3.177 +
3.178 +\end{verbatim}
3.179 +
3.180 +\section{Stundenliste}
3.181 +
3.182 +\subsection*{Voraussetzungen zum Arbeitsbeginn schaffen}
3.183 +\begin{tabular}[t]{lll}
3.184 + {\bf Datum} & {\bf Stunden} & {\bf Beschreibung} \\
3.185 + 10.02.2011 & 2:00 & Besprechung der Problemstellung \\
3.186 +\end{tabular}
3.187 +
3.188 +
3.189 +\end{document}
4.1 --- a/doc-src/isac/jrocnik/vorlage-bakkarbeit.tex Mon Sep 05 14:43:38 2011 +0200
4.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
4.3 @@ -1,186 +0,0 @@
4.4 -\documentclass[a4paper,12pt]{article}
4.5 -
4.6 -\usepackage[german]{babel}
4.7 -\usepackage[T1]{fontenc}
4.8 -\usepackage[latin1]{inputenc}
4.9 -
4.10 -\usepackage{graphicx}
4.11 -
4.12 -\bibliographystyle{alpha}
4.13 -
4.14 -\def\isac{${\cal I}\mkern-2mu{\cal S}\mkern-5mu{\cal AC}$}
4.15 -\def\sisac{\footnotesize${\cal I}\mkern-2mu{\cal S}\mkern-5mu{\cal AC}$}
4.16 -
4.17 -\begin{document}
4.18 -
4.19 -\title{
4.20 - \Large{
4.21 - \bf Interactive Course Material for Signal Processing based on Isabelle/\isac\\~\\
4.22 - }
4.23 - \sisac-Projektteam des Instituts für Softwaretechnologie,\\Technische Universität Graz\\
4.24 - \vspace{0.7cm}
4.25 - \large{
4.26 - Betreuer: Dr. Walther Neuper
4.27 - }
4.28 -}
4.29 -\author{Jan Simon Rocnik\\{\tt jan.rocnik@student.tugraz.at}}
4.30 -
4.31 -\date{\today}
4.32 -\maketitle
4.33 -\clearpage
4.34 -\tableofcontents
4.35 -\clearpage
4.36 -
4.37 -
4.38 -\section{Introduction}
4.39 -
4.40 -todo
4.41 -
4.42 -motivation from \textbf{practice of mathematics learning} ... STEOP
4.43 -
4.44 -\textbf{mathematics applied} in signal processing (SP)
4.45 -
4.46 -mathematics mechanized in Computer Theorem Provers \textbf{(CTP)} ... (almost) traditional mathematical notation (predicate calculus) for axioms, definitions, lemmas, theorems. Recent developments provide also proofs in a humand readable format \cite{TODO}
4.47 -
4.48 -This thesis tries to \textbf{connect these two worlds} ... this trial is one of the first; others see related work
4.49 -
4.50 -\subsection{Mechanization of Mathematics}
4.51 -
4.52 -todo
4.53 -
4.54 -hughe theories of mathematics
4.55 -
4.56 -still a hugh gap between these theories and ``real applications'' e.g. in SP
4.57 -
4.58 -? academic engineering starts from physics (experimentation, measurement) and then proceeds to mathematical modelling --- mechanized math starts from mathematical models and (hopefully !) proceeds to match physics.
4.59 -
4.60 -CTP Isabelle ... survey of knowledge, links to knowledge
4.61 -
4.62 -\paragraph{\sisac}
4.63 -TODO
4.64 -
4.65 -adds an ``application'' axis (formal specifications of problems) and an ``algorithmic'' axis to the ``deductive'' axis of knowledge ... 3-dimensional universe of mathematics.
4.66 -
4.67 -\subsection{Goals of the Thesis}
4.68 -
4.69 -todo
4.70 -
4.71 -\subsection{Milestones for the Project}
4.72 -Die Planung des Projekts teilt sich in folgende Iterationen:
4.73 -\begin{enumerate}
4.74 -\item \textbf{Sammeln von Informationen über Themengebiete und deren Realisierbarkeit } (29.06. -- 27.07.)
4.75 -identify problems relevant for certain SP lectures
4.76 -
4.77 -estimate chances to realized them within the scope of this thesis
4.78 -
4.79 -order for implementing the problems negotiated with lecturers
4.80 -
4.81 -
4.82 -\item \textbf{1. Präsentation - Auswählen der realisierbaren Themengebiete} (27.07.)
4.83 -\item \textbf{Ausarbeiten der Aufgaben in \isac} (01.09. -- 11.11.)
4.84 -\item \textbf{Dokumentation der Aufgaben} (14.11. -- 02.12.)
4.85 -\item \textbf{Ausarbeitung in Latex, Bakkarbeit} (05.12. -- todo)
4.86 -\item \textbf{2. Präsentation - Abschluss der Arbeit} (todo)
4.87 -\end{enumerate}
4.88 -
4.89 -\subsection{Structure of the Thesis}
4.90 -
4.91 -todo
4.92 -
4.93 -\section{Mechanization of Mathematics for SP Problems}
4.94 -todo
4.95 -
4.96 -\subsection{Relevant Knowledge available in Isabelle}
4.97 -todo
4.98 -
4.99 -\paragraph{example FFT}, describe in detail !!!!
4.100 -
4.101 -? different meaning: FFT in Maple
4.102 -
4.103 -gap between what is available and what is required (@)!
4.104 -
4.105 -traditional notation ?
4.106 -
4.107 -\subsection{Relevant Knowledge available in \isac}
4.108 -todo
4.109 -
4.110 -specifications (``application axis'') and methods (``algorithmic axis'')
4.111 -
4.112 -partial fractions, cancellation of multivariate rational terms, ...
4.113 -
4.114 -\subsection{Survey: Available Knowledge and Selected Problems}
4.115 -todo
4.116 -
4.117 -estimate gap (@) for each problem (tables)
4.118 -
4.119 -conclusion: following order for implementing the problems ...
4.120 -
4.121 -\subsection{Formalization of missing knowledge in Isabelle}
4.122 -todo
4.123 -
4.124 -axiomatization ... where ... and
4.125 -
4.126 -\subsection{Notes on Problems with Traditional Notation}
4.127 -todo
4.128 -
4.129 -u[n] !!
4.130 -
4.131 -f x = why not f(x) ?!?!
4.132 -
4.133 -...
4.134 -
4.135 -\section{Implementation of Certain SP Problems}
4.136 -todo
4.137 -
4.138 -\subsection{Formal Specification of Problems}
4.139 -todo
4.140 -
4.141 -\subsection{Methods Solving the Problems}
4.142 -todo
4.143 -
4.144 -\subsection{Integration of Subproblems available in \isac}
4.145 -todo
4.146 -
4.147 -\subsection{Examples and Multimedia Content}
4.148 -todo
4.149 -
4.150 -
4.151 -\section{Related Work and Open Questions}
4.152 -todo
4.153 -
4.154 -See ``introduction'': This thesis tries to connect these two worlds ... this trial is one of the first; others see related work
4.155 -
4.156 -
4.157 -
4.158 -\section{Beschreibung der Meilensteine}\label{ms-desc}
4.159 -todo
4.160 -\section{Bericht zum Projektverlauf}
4.161 -todo
4.162 -\section{Abschliesende Bemerkungen}
4.163 -todo
4.164 -
4.165 -\clearpage
4.166 -
4.167 -\bibliography{bib}
4.168 -
4.169 -\clearpage
4.170 -
4.171 -\appendix
4.172 -%\section*{Anhang}
4.173 -\section{Demobeispiel}\label{demo-code}
4.174 -\begin{verbatim}
4.175 -
4.176 -bsp
4.177 -
4.178 -\end{verbatim}
4.179 -
4.180 -\section{Stundenliste}
4.181 -
4.182 -\subsection*{Voraussetzungen zum Arbeitsbeginn schaffen}
4.183 -\begin{tabular}[t]{lll}
4.184 - {\bf Datum} & {\bf Stunden} & {\bf Beschreibung} \\
4.185 - 10.02.2011 & 2:00 & Besprechung der Problemstellung \\
4.186 -\end{tabular}
4.187 -
4.188 -
4.189 -\end{document}