jan@42245
|
1 |
(* Title: Fast Fourier Transform
|
jan@42245
|
2 |
Author: Clemens Ballarin <ballarin at in.tum.de>, started 12 April 2005
|
jan@42245
|
3 |
Maintainer: Clemens Ballarin <ballarin at in.tum.de>
|
jan@42245
|
4 |
*)
|
jan@42245
|
5 |
|
jan@42245
|
6 |
theory FFT
|
jan@42245
|
7 |
imports Complex_Main
|
jan@42245
|
8 |
begin
|
jan@42245
|
9 |
|
jan@42245
|
10 |
text {* We formalise a functional implementation of the FFT algorithm
|
jan@42245
|
11 |
over the complex numbers, and its inverse. Both are shown
|
jan@42245
|
12 |
equivalent to the usual definitions
|
jan@42245
|
13 |
of these operations through Vandermonde matrices. They are also
|
jan@42245
|
14 |
shown to be inverse to each other, more precisely, that composition
|
jan@42245
|
15 |
of the inverse and the transformation yield the identity up to a
|
jan@42245
|
16 |
scalar.
|
jan@42245
|
17 |
|
jan@42245
|
18 |
The presentation closely follows Section 30.2 of Cormen \textit{et
|
jan@42245
|
19 |
al.}, \emph{Introduction to Algorithms}, 2nd edition, MIT Press,
|
jan@42245
|
20 |
2003. *}
|
jan@42245
|
21 |
|
jan@42245
|
22 |
|
jan@42245
|
23 |
section {* Preliminaries *}
|
jan@42245
|
24 |
|
jan@42245
|
25 |
lemma of_nat_cplx:
|
jan@42245
|
26 |
"of_nat n = Complex (of_nat n) 0"
|
jan@42245
|
27 |
by (induct n) (simp_all add: complex_one_def)
|
jan@42245
|
28 |
|
jan@42245
|
29 |
|
jan@42245
|
30 |
text {* The following two lemmas are useful for experimenting with the
|
jan@42245
|
31 |
transformations, at a vector length of four. *}
|
jan@42245
|
32 |
|
jan@42245
|
33 |
lemma Ivl4:
|
jan@42245
|
34 |
"{0..<4::nat} = {0, 1, 2, 3}"
|
jan@42245
|
35 |
proof -
|
jan@42245
|
36 |
have "{0..<4::nat} = {0..<Suc (Suc (Suc (Suc 0)))}" by (simp add: eval_nat_numeral)
|
jan@42245
|
37 |
also have "... = {0, 1, 2, 3}"
|
jan@42245
|
38 |
by (simp add: atLeastLessThanSuc eval_nat_numeral insert_commute)
|
jan@42245
|
39 |
finally show ?thesis .
|
jan@42245
|
40 |
qed
|
jan@42245
|
41 |
|
jan@42245
|
42 |
lemma Sum4:
|
jan@42245
|
43 |
"(\<Sum>i=0..<4::nat. x i) = x 0 + x 1 + x 2 + x 3"
|
jan@42245
|
44 |
by (simp add: Ivl4 eval_nat_numeral)
|
jan@42245
|
45 |
|
jan@42245
|
46 |
|
jan@42245
|
47 |
text {* A number of specialised lemmas for the summation operator,
|
jan@42245
|
48 |
where the index set is the natural numbers *}
|
jan@42245
|
49 |
|
jan@42245
|
50 |
lemma setsum_add_nat_ivl_singleton:
|
jan@42245
|
51 |
assumes less: "m < (n::nat)"
|
jan@42245
|
52 |
shows "f m + setsum f {m<..<n} = setsum f {m..<n}"
|
jan@42245
|
53 |
proof -
|
jan@42245
|
54 |
have "f m + setsum f {m<..<n} = setsum f ({m} \<union> {m<..<n})"
|
jan@42245
|
55 |
by (simp add: setsum_Un_disjoint ivl_disj_int)
|
jan@42245
|
56 |
also from less have "... = setsum f {m..<n}"
|
jan@42245
|
57 |
by (simp only: ivl_disj_un)
|
jan@42245
|
58 |
finally show ?thesis .
|
jan@42245
|
59 |
qed
|
jan@42245
|
60 |
|
jan@42245
|
61 |
lemma setsum_add_split_nat_ivl_singleton:
|
jan@42245
|
62 |
assumes less: "m < (n::nat)"
|
jan@42245
|
63 |
and g: "!!i. [| m < i; i < n |] ==> g i = f i"
|
jan@42245
|
64 |
shows "f m + setsum g {m<..<n} = setsum f {m..<n}"
|
jan@42245
|
65 |
using less g
|
jan@42245
|
66 |
by(simp add: setsum_add_nat_ivl_singleton cong: strong_setsum_cong)
|
jan@42245
|
67 |
|
jan@42245
|
68 |
lemma setsum_add_split_nat_ivl:
|
jan@42245
|
69 |
assumes le: "m <= (k::nat)" "k <= n"
|
jan@42245
|
70 |
and g: "!!i. [| m <= i; i < k |] ==> g i = f i"
|
jan@42245
|
71 |
and h: "!!i. [| k <= i; i < n |] ==> h i = f i"
|
jan@42245
|
72 |
shows "setsum g {m..<k} + setsum h {k..<n} = setsum f {m..<n}"
|
jan@42245
|
73 |
using le g h by (simp add: setsum_add_nat_ivl cong: strong_setsum_cong)
|
jan@42245
|
74 |
|
jan@42245
|
75 |
lemma ivl_splice_Un:
|
jan@42245
|
76 |
"{0..<2*n::nat} = (op * 2 ` {0..<n}) \<union> ((%i. Suc (2*i)) ` {0..<n})"
|
jan@42245
|
77 |
apply (unfold image_def Bex_def)
|
jan@42245
|
78 |
apply auto
|
jan@42245
|
79 |
apply arith
|
jan@42245
|
80 |
done
|
jan@42245
|
81 |
|
jan@42245
|
82 |
lemma ivl_splice_Int:
|
jan@42245
|
83 |
"(op * 2 ` {0..<n}) \<inter> ((%i. Suc (2*i)) ` {0..<n}) = {}"
|
jan@42245
|
84 |
by auto arith
|
jan@42245
|
85 |
|
jan@42245
|
86 |
lemma double_inj_on:
|
jan@42245
|
87 |
"inj_on (%i. 2*i::nat) A"
|
jan@42245
|
88 |
by (simp add: inj_onI)
|
jan@42245
|
89 |
|
jan@42245
|
90 |
lemma Suc_double_inj_on:
|
jan@42245
|
91 |
"inj_on (%i. Suc (2*i)) A"
|
jan@42245
|
92 |
by (rule inj_onI) simp
|
jan@42245
|
93 |
|
jan@42245
|
94 |
lemma setsum_splice:
|
jan@42245
|
95 |
"(\<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
|
96 |
proof -
|
jan@42245
|
97 |
have "(\<Sum>i::nat = 0..<2*n. f i) =
|
jan@42245
|
98 |
setsum f (op * 2 ` {0..<n}) + setsum f ((%i. 2*i+1) ` {0..<n})"
|
jan@42245
|
99 |
by (simp add: ivl_splice_Un ivl_splice_Int setsum_Un_disjoint)
|
jan@42245
|
100 |
also have "... = (\<Sum>i = 0..<n. f (2*i)) + (\<Sum>i = 0..<n. f (2*i+1))"
|
jan@42245
|
101 |
by (simp add: setsum_reindex [OF double_inj_on]
|
jan@42245
|
102 |
setsum_reindex [OF Suc_double_inj_on])
|
jan@42245
|
103 |
finally show ?thesis .
|
jan@42245
|
104 |
qed
|
jan@42245
|
105 |
|
jan@42245
|
106 |
|
jan@42245
|
107 |
section {* Complex Roots of Unity *}
|
jan@42245
|
108 |
|
jan@42245
|
109 |
text {* The function @{term cis} from the complex library returns the
|
jan@42245
|
110 |
point on the unity circle corresponding to the argument angle. It
|
jan@42245
|
111 |
is the base for our definition of @{text root}. The main property,
|
jan@42245
|
112 |
De Moirve's formula is already there in the library. *}
|
jan@42245
|
113 |
|
jan@42245
|
114 |
definition root :: "nat => complex" where
|
jan@42245
|
115 |
"root n == cis (2*pi/(real (n::nat)))"
|
jan@42245
|
116 |
|
jan@42245
|
117 |
lemma sin_periodic_pi_diff [simp]: "sin (x - pi) = - sin x"
|
jan@42245
|
118 |
by (simp add: sin_diff)
|
jan@42245
|
119 |
|
jan@42245
|
120 |
lemma sin_cos_between_zero_two_pi:
|
jan@42245
|
121 |
assumes 0: "0 < x" and pi: "x < 2 * pi"
|
jan@42245
|
122 |
shows "sin x \<noteq> 0 \<or> cos x \<noteq> 1"
|
jan@42245
|
123 |
proof -
|
jan@42245
|
124 |
{ assume "0 < x" and "x < pi"
|
jan@42245
|
125 |
then have "sin x \<noteq> 0" by (auto dest: sin_gt_zero_pi) }
|
jan@42245
|
126 |
moreover
|
jan@42245
|
127 |
{ assume "x = pi"
|
jan@42245
|
128 |
then have "cos x \<noteq> 1" by simp }
|
jan@42245
|
129 |
moreover
|
jan@42245
|
130 |
{ assume pi1: "pi < x" and pi2: "x < 2 * pi"
|
jan@42245
|
131 |
then have "0 < x - pi" and "x - pi < pi" by arith+
|
jan@42245
|
132 |
then have "sin (x - pi) \<noteq> 0" by (auto dest: sin_gt_zero_pi)
|
jan@42245
|
133 |
with pi1 pi2 have "sin x \<noteq> 0" by simp }
|
jan@42245
|
134 |
ultimately show ?thesis using 0 pi by arith
|
jan@42245
|
135 |
qed
|
jan@42245
|
136 |
|
jan@42245
|
137 |
|
jan@42245
|
138 |
subsection {* Basic Lemmas *}
|
jan@42245
|
139 |
|
jan@42245
|
140 |
lemma root_nonzero:
|
jan@42245
|
141 |
"root n ~= 0"
|
jan@42245
|
142 |
apply (unfold root_def)
|
jan@42245
|
143 |
apply (unfold cis_def)
|
jan@42245
|
144 |
apply auto
|
jan@42245
|
145 |
apply (drule sin_zero_abs_cos_one)
|
jan@42245
|
146 |
apply arith
|
jan@42245
|
147 |
done
|
jan@42245
|
148 |
|
jan@42245
|
149 |
lemma root_unity:
|
jan@42245
|
150 |
"root n ^ n = 1"
|
jan@42245
|
151 |
apply (unfold root_def)
|
jan@42245
|
152 |
apply (simp add: DeMoivre)
|
jan@42245
|
153 |
apply (simp add: cis_def)
|
jan@42245
|
154 |
done
|
jan@42245
|
155 |
|
jan@42245
|
156 |
lemma root_cancel:
|
jan@42245
|
157 |
"0 < d ==> root (d * n) ^ (d * k) = root n ^ k"
|
jan@42245
|
158 |
apply (unfold root_def)
|
jan@42245
|
159 |
apply (simp add: DeMoivre)
|
jan@42245
|
160 |
done
|
jan@42245
|
161 |
|
jan@42245
|
162 |
lemma root_summation:
|
jan@42245
|
163 |
assumes k: "0 < k" "k < n"
|
jan@42245
|
164 |
shows "(\<Sum>i=0..<n. (root n ^ k) ^ i) = 0"
|
jan@42245
|
165 |
proof -
|
jan@42245
|
166 |
from k have real0: "0 < real k * (2 * pi) / real n"
|
jan@42245
|
167 |
by (simp add: zero_less_divide_iff
|
jan@42245
|
168 |
mult_strict_right_mono [where a = 0, simplified])
|
jan@42245
|
169 |
from k mult_strict_right_mono [where a = "real k" and
|
jan@42245
|
170 |
b = "real n" and c = "2 * pi / real n", simplified]
|
jan@42245
|
171 |
have realk: "real k * (2 * pi) / real n < 2 * pi"
|
jan@42245
|
172 |
by (simp add: zero_less_divide_iff)
|
jan@42245
|
173 |
txt {* Main part of the proof *}
|
jan@42245
|
174 |
have "(\<Sum>i=0..<n. (root n ^ k) ^ i) =
|
jan@42245
|
175 |
((root n ^ k) ^ n - 1) / (root n ^ k - 1)"
|
jan@42245
|
176 |
apply (rule geometric_sum)
|
jan@42245
|
177 |
apply (unfold root_def)
|
jan@42245
|
178 |
apply (simp add: DeMoivre)
|
jan@42245
|
179 |
using real0 realk sin_cos_between_zero_two_pi
|
jan@42245
|
180 |
apply (auto simp add: cis_def complex_one_def)
|
jan@42245
|
181 |
done
|
jan@42245
|
182 |
also have "... = ((root n ^ n) ^ k - 1) / (root n ^ k - 1)"
|
jan@42245
|
183 |
by (simp add: power_mult [THEN sym] mult_ac)
|
jan@42245
|
184 |
also have "... = 0"
|
jan@42245
|
185 |
by (simp add: root_unity)
|
jan@42245
|
186 |
finally show ?thesis .
|
jan@42245
|
187 |
qed
|
jan@42245
|
188 |
|
jan@42245
|
189 |
lemma root_summation_inv:
|
jan@42245
|
190 |
assumes k: "0 < k" "k < n"
|
jan@42245
|
191 |
shows "(\<Sum>i=0..<n. ((1 / root n) ^ k) ^ i) = 0"
|
jan@42245
|
192 |
proof -
|
jan@42245
|
193 |
from k have real0: "0 < real k * (2 * pi) / real n"
|
jan@42245
|
194 |
by (simp add: zero_less_divide_iff
|
jan@42245
|
195 |
mult_strict_right_mono [where a = 0, simplified])
|
jan@42245
|
196 |
from k mult_strict_right_mono [where a = "real k" and
|
jan@42245
|
197 |
b = "real n" and c = "2 * pi / real n", simplified]
|
jan@42245
|
198 |
have realk: "real k * (2 * pi) / real n < 2 * pi"
|
jan@42245
|
199 |
by (simp add: zero_less_divide_iff)
|
jan@42245
|
200 |
txt {* Main part of the proof *}
|
jan@42245
|
201 |
have "(\<Sum>i=0..<n. ((1 / root n) ^ k) ^ i) =
|
jan@42245
|
202 |
(((1 / root n) ^ k) ^ n - 1) / ((1 / root n) ^ k - 1)"
|
jan@42245
|
203 |
apply (rule geometric_sum)
|
jan@42245
|
204 |
apply (simp add: nonzero_inverse_eq_divide [THEN sym] root_nonzero)
|
jan@42245
|
205 |
apply (unfold root_def)
|
jan@42245
|
206 |
apply (simp add: DeMoivre)
|
jan@42245
|
207 |
using real0 realk sin_cos_between_zero_two_pi
|
jan@42245
|
208 |
apply (auto simp add: cis_def complex_one_def)
|
jan@42245
|
209 |
done
|
jan@42245
|
210 |
also have "... = (((1 / root n) ^ n) ^ k - 1) / ((1 / root n) ^ k - 1)"
|
jan@42245
|
211 |
by (simp add: power_mult [THEN sym] mult_ac)
|
jan@42245
|
212 |
also have "... = 0"
|
jan@42245
|
213 |
by (simp add: power_divide root_unity)
|
jan@42245
|
214 |
finally show ?thesis .
|
jan@42245
|
215 |
qed
|
jan@42245
|
216 |
|
jan@42245
|
217 |
lemma root0 [simp]:
|
jan@42245
|
218 |
"root 0 = 1"
|
jan@42245
|
219 |
by (simp add: root_def cis_def)
|
jan@42245
|
220 |
|
jan@42245
|
221 |
lemma root1 [simp]:
|
jan@42245
|
222 |
"root 1 = 1"
|
jan@42245
|
223 |
by (simp add: root_def cis_def)
|
jan@42245
|
224 |
|
jan@42245
|
225 |
lemma root2 [simp]:
|
jan@42245
|
226 |
"root 2 = Complex -1 0"
|
jan@42245
|
227 |
by (simp add: root_def cis_def)
|
jan@42245
|
228 |
|
jan@42245
|
229 |
lemma root4 [simp]:
|
jan@42245
|
230 |
"root 4 = ii"
|
jan@42245
|
231 |
by (simp add: root_def cis_def)
|
jan@42245
|
232 |
|
jan@42245
|
233 |
|
jan@42245
|
234 |
subsection {* Derived Lemmas *}
|
jan@42245
|
235 |
|
jan@42245
|
236 |
lemma root_cancel1:
|
jan@42245
|
237 |
"root (2 * m) ^ (i * (2 * j)) = root m ^ (i * j)"
|
jan@42245
|
238 |
proof -
|
jan@42245
|
239 |
have "root (2 * m) ^ (i * (2 * j)) = root (2 * m) ^ (2 * (i * j))"
|
jan@42245
|
240 |
by (simp add: mult_ac)
|
jan@42245
|
241 |
also have "... = root m ^ (i * j)"
|
jan@42245
|
242 |
by (simp add: root_cancel)
|
jan@42245
|
243 |
finally show ?thesis .
|
jan@42245
|
244 |
qed
|
jan@42245
|
245 |
|
jan@42245
|
246 |
lemma root_cancel2:
|
jan@42245
|
247 |
"0 < n ==> root (2 * n) ^ n = - 1"
|
jan@42245
|
248 |
txt {* Note the space between @{text "-"} and @{text "1"}. *}
|
jan@42245
|
249 |
using root_cancel [where n = 2 and k = 1]
|
jan@42245
|
250 |
apply (simp only: mult_ac)
|
jan@42245
|
251 |
apply (simp add: complex_one_def)
|
jan@42245
|
252 |
done
|
jan@42245
|
253 |
|
jan@42245
|
254 |
|
jan@42245
|
255 |
section {* Discrete Fourier Transformation *}
|
jan@42245
|
256 |
|
jan@42245
|
257 |
text {*
|
jan@42245
|
258 |
We define operations @{text DFT} and @{text IDFT} for the discrete
|
jan@42245
|
259 |
Fourier Transform and its inverse. Vectors are simply functions of
|
jan@42245
|
260 |
type @{text "nat => complex"}. *}
|
jan@42245
|
261 |
|
jan@42245
|
262 |
text {*
|
jan@42245
|
263 |
@{text "DFT n a"} is the transform of vector @{text a}
|
jan@42245
|
264 |
of length @{text n}, @{text IDFT} its inverse. *}
|
jan@42245
|
265 |
|
jan@42245
|
266 |
definition DFT :: "nat => (nat => complex) => (nat => complex)" where
|
jan@42245
|
267 |
"DFT n a == (%i. \<Sum>j=0..<n. (root n) ^ (i * j) * (a j))"
|
jan@42245
|
268 |
|
jan@42245
|
269 |
definition IDFT :: "nat => (nat => complex) => (nat => complex)" where
|
jan@42245
|
270 |
"IDFT n a == (%i. (\<Sum>k=0..<n. (a k) / (root n) ^ (i * k)))"
|
jan@42245
|
271 |
|
jan@42245
|
272 |
schematic_lemma "map (DFT 4 a) [0, 1, 2, 3] = ?x"
|
jan@42245
|
273 |
by(simp add: DFT_def Sum4)
|
jan@42245
|
274 |
|
jan@42245
|
275 |
text {* Lemmas for the correctness proof. *}
|
jan@42245
|
276 |
|
jan@42245
|
277 |
lemma DFT_lower:
|
jan@42245
|
278 |
"DFT (2 * m) a i =
|
jan@42245
|
279 |
DFT m (%i. a (2 * i)) i +
|
jan@42245
|
280 |
(root (2 * m)) ^ i * DFT m (%i. a (2 * i + 1)) i"
|
jan@42245
|
281 |
proof (unfold DFT_def)
|
jan@42245
|
282 |
have "(\<Sum>j = 0..<2 * m. root (2 * m) ^ (i * j) * a j) =
|
jan@42245
|
283 |
(\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j)) * a (2 * j)) +
|
jan@42245
|
284 |
(\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j + 1)) * a (2 * j + 1))"
|
jan@42245
|
285 |
(is "?s = _")
|
jan@42245
|
286 |
by (simp add: setsum_splice)
|
jan@42245
|
287 |
also have "... = (\<Sum>j = 0..<m. root m ^ (i * j) * a (2 * j)) +
|
jan@42245
|
288 |
root (2 * m) ^ i *
|
jan@42245
|
289 |
(\<Sum>j = 0..<m. root m ^ (i * j) * a (2 * j + 1))"
|
jan@42245
|
290 |
(is "_ = ?t")
|
jan@42245
|
291 |
txt {* First pair of sums *}
|
jan@42245
|
292 |
apply (simp add: root_cancel1)
|
jan@42245
|
293 |
txt {* Second pair of sums *}
|
jan@42245
|
294 |
apply (simp add: setsum_right_distrib)
|
jan@42245
|
295 |
apply (simp add: power_add)
|
jan@42245
|
296 |
apply (simp add: root_cancel1)
|
jan@42245
|
297 |
apply (simp add: mult_ac)
|
jan@42245
|
298 |
done
|
jan@42245
|
299 |
finally show "?s = ?t" .
|
jan@42245
|
300 |
qed
|
jan@42245
|
301 |
|
jan@42245
|
302 |
lemma DFT_upper:
|
jan@42245
|
303 |
assumes mbound: "0 < m" and ibound: "m <= i"
|
jan@42245
|
304 |
shows "DFT (2 * m) a i =
|
jan@42245
|
305 |
DFT m (%i. a (2 * i)) (i - m) -
|
jan@42245
|
306 |
root (2 * m) ^ (i - m) * DFT m (%i. a (2 * i + 1)) (i - m)"
|
jan@42245
|
307 |
proof (unfold DFT_def)
|
jan@42245
|
308 |
have "(\<Sum>j = 0..<2 * m. root (2 * m) ^ (i * j) * a j) =
|
jan@42245
|
309 |
(\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j)) * a (2 * j)) +
|
jan@42245
|
310 |
(\<Sum>j = 0..<m. root (2 * m) ^ (i * (2 * j + 1)) * a (2 * j + 1))"
|
jan@42245
|
311 |
(is "?s = _")
|
jan@42245
|
312 |
by (simp add: setsum_splice)
|
jan@42245
|
313 |
also have "... =
|
jan@42245
|
314 |
(\<Sum>j = 0..<m. root m ^ ((i - m) * j) * a (2 * j)) -
|
jan@42245
|
315 |
root (2 * m) ^ (i - m) *
|
jan@42245
|
316 |
(\<Sum>j = 0..<m. root m ^ ((i - m) * j) * a (2 * j + 1))"
|
jan@42245
|
317 |
(is "_ = ?t")
|
jan@42245
|
318 |
txt {* First pair of sums *}
|
jan@42245
|
319 |
apply (simp add: root_cancel1)
|
jan@42245
|
320 |
apply (simp add: root_unity ibound root_nonzero power_diff power_mult)
|
jan@42245
|
321 |
txt {* Second pair of sums *}
|
jan@42245
|
322 |
apply (simp add: mbound root_cancel2)
|
jan@42245
|
323 |
apply (simp add: setsum_right_distrib)
|
jan@42245
|
324 |
apply (simp add: power_add)
|
jan@42245
|
325 |
apply (simp add: root_cancel1)
|
jan@42245
|
326 |
apply (simp add: power_mult)
|
jan@42245
|
327 |
apply (simp add: mult_ac)
|
jan@42245
|
328 |
done
|
jan@42245
|
329 |
finally show "?s = ?t" .
|
jan@42245
|
330 |
qed
|
jan@42245
|
331 |
|
jan@42245
|
332 |
lemma IDFT_lower:
|
jan@42245
|
333 |
"IDFT (2 * m) a i =
|
jan@42245
|
334 |
IDFT m (%i. a (2 * i)) i +
|
jan@42245
|
335 |
(1 / root (2 * m)) ^ i * IDFT m (%i. a (2 * i + 1)) i"
|
jan@42245
|
336 |
proof (unfold IDFT_def)
|
jan@42245
|
337 |
have "(\<Sum>j = 0..<2 * m. a j / root (2 * m) ^ (i * j)) =
|
jan@42245
|
338 |
(\<Sum>j = 0..<m. a (2 * j) / root (2 * m) ^ (i * (2 * j))) +
|
jan@42245
|
339 |
(\<Sum>j = 0..<m. a (2 * j + 1) / root (2 * m) ^ (i * (2 * j + 1)))"
|
jan@42245
|
340 |
(is "?s = _")
|
jan@42245
|
341 |
by (simp add: setsum_splice)
|
jan@42245
|
342 |
also have "... = (\<Sum>j = 0..<m. a (2 * j) / root m ^ (i * j)) +
|
jan@42245
|
343 |
(1 / root (2 * m)) ^ i *
|
jan@42245
|
344 |
(\<Sum>j = 0..<m. a (2 * j + 1) / root m ^ (i * j))"
|
jan@42245
|
345 |
(is "_ = ?t")
|
jan@42245
|
346 |
txt {* First pair of sums *}
|
jan@42245
|
347 |
apply (simp add: root_cancel1)
|
jan@42245
|
348 |
txt {* Second pair of sums *}
|
jan@42245
|
349 |
apply (simp add: setsum_right_distrib)
|
jan@42245
|
350 |
apply (simp add: power_add)
|
jan@42245
|
351 |
apply (simp add: nonzero_power_divide root_nonzero)
|
jan@42245
|
352 |
apply (simp add: root_cancel1)
|
jan@42245
|
353 |
done
|
jan@42245
|
354 |
finally show "?s = ?t" .
|
jan@42245
|
355 |
qed
|
jan@42245
|
356 |
|
jan@42245
|
357 |
lemma IDFT_upper:
|
jan@42245
|
358 |
assumes mbound: "0 < m" and ibound: "m <= i"
|
jan@42245
|
359 |
shows "IDFT (2 * m) a i =
|
jan@42245
|
360 |
IDFT m (%i. a (2 * i)) (i - m) -
|
jan@42245
|
361 |
(1 / root (2 * m)) ^ (i - m) *
|
jan@42245
|
362 |
IDFT m (%i. a (2 * i + 1)) (i - m)"
|
jan@42245
|
363 |
proof (unfold IDFT_def)
|
jan@42245
|
364 |
have "(\<Sum>j = 0..<2 * m. a j / root (2 * m) ^ (i * j)) =
|
jan@42245
|
365 |
(\<Sum>j = 0..<m. a (2 * j) / root (2 * m) ^ (i * (2 * j))) +
|
jan@42245
|
366 |
(\<Sum>j = 0..<m. a (2 * j + 1) / root (2 * m) ^ (i * (2 * j + 1)))"
|
jan@42245
|
367 |
(is "?s = _")
|
jan@42245
|
368 |
by (simp add: setsum_splice)
|
jan@42245
|
369 |
also have "... =
|
jan@42245
|
370 |
(\<Sum>j = 0..<m. a (2 * j) / root m ^ ((i - m) * j)) -
|
jan@42245
|
371 |
(1 / root (2 * m)) ^ (i - m) *
|
jan@42245
|
372 |
(\<Sum>j = 0..<m. a (2 * j + 1) / root m ^ ((i - m) * j))"
|
jan@42245
|
373 |
(is "_ = ?t")
|
jan@42245
|
374 |
txt {* First pair of sums *}
|
jan@42245
|
375 |
apply (simp add: root_cancel1)
|
jan@42245
|
376 |
apply (simp add: root_unity ibound root_nonzero power_diff power_mult)
|
jan@42245
|
377 |
txt {* Second pair of sums *}
|
jan@42245
|
378 |
apply (simp add: nonzero_power_divide root_nonzero)
|
jan@42245
|
379 |
apply (simp add: mbound root_cancel2)
|
jan@42245
|
380 |
apply (simp add: setsum_divide_distrib)
|
jan@42245
|
381 |
apply (simp add: power_add)
|
jan@42245
|
382 |
apply (simp add: root_cancel1)
|
jan@42245
|
383 |
apply (simp add: power_mult)
|
jan@42245
|
384 |
apply (simp add: mult_ac)
|
jan@42245
|
385 |
done
|
jan@42245
|
386 |
finally show "?s = ?t" .
|
jan@42245
|
387 |
qed
|
jan@42245
|
388 |
|
jan@42245
|
389 |
text {* @{text DFT} und @{text IDFT} are inverses. *}
|
jan@42245
|
390 |
|
jan@42245
|
391 |
declare divide_divide_eq_right [simp del]
|
jan@42245
|
392 |
divide_divide_eq_left [simp del]
|
jan@42245
|
393 |
|
jan@42245
|
394 |
lemma power_diff_inverse:
|
jan@42245
|
395 |
assumes nz: "(a::'a::field) ~= 0"
|
jan@42245
|
396 |
shows "m <= n ==> (inverse a) ^ (n-m) = (a^m) / (a^n)"
|
jan@42245
|
397 |
apply (induct n m rule: diff_induct)
|
jan@42245
|
398 |
apply (simp add: nonzero_power_inverse
|
jan@42245
|
399 |
nonzero_inverse_eq_divide [THEN sym] nz)
|
jan@42245
|
400 |
apply simp
|
jan@42245
|
401 |
apply (simp add: nz)
|
jan@42245
|
402 |
done
|
jan@42245
|
403 |
|
jan@42245
|
404 |
lemma power_diff_rev_if:
|
jan@42245
|
405 |
assumes nz: "(a::'a::field) ~= 0"
|
jan@42245
|
406 |
shows "(a^m) / (a^n) = (if n <= m then a ^ (m-n) else (1/a) ^ (n-m))"
|
jan@42245
|
407 |
proof (cases "n <= m")
|
jan@42245
|
408 |
case True with nz show ?thesis
|
jan@42245
|
409 |
by (simp add: power_diff)
|
jan@42245
|
410 |
next
|
jan@42245
|
411 |
case False with nz show ?thesis
|
jan@42245
|
412 |
by (simp add: power_diff_inverse nonzero_inverse_eq_divide [THEN sym])
|
jan@42245
|
413 |
qed
|
jan@42245
|
414 |
|
jan@42245
|
415 |
lemma power_divides_special:
|
jan@42245
|
416 |
"(a::'a::field) ~= 0 ==>
|
jan@42245
|
417 |
a ^ (i * j) / a ^ (k * i) = (a ^ j / a ^ k) ^ i"
|
jan@42245
|
418 |
by (simp add: nonzero_power_divide power_mult [THEN sym] mult_ac)
|
jan@42245
|
419 |
|
jan@42245
|
420 |
theorem DFT_inverse:
|
jan@42245
|
421 |
assumes i_less: "i < n"
|
jan@42245
|
422 |
shows "IDFT n (DFT n a) i = of_nat n * a i"
|
jan@42245
|
423 |
using [[linarith_split_limit = 0]]
|
jan@42245
|
424 |
apply (unfold DFT_def IDFT_def)
|
jan@42245
|
425 |
apply (simp add: setsum_divide_distrib)
|
jan@42245
|
426 |
apply (subst setsum_commute)
|
jan@42245
|
427 |
apply (simp only: times_divide_eq_left [THEN sym])
|
jan@42245
|
428 |
apply (simp only: power_divides_special [OF root_nonzero])
|
jan@42245
|
429 |
apply (simp add: power_diff_rev_if root_nonzero)
|
jan@42245
|
430 |
apply (simp add: setsum_divide_distrib [THEN sym]
|
jan@42245
|
431 |
setsum_left_distrib [THEN sym])
|
jan@42245
|
432 |
proof -
|
jan@42245
|
433 |
from i_less have i_diff: "!!k. i - k < n" by arith
|
jan@42245
|
434 |
have diff_i: "!!k. k < n ==> k - i < n" by arith
|
jan@42245
|
435 |
|
jan@42245
|
436 |
let ?sum = "%i j n. setsum (op ^ (if i <= j then root n ^ (j - i)
|
jan@42245
|
437 |
else (1 / root n) ^ (i - j))) {0..<n} * a j"
|
jan@42245
|
438 |
let ?sum1 = "%i j n. setsum (op ^ (root n ^ (j - i))) {0..<n} * a j"
|
jan@42245
|
439 |
let ?sum2 = "%i j n. setsum (op ^ ((1 / root n) ^ (i - j))) {0..<n} * a j"
|
jan@42245
|
440 |
|
jan@42245
|
441 |
from i_less have "(\<Sum>j = 0..<n. ?sum i j n) =
|
jan@42245
|
442 |
(\<Sum>j = 0..<i. ?sum2 i j n) + (\<Sum>j = i..<n. ?sum1 i j n)"
|
jan@42245
|
443 |
(is "?s = _")
|
jan@42245
|
444 |
by (simp add: root_summation_inv nat_dvd_not_less
|
jan@42245
|
445 |
setsum_add_split_nat_ivl [where f = "%j. ?sum i j n"])
|
jan@42245
|
446 |
also from i_less i_diff
|
jan@42245
|
447 |
have "... = (\<Sum>j = i..<n. ?sum1 i j n)"
|
jan@42245
|
448 |
by (simp add: root_summation_inv nat_dvd_not_less)
|
jan@42245
|
449 |
also from i_less have "... =
|
jan@42245
|
450 |
(\<Sum>j\<in>{i} \<union> {i<..<n}. ?sum1 i j n)"
|
jan@42245
|
451 |
by (simp only: ivl_disj_un)
|
jan@42245
|
452 |
also have "... =
|
jan@42245
|
453 |
(?sum1 i i n + (\<Sum>j\<in>{i<..<n}. ?sum1 i j n))"
|
jan@42245
|
454 |
by (simp add: setsum_Un_disjoint ivl_disj_int)
|
jan@42245
|
455 |
also from i_less diff_i have "... = ?sum1 i i n"
|
jan@42245
|
456 |
by (simp add: root_summation nat_dvd_not_less)
|
jan@42245
|
457 |
also from i_less have "... = of_nat n * a i" (is "_ = ?t")
|
jan@42245
|
458 |
by (simp add: of_nat_cplx)
|
jan@42245
|
459 |
finally show "?s = ?t" .
|
jan@42245
|
460 |
qed
|
jan@42245
|
461 |
|
jan@42245
|
462 |
|
jan@42245
|
463 |
section {* Discrete, Fast Fourier Transformation *}
|
jan@42245
|
464 |
|
jan@42245
|
465 |
text {* @{text "FFT k a"} is the transform of vector @{text a}
|
jan@42245
|
466 |
of length @{text "2 ^ k"}, @{text IFFT} its inverse. *}
|
jan@42245
|
467 |
|
jan@42245
|
468 |
primrec FFT :: "nat => (nat => complex) => (nat => complex)" where
|
jan@42245
|
469 |
"FFT 0 a = a"
|
jan@42245
|
470 |
| "FFT (Suc k) a =
|
jan@42245
|
471 |
(let (x, y) = (FFT k (%i. a (2*i)), FFT k (%i. a (2*i+1)))
|
jan@42245
|
472 |
in (%i. if i < 2^k
|
jan@42245
|
473 |
then x i + (root (2 ^ (Suc k))) ^ i * y i
|
jan@42245
|
474 |
else x (i- 2^k) - (root (2 ^ (Suc k))) ^ (i- 2^k) * y (i- 2^k)))"
|
jan@42245
|
475 |
|
jan@42245
|
476 |
primrec IFFT :: "nat => (nat => complex) => (nat => complex)" where
|
jan@42245
|
477 |
"IFFT 0 a = a"
|
jan@42245
|
478 |
| "IFFT (Suc k) a =
|
jan@42245
|
479 |
(let (x, y) = (IFFT k (%i. a (2*i)), IFFT k (%i. a (2*i+1)))
|
jan@42245
|
480 |
in (%i. if i < 2^k
|
jan@42245
|
481 |
then x i + (1 / root (2 ^ (Suc k))) ^ i * y i
|
jan@42245
|
482 |
else x (i - 2^k) -
|
jan@42245
|
483 |
(1 / root (2 ^ (Suc k))) ^ (i - 2^k) * y (i - 2^k)))"
|
jan@42245
|
484 |
|
jan@42245
|
485 |
text {* Finally, for vectors of length @{text "2 ^ k"},
|
jan@42245
|
486 |
@{text DFT} and @{text FFT}, and @{text IDFT} and
|
jan@42245
|
487 |
@{text IFFT} are equivalent. *}
|
jan@42245
|
488 |
|
jan@42245
|
489 |
theorem DFT_FFT:
|
jan@42245
|
490 |
"!!a i. i < 2 ^ k ==> DFT (2 ^ k) a i = FFT k a i"
|
jan@42245
|
491 |
proof (induct k)
|
jan@42245
|
492 |
case 0
|
jan@42245
|
493 |
then show ?case by (simp add: DFT_def)
|
jan@42245
|
494 |
next
|
jan@42245
|
495 |
case (Suc k)
|
jan@42245
|
496 |
assume i: "i < 2 ^ Suc k"
|
jan@42245
|
497 |
show ?case proof (cases "i < 2 ^ k")
|
jan@42245
|
498 |
case True
|
jan@42245
|
499 |
then show ?thesis apply simp apply (simp add: DFT_lower)
|
jan@42245
|
500 |
apply (simp add: Suc) done
|
jan@42245
|
501 |
next
|
jan@42245
|
502 |
case False
|
jan@42245
|
503 |
from i have "i - 2 ^ k < 2 ^ k" by simp
|
jan@42245
|
504 |
with False i show ?thesis apply simp apply (simp add: DFT_upper)
|
jan@42245
|
505 |
apply (simp add: Suc) done
|
jan@42245
|
506 |
qed
|
jan@42245
|
507 |
qed
|
jan@42245
|
508 |
|
jan@42245
|
509 |
theorem IDFT_IFFT:
|
jan@42245
|
510 |
"!!a i. i < 2 ^ k ==> IDFT (2 ^ k) a i = IFFT k a i"
|
jan@42245
|
511 |
proof (induct k)
|
jan@42245
|
512 |
case 0
|
jan@42245
|
513 |
then show ?case by (simp add: IDFT_def)
|
jan@42245
|
514 |
next
|
jan@42245
|
515 |
case (Suc k)
|
jan@42245
|
516 |
assume i: "i < 2 ^ Suc k"
|
jan@42245
|
517 |
show ?case proof (cases "i < 2 ^ k")
|
jan@42245
|
518 |
case True
|
jan@42245
|
519 |
then show ?thesis apply simp apply (simp add: IDFT_lower)
|
jan@42245
|
520 |
apply (simp add: Suc) done
|
jan@42245
|
521 |
next
|
jan@42245
|
522 |
case False
|
jan@42245
|
523 |
from i have "i - 2 ^ k < 2 ^ k" by simp
|
jan@42245
|
524 |
with False i show ?thesis apply simp apply (simp add: IDFT_upper)
|
jan@42245
|
525 |
apply (simp add: Suc) done
|
jan@42245
|
526 |
qed
|
jan@42245
|
527 |
qed
|
jan@42245
|
528 |
|
jan@42245
|
529 |
schematic_lemma "map (FFT (Suc (Suc 0)) a) [0, 1, 2, 3] = ?x"
|
jan@42245
|
530 |
by simp
|
jan@42245
|
531 |
|
jan@42245
|
532 |
end
|