jan@42245: (*  Title:      Fast Fourier Transform
jan@42245:     Author:     Clemens Ballarin <ballarin at in.tum.de>, started 12 April 2005
jan@42245:     Maintainer: Clemens Ballarin <ballarin at in.tum.de>
jan@42245: *)
jan@42245: 
jan@42245: theory FFT
jan@42245: imports Complex_Main
jan@42245: begin
jan@42245: 
jan@42245: text {* We formalise a functional implementation of the FFT algorithm
jan@42245:   over the complex numbers, and its inverse.  Both are shown
jan@42245:   equivalent to the usual definitions
jan@42245:   of these operations through Vandermonde matrices.  They are also
jan@42245:   shown to be inverse to each other, more precisely, that composition
jan@42245:   of the inverse and the transformation yield the identity up to a
jan@42245:   scalar.
jan@42245: 
jan@42245:   The presentation closely follows Section 30.2 of Cormen \textit{et
jan@42245:   al.}, \emph{Introduction to Algorithms}, 2nd edition, MIT Press,
jan@42245:   2003. *}
jan@42245: 
jan@42245: 
jan@42245: section {* Preliminaries *}
jan@42245: 
jan@42245: lemma of_nat_cplx:
jan@42245:   "of_nat n = Complex (of_nat n) 0"
jan@42245:   by (induct n) (simp_all add: complex_one_def)
jan@42245: 
jan@42245: 
jan@42245: text {* The following two lemmas are useful for experimenting with the
jan@42245:   transformations, at a vector length of four. *}
jan@42245: 
jan@42245: lemma Ivl4:
jan@42245:   "{0..<4::nat} = {0, 1, 2, 3}"
jan@42245: proof -
jan@42245:   have "{0..<4::nat} = {0..<Suc (Suc (Suc (Suc 0)))}" by (simp add: eval_nat_numeral)
jan@42245:   also have "... = {0, 1, 2, 3}"
jan@42245:     by (simp add: atLeastLessThanSuc eval_nat_numeral insert_commute)
jan@42245:   finally show ?thesis .
jan@42245: qed
jan@42245: 
jan@42245: lemma Sum4:
jan@42245:   "(\<Sum>i=0..<4::nat. x i) = x 0 + x 1 + x 2 + x 3"
jan@42245:   by (simp add: Ivl4 eval_nat_numeral)
jan@42245: 
jan@42245: 
jan@42245: text {* A number of specialised lemmas for the summation operator,
jan@42245:   where the index set is the natural numbers *}
jan@42245: 
jan@42245: lemma setsum_add_nat_ivl_singleton:
jan@42245:   assumes less: "m < (n::nat)"
jan@42245:   shows "f m + setsum f {m<..<n} = setsum f {m..<n}"
jan@42245: proof -
jan@42245:   have "f m + setsum f {m<..<n} = setsum f ({m} \<union> {m<..<n})"
jan@42245:     by (simp add: setsum_Un_disjoint ivl_disj_int)
jan@42245:   also from less have "... = setsum f {m..<n}"
jan@42245:     by (simp only: ivl_disj_un)
jan@42245:   finally show ?thesis .
jan@42245: qed
jan@42245: 
jan@42245: lemma setsum_add_split_nat_ivl_singleton:
jan@42245:   assumes less: "m < (n::nat)"
jan@42245:     and g: "!!i. [| m < i; i < n |] ==> g i = f i"
jan@42245:   shows "f m + setsum g {m<..<n} = setsum f {m..<n}"
jan@42245:   using less g
jan@42245:   by(simp add: setsum_add_nat_ivl_singleton cong: strong_setsum_cong)
jan@42245: 
jan@42245: lemma setsum_add_split_nat_ivl:
jan@42245:   assumes le: "m <= (k::nat)" "k <= n"
jan@42245:     and g: "!!i. [| m <= i; i < k |] ==> g i = f i"
jan@42245:     and h: "!!i. [| k <= i; i < n |] ==> h i = f i"
jan@42245:   shows "setsum g {m..<k} + setsum h {k..<n} = setsum f {m..<n}"
jan@42245:   using le g h by (simp add: setsum_add_nat_ivl cong: strong_setsum_cong)
jan@42245: 
jan@42245: lemma ivl_splice_Un:
jan@42245:   "{0..<2*n::nat} = (op * 2 ` {0..<n}) \<union> ((%i. Suc (2*i)) ` {0..<n})"
jan@42245:   apply (unfold image_def Bex_def)
jan@42245:   apply auto
jan@42245:   apply arith
jan@42245:   done
jan@42245: 
jan@42245: lemma ivl_splice_Int:
jan@42245:   "(op * 2 ` {0..<n}) \<inter> ((%i. Suc (2*i)) ` {0..<n}) = {}"
jan@42245:   by auto arith
jan@42245: 
jan@42245: lemma double_inj_on:
jan@42245:   "inj_on (%i. 2*i::nat) A"
jan@42245:   by (simp add: inj_onI)
jan@42245: 
jan@42245: lemma Suc_double_inj_on:
jan@42245:   "inj_on (%i. Suc (2*i)) A"
jan@42245:   by (rule inj_onI) simp
jan@42245: 
jan@42245: lemma setsum_splice:
jan@42245:   "(\<Sum>i::nat = 0..<2*n. f i) = (\<Sum>i = 0..<n. f (2*i)) + (\<Sum>i = 0..<n. f (2*i+1))"
jan@42245: proof -
jan@42245:   have "(\<Sum>i::nat = 0..<2*n. f i) =
jan@42245:     setsum f (op * 2 ` {0..<n}) + setsum f ((%i. 2*i+1) ` {0..<n})"
jan@42245:     by (simp add: ivl_splice_Un ivl_splice_Int setsum_Un_disjoint)
jan@42245:   also have "... = (\<Sum>i = 0..<n. f (2*i)) + (\<Sum>i = 0..<n. f (2*i+1))"
jan@42245:     by (simp add: setsum_reindex [OF double_inj_on]
jan@42245:       setsum_reindex [OF Suc_double_inj_on])
jan@42245:   finally show ?thesis .
jan@42245: qed
jan@42245: 
jan@42245: 
jan@42245: section {* Complex Roots of Unity *}
jan@42245: 
jan@42245: text {* The function @{term cis} from the complex library returns the
jan@42245:   point on the unity circle corresponding to the argument angle.  It
jan@42245:   is the base for our definition of @{text root}.  The main property,
jan@42245:   De Moirve's formula is already there in the library. *}
jan@42245: 
jan@42245: definition root :: "nat => complex" where
jan@42245:   "root n == cis (2*pi/(real (n::nat)))"
jan@42245: 
jan@42245: lemma sin_periodic_pi_diff [simp]: "sin (x - pi) = - sin x"
jan@42245:   by (simp add: sin_diff)
jan@42245: 
jan@42245: lemma sin_cos_between_zero_two_pi:
jan@42245:   assumes 0: "0 < x" and pi: "x < 2 * pi"
jan@42245:   shows "sin x \<noteq> 0 \<or> cos x \<noteq> 1"
jan@42245: proof -
jan@42245:   { assume "0 < x" and "x < pi"
jan@42245:     then have "sin x \<noteq> 0" by (auto dest: sin_gt_zero_pi) }
jan@42245:   moreover
jan@42245:   { assume "x = pi"
jan@42245:     then have "cos x \<noteq> 1" by simp }
jan@42245:   moreover
jan@42245:   { assume pi1: "pi < x" and pi2: "x < 2 * pi"
jan@42245:     then have "0 < x - pi" and "x - pi < pi" by arith+
jan@42245:     then have "sin (x - pi) \<noteq> 0" by (auto dest: sin_gt_zero_pi)
jan@42245:     with pi1 pi2 have "sin x \<noteq> 0" by simp }
jan@42245:   ultimately show ?thesis using 0 pi by arith
jan@42245: qed
jan@42245: 
jan@42245: 
jan@42245: subsection {* Basic Lemmas *}
jan@42245: 
jan@42245: lemma root_nonzero:
jan@42245:   "root n ~= 0"
jan@42245:   apply (unfold root_def)
jan@42245:   apply (unfold cis_def)
jan@42245:   apply auto
jan@42245:   apply (drule sin_zero_abs_cos_one)
jan@42245:   apply arith
jan@42245:   done
jan@42245: 
jan@42245: lemma root_unity:
jan@42245:   "root n ^ n = 1"
jan@42245:   apply (unfold root_def)
jan@42245:   apply (simp add: DeMoivre)
jan@42245:   apply (simp add: cis_def)
jan@42245:   done
jan@42245: 
jan@42245: lemma root_cancel:
jan@42245:   "0 < d ==> root (d * n) ^ (d * k) = root n ^ k"
jan@42245:   apply (unfold root_def)
jan@42245:   apply (simp add: DeMoivre)
jan@42245:   done
jan@42245: 
jan@42245: lemma root_summation:
jan@42245:   assumes k: "0 < k" "k < n"
jan@42245:   shows "(\<Sum>i=0..<n. (root n ^ k) ^ i) = 0"
jan@42245: proof -
jan@42245:   from k have real0: "0 < real k * (2 * pi) / real n"
jan@42245:     by (simp add: zero_less_divide_iff
jan@42245:       mult_strict_right_mono [where a = 0, simplified])
jan@42245:   from k mult_strict_right_mono [where a = "real k" and
jan@42245:     b = "real n" and c = "2 * pi / real n", simplified]
jan@42245:   have realk: "real k * (2 * pi) / real n < 2 * pi"
jan@42245:     by (simp add: zero_less_divide_iff)
jan@42245:   txt {* Main part of the proof *}
jan@42245:   have "(\<Sum>i=0..<n. (root n ^ k) ^ i) =
jan@42245:     ((root n ^ k) ^ n - 1) / (root n ^ k - 1)"
jan@42245:     apply (rule geometric_sum)
jan@42245:     apply (unfold root_def)
jan@42245:     apply (simp add: DeMoivre)
jan@42245:     using real0 realk sin_cos_between_zero_two_pi 
jan@42245:     apply (auto simp add: cis_def complex_one_def)
jan@42245:     done
jan@42245:   also have "... = ((root n ^ n) ^ k - 1) / (root n ^ k - 1)"
jan@42245:     by (simp add: power_mult [THEN sym] mult_ac)
jan@42245:   also have "... = 0"
jan@42245:     by (simp add: root_unity)
jan@42245:   finally show ?thesis .
jan@42245: qed
jan@42245: 
jan@42245: lemma root_summation_inv:
jan@42245:   assumes k: "0 < k" "k < n"
jan@42245:   shows "(\<Sum>i=0..<n. ((1 / root n) ^ k) ^ i) = 0"
jan@42245: proof -
jan@42245:   from k have real0: "0 < real k * (2 * pi) / real n"
jan@42245:     by (simp add: zero_less_divide_iff
jan@42245:       mult_strict_right_mono [where a = 0, simplified])
jan@42245:   from k mult_strict_right_mono [where a = "real k" and
jan@42245:     b = "real n" and c = "2 * pi / real n", simplified]
jan@42245:   have realk: "real k * (2 * pi) / real n < 2 * pi"
jan@42245:     by (simp add: zero_less_divide_iff)
jan@42245:   txt {* Main part of the proof *}
jan@42245:   have "(\<Sum>i=0..<n. ((1 / root n) ^ k) ^ i) =
jan@42245:     (((1 / root n) ^ k) ^ n - 1) / ((1 / root n) ^ k - 1)"
jan@42245:     apply (rule geometric_sum)
jan@42245:     apply (simp add: nonzero_inverse_eq_divide [THEN sym] root_nonzero)
jan@42245:     apply (unfold root_def)
jan@42245:     apply (simp add: DeMoivre)
jan@42245:     using real0 realk sin_cos_between_zero_two_pi
jan@42245:     apply (auto simp add: cis_def complex_one_def)
jan@42245:     done
jan@42245:   also have "... = (((1 / root n) ^ n) ^ k - 1) / ((1 / root n) ^ k - 1)"
jan@42245:     by (simp add: power_mult [THEN sym] mult_ac)
jan@42245:   also have "... = 0"
jan@42245:     by (simp add: power_divide root_unity)
jan@42245:   finally show ?thesis .
jan@42245: qed
jan@42245: 
jan@42245: lemma root0 [simp]:
jan@42245:   "root 0 = 1"
jan@42245:   by (simp add: root_def cis_def)
jan@42245: 
jan@42245: lemma root1 [simp]:
jan@42245:   "root 1 = 1"
jan@42245:   by (simp add: root_def cis_def)
jan@42245: 
jan@42245: lemma root2 [simp]:
jan@42245:   "root 2 = Complex -1 0"
jan@42245:   by (simp add: root_def cis_def)
jan@42245: 
jan@42245: lemma root4 [simp]:
jan@42245:   "root 4 = ii"
jan@42245:   by (simp add: root_def cis_def)
jan@42245: 
jan@42245: 
jan@42245: subsection {* Derived Lemmas *}
jan@42245: 
jan@42245: lemma root_cancel1:
jan@42245:   "root (2 * m) ^ (i * (2 * j)) = root m ^ (i * j)"
jan@42245: proof -
jan@42245:   have "root (2 * m) ^ (i * (2 * j)) = root (2 * m) ^ (2 * (i * j))"
jan@42245:     by (simp add: mult_ac)
jan@42245:   also have "... = root m ^ (i * j)"
jan@42245:     by (simp add: root_cancel)
jan@42245:   finally show ?thesis .
jan@42245: qed
jan@42245: 
jan@42245: lemma root_cancel2:
jan@42245:   "0 < n ==> root (2 * n) ^ n = - 1"
jan@42245:   txt {* Note the space between @{text "-"} and @{text "1"}. *}
jan@42245:   using root_cancel [where n = 2 and k = 1]
jan@42245:   apply (simp only: mult_ac)
jan@42245:   apply (simp add: complex_one_def)
jan@42245:   done
jan@42245: 
jan@42245: 
jan@42245: section {* Discrete Fourier Transformation *}
jan@42245: 
jan@42245: text {*
jan@42245:   We define operations  @{text DFT} and @{text IDFT} for the discrete
jan@42245:   Fourier Transform and its inverse.  Vectors are simply functions of
jan@42245:   type @{text "nat => complex"}. *}
jan@42245: 
jan@42245: text {*
jan@42245:   @{text "DFT n a"} is the transform of vector @{text a}
jan@42245:   of length @{text n}, @{text IDFT} its inverse. *}
jan@42245: 
jan@42245: definition DFT :: "nat => (nat => complex) => (nat => complex)" where
jan@42245:   "DFT n a == (%i. \<Sum>j=0..<n. (root n) ^ (i * j) * (a j))"
jan@42245: 
jan@42245: definition IDFT :: "nat => (nat => complex) => (nat => complex)" where
jan@42245:   "IDFT n a == (%i. (\<Sum>k=0..<n. (a k) / (root n) ^ (i * k)))"
jan@42245: 
jan@42245: schematic_lemma "map (DFT 4 a) [0, 1, 2, 3] = ?x"
jan@42245:   by(simp add: DFT_def Sum4)
jan@42245: 
jan@42245: text {* Lemmas for the correctness proof. *}
jan@42245: 
jan@42245: lemma DFT_lower:
jan@42245:   "DFT (2 * m) a i =
jan@42245:   DFT m (%i. a (2 * i)) i +
jan@42245:   (root (2 * m)) ^ i * DFT m (%i. a (2 * i + 1)) i"
jan@42245: proof (unfold DFT_def)
jan@42245:   have "(\<Sum>j = 0..<2 * m. root (2 * m) ^ (i * j) * a j) =
jan@42245:     (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j)) * a (2 * j)) +
jan@42245:     (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j + 1)) * a (2 * j + 1))"
jan@42245:     (is "?s = _")
jan@42245:     by (simp add: setsum_splice)
jan@42245:   also have "... = (\<Sum>j = 0..<m. root m ^ (i * j) * a (2 * j)) +
jan@42245:     root (2 * m) ^ i *
jan@42245:     (\<Sum>j = 0..<m. root m ^ (i * j) * a (2 * j + 1))"
jan@42245:     (is "_ = ?t")
jan@42245:     txt {* First pair of sums *}
jan@42245:     apply (simp add: root_cancel1)
jan@42245:     txt {* Second pair of sums *}
jan@42245:     apply (simp add: setsum_right_distrib)
jan@42245:     apply (simp add: power_add)
jan@42245:     apply (simp add: root_cancel1)
jan@42245:     apply (simp add: mult_ac)
jan@42245:     done
jan@42245:   finally show "?s = ?t" .
jan@42245: qed
jan@42245: 
jan@42245: lemma DFT_upper:
jan@42245:   assumes mbound: "0 < m" and ibound: "m <= i"
jan@42245:   shows "DFT (2 * m) a i =
jan@42245:     DFT m (%i. a (2 * i)) (i - m) -
jan@42245:     root (2 * m) ^ (i - m) * DFT m (%i. a (2 * i + 1)) (i - m)"
jan@42245: proof (unfold DFT_def)
jan@42245:   have "(\<Sum>j = 0..<2 * m. root (2 * m) ^ (i * j) * a j) =
jan@42245:     (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j)) * a (2 * j)) +
jan@42245:     (\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j + 1)) * a (2 * j + 1))"
jan@42245:     (is "?s = _")
jan@42245:     by (simp add: setsum_splice)
jan@42245:   also have "... =
jan@42245:     (\<Sum>j = 0..<m. root m ^ ((i - m) * j) * a (2 * j)) -
jan@42245:     root (2 * m) ^ (i - m) *
jan@42245:     (\<Sum>j = 0..<m. root m ^ ((i - m) * j) * a (2 * j + 1))"
jan@42245:     (is "_ = ?t")
jan@42245:     txt {* First pair of sums *}
jan@42245:     apply (simp add: root_cancel1)
jan@42245:     apply (simp add: root_unity ibound root_nonzero power_diff power_mult)
jan@42245:     txt {* Second pair of sums *}
jan@42245:     apply (simp add: mbound root_cancel2)
jan@42245:     apply (simp add: setsum_right_distrib)
jan@42245:     apply (simp add: power_add)
jan@42245:     apply (simp add: root_cancel1)
jan@42245:     apply (simp add: power_mult)
jan@42245:     apply (simp add: mult_ac)
jan@42245:     done
jan@42245:   finally show "?s = ?t" .
jan@42245: qed
jan@42245: 
jan@42245: lemma IDFT_lower:
jan@42245:   "IDFT (2 * m) a i =
jan@42245:   IDFT m (%i. a (2 * i)) i +
jan@42245:   (1 / root (2 * m)) ^ i * IDFT m (%i. a (2 * i + 1)) i"
jan@42245: proof (unfold IDFT_def)
jan@42245:   have "(\<Sum>j = 0..<2 * m. a j / root (2 * m) ^ (i * j)) =
jan@42245:     (\<Sum>j = 0..<m. a (2 * j) / root (2 * m) ^ (i * (2 * j))) +
jan@42245:     (\<Sum>j = 0..<m. a (2 * j + 1) / root (2 * m) ^ (i * (2 * j + 1)))"
jan@42245:     (is "?s = _")
jan@42245:     by (simp add: setsum_splice)
jan@42245:   also have "... = (\<Sum>j = 0..<m. a (2 * j) / root m ^ (i * j)) +
jan@42245:     (1 / root (2 * m)) ^ i *
jan@42245:     (\<Sum>j = 0..<m. a (2 * j + 1) / root m ^ (i * j))"
jan@42245:     (is "_ = ?t")
jan@42245:     txt {* First pair of sums *}
jan@42245:     apply (simp add: root_cancel1)
jan@42245:     txt {* Second pair of sums *}
jan@42245:     apply (simp add: setsum_right_distrib)
jan@42245:     apply (simp add: power_add)
jan@42245:     apply (simp add: nonzero_power_divide root_nonzero)
jan@42245:     apply (simp add: root_cancel1)
jan@42245:     done
jan@42245:   finally show "?s = ?t" .
jan@42245: qed
jan@42245: 
jan@42245: lemma IDFT_upper:
jan@42245:   assumes mbound: "0 < m" and ibound: "m <= i"
jan@42245:   shows "IDFT (2 * m) a i =
jan@42245:     IDFT m (%i. a (2 * i)) (i - m) -
jan@42245:     (1 / root (2 * m)) ^ (i - m) *
jan@42245:     IDFT m (%i. a (2 * i + 1)) (i - m)"
jan@42245: proof (unfold IDFT_def)
jan@42245:   have "(\<Sum>j = 0..<2 * m. a j / root (2 * m) ^ (i * j)) =
jan@42245:     (\<Sum>j = 0..<m. a (2 * j) / root (2 * m) ^ (i * (2 * j))) +
jan@42245:     (\<Sum>j = 0..<m. a (2 * j + 1) / root (2 * m) ^ (i * (2 * j + 1)))"
jan@42245:     (is "?s = _")
jan@42245:     by (simp add: setsum_splice)
jan@42245:   also have "... =
jan@42245:     (\<Sum>j = 0..<m. a (2 * j) / root m ^ ((i - m) * j)) -
jan@42245:     (1 / root (2 * m)) ^ (i - m) *
jan@42245:     (\<Sum>j = 0..<m. a (2 * j + 1) / root m ^ ((i - m) * j))"
jan@42245:     (is "_ = ?t")
jan@42245:     txt {* First pair of sums *}
jan@42245:     apply (simp add: root_cancel1)
jan@42245:     apply (simp add: root_unity ibound root_nonzero power_diff power_mult)
jan@42245:     txt {* Second pair of sums *}
jan@42245:     apply (simp add: nonzero_power_divide root_nonzero)
jan@42245:     apply (simp add: mbound root_cancel2)
jan@42245:     apply (simp add: setsum_divide_distrib)
jan@42245:     apply (simp add: power_add)
jan@42245:     apply (simp add: root_cancel1)
jan@42245:     apply (simp add: power_mult)
jan@42245:     apply (simp add: mult_ac)
jan@42245:     done
jan@42245:   finally show "?s = ?t" .
jan@42245: qed
jan@42245: 
jan@42245: text {* @{text DFT} und @{text IDFT} are inverses. *}
jan@42245: 
jan@42245: declare divide_divide_eq_right [simp del]
jan@42245:   divide_divide_eq_left [simp del]
jan@42245: 
jan@42245: lemma power_diff_inverse:
jan@42245:   assumes nz: "(a::'a::field) ~= 0"
jan@42245:   shows "m <= n ==> (inverse a) ^ (n-m) = (a^m) / (a^n)"
jan@42245:   apply (induct n m rule: diff_induct)
jan@42245:     apply (simp add: nonzero_power_inverse
jan@42245:       nonzero_inverse_eq_divide [THEN sym] nz)
jan@42245:    apply simp
jan@42245:   apply (simp add: nz)
jan@42245:   done
jan@42245: 
jan@42245: lemma power_diff_rev_if:
jan@42245:   assumes nz: "(a::'a::field) ~= 0"
jan@42245:   shows "(a^m) / (a^n) = (if n <= m then a ^ (m-n) else (1/a) ^ (n-m))"
jan@42245: proof (cases "n <= m")
jan@42245:   case True with nz show ?thesis
jan@42245:     by (simp add: power_diff)
jan@42245: next
jan@42245:   case False with nz show ?thesis
jan@42245:     by (simp add: power_diff_inverse nonzero_inverse_eq_divide [THEN sym])
jan@42245: qed
jan@42245: 
jan@42245: lemma power_divides_special:
jan@42245:   "(a::'a::field) ~= 0 ==>
jan@42245:   a ^ (i * j) / a ^ (k * i) = (a ^ j / a ^ k) ^ i"
jan@42245:   by (simp add: nonzero_power_divide power_mult [THEN sym] mult_ac)
jan@42245: 
jan@42245: theorem DFT_inverse:
jan@42245:   assumes i_less: "i < n"
jan@42245:   shows  "IDFT n (DFT n a) i = of_nat n * a i"
jan@42245:   using [[linarith_split_limit = 0]]
jan@42245:   apply (unfold DFT_def IDFT_def)
jan@42245:   apply (simp add: setsum_divide_distrib)
jan@42245:   apply (subst setsum_commute)
jan@42245:   apply (simp only: times_divide_eq_left [THEN sym])
jan@42245:   apply (simp only: power_divides_special [OF root_nonzero])
jan@42245:   apply (simp add: power_diff_rev_if root_nonzero)
jan@42245:   apply (simp add: setsum_divide_distrib [THEN sym]
jan@42245:     setsum_left_distrib [THEN sym])
jan@42245:   proof -
jan@42245:     from i_less have i_diff: "!!k. i - k < n" by arith
jan@42245:     have diff_i: "!!k. k < n ==> k - i < n" by arith
jan@42245: 
jan@42245:     let ?sum = "%i j n. setsum (op ^ (if i <= j then root n ^ (j - i)
jan@42245:                   else (1 / root n) ^ (i - j))) {0..<n} * a j"
jan@42245:     let ?sum1 = "%i j n. setsum (op ^ (root n ^ (j - i))) {0..<n} * a j"
jan@42245:     let ?sum2 = "%i j n. setsum (op ^ ((1 / root n) ^ (i - j))) {0..<n} * a j"
jan@42245: 
jan@42245:     from i_less have "(\<Sum>j = 0..<n. ?sum i j n) =
jan@42245:       (\<Sum>j = 0..<i. ?sum2 i j n) + (\<Sum>j = i..<n. ?sum1 i j n)"
jan@42245:       (is "?s = _")
jan@42245:       by (simp add: root_summation_inv nat_dvd_not_less
jan@42245:         setsum_add_split_nat_ivl [where f = "%j. ?sum i j n"])
jan@42245:     also from i_less i_diff
jan@42245:     have "... = (\<Sum>j = i..<n. ?sum1 i j n)"
jan@42245:       by (simp add: root_summation_inv nat_dvd_not_less)
jan@42245:     also from i_less have "... =
jan@42245:       (\<Sum>j\<in>{i} \<union> {i<..<n}. ?sum1 i j n)"
jan@42245:       by (simp only: ivl_disj_un)
jan@42245:     also have "... =
jan@42245:       (?sum1 i i n + (\<Sum>j\<in>{i<..<n}. ?sum1 i j n))"
jan@42245:       by (simp add: setsum_Un_disjoint ivl_disj_int)
jan@42245:     also from i_less diff_i have "... = ?sum1 i i n"
jan@42245:       by (simp add: root_summation nat_dvd_not_less)
jan@42245:     also from i_less have "... = of_nat n * a i" (is "_ = ?t")
jan@42245:       by (simp add: of_nat_cplx)
jan@42245:     finally show "?s = ?t" .
jan@42245:   qed
jan@42245: 
jan@42245: 
jan@42245: section {* Discrete, Fast Fourier Transformation *}
jan@42245: 
jan@42245: text {* @{text "FFT k a"} is the transform of vector @{text a}
jan@42245:   of length @{text "2 ^ k"}, @{text IFFT} its inverse. *}
jan@42245: 
jan@42245: primrec FFT :: "nat => (nat => complex) => (nat => complex)" where
jan@42245:   "FFT 0 a = a"
jan@42245: | "FFT (Suc k) a =
jan@42245:      (let (x, y) = (FFT k (%i. a (2*i)), FFT k (%i. a (2*i+1)))
jan@42245:       in (%i. if i < 2^k
jan@42245:             then x i + (root (2 ^ (Suc k))) ^ i * y i
jan@42245:             else x (i- 2^k) - (root (2 ^ (Suc k))) ^ (i- 2^k) * y (i- 2^k)))"
jan@42245: 
jan@42245: primrec IFFT :: "nat => (nat => complex) => (nat => complex)" where
jan@42245:   "IFFT 0 a = a"
jan@42245: | "IFFT (Suc k) a =
jan@42245:      (let (x, y) = (IFFT k (%i. a (2*i)), IFFT k (%i. a (2*i+1)))
jan@42245:       in (%i. if i < 2^k
jan@42245:             then x i + (1 / root (2 ^ (Suc k))) ^ i * y i
jan@42245:             else x (i - 2^k) -
jan@42245:               (1 / root (2 ^ (Suc k))) ^ (i - 2^k) * y (i - 2^k)))"
jan@42245: 
jan@42245: text {* Finally, for vectors of length @{text "2 ^ k"},
jan@42245:   @{text DFT} and @{text FFT}, and @{text IDFT} and
jan@42245:   @{text IFFT} are equivalent. *}
jan@42245: 
jan@42245: theorem DFT_FFT:
jan@42245:   "!!a i. i < 2 ^ k ==> DFT (2 ^ k) a i = FFT k a i"
jan@42245: proof (induct k)
jan@42245:   case 0
jan@42245:   then show ?case by (simp add: DFT_def)
jan@42245: next
jan@42245:   case (Suc k)
jan@42245:   assume i: "i < 2 ^ Suc k"
jan@42245:   show ?case proof (cases "i < 2 ^ k")
jan@42245:     case True
jan@42245:     then show ?thesis apply simp apply (simp add: DFT_lower)
jan@42245:       apply (simp add: Suc) done
jan@42245:   next
jan@42245:     case False
jan@42245:     from i have "i - 2 ^ k < 2 ^ k" by simp
jan@42245:     with False i show ?thesis apply simp apply (simp add: DFT_upper)
jan@42245:       apply (simp add: Suc) done
jan@42245:   qed
jan@42245: qed
jan@42245: 
jan@42245: theorem IDFT_IFFT:
jan@42245:   "!!a i. i < 2 ^ k ==> IDFT (2 ^ k) a i = IFFT k a i"
jan@42245: proof (induct k)
jan@42245:   case 0
jan@42245:   then show ?case by (simp add: IDFT_def)
jan@42245: next
jan@42245:   case (Suc k)
jan@42245:   assume i: "i < 2 ^ Suc k"
jan@42245:   show ?case proof (cases "i < 2 ^ k")
jan@42245:     case True
jan@42245:     then show ?thesis apply simp apply (simp add: IDFT_lower)
jan@42245:       apply (simp add: Suc) done
jan@42245:   next
jan@42245:     case False
jan@42245:     from i have "i - 2 ^ k < 2 ^ k" by simp
jan@42245:     with False i show ?thesis apply simp apply (simp add: IDFT_upper)
jan@42245:       apply (simp add: Suc) done
jan@42245:   qed
jan@42245: qed
jan@42245: 
jan@42245: schematic_lemma "map (FFT (Suc (Suc 0)) a) [0, 1, 2, 3] = ?x"
jan@42245:   by simp
jan@42245: 
jan@42245: end