VIEWS: 13 PAGES: 23 CATEGORY: Education POSTED ON: 6/24/2010
MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk University of Leicester November 5, 2009 MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 1 / 23 Introduction We receive some sort of signal MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 2 / 23 Fourier approximation In signal processing the need to rapidly compute approximate Fourier series: SN f (x) = αk (f )exp(ikx), |k|≤N with 2π 1 αk (f ) = f (x) exp(−ikx)dx, 2π 0 is fundamental. Approximate using all Newton-Cotes rules 2π 2N−1 2π 2πj 2πijk f (x) exp(−ikx)dx ≈ f exp − , |k| ≤ N. 0 2N 2N 2N j=0 MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 3 / 23 Discrete Fourier transform Thus it is important to compute the ˆ Deﬁnition 9.1 Discrete Fourier transform (DFT) v of a vector v of length N: N−1 2πijk ˆ vk = vj exp − , k = 0, 1, · · · , N − 1. N j=0 Remark 9.2 We will view the input v as being periodic (vk = vN+k ). History: Famous fast Fourier transform (FFT) algorithm by Cooley and Tukey 1965. First used by Gauss (1805). MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 4 / 23 In this lecture we will look at: 1 Properties of the DFT; 2 FFT algorithm. 3 Smoothing and compression. 4 What is wrong with Fourier series. 5 The wavelet paradigm. MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 5 / 23 Fourier basis The Fourier basis, 2πijk fj (k) = exp , k = 0, · · · , N − 1, N j = 0, 1, · · · , N − 1. Lemma 9.3 The sum of the roots of unity is 0. Proposition 9.4 The Fourier basis is orthogonal. Proof N−1 2πijk 2πilk f j , fl = exp exp − N N k=0 N−1 2πi(j − l)k = exp N k=0 N/d−1 2πik(j − l)/d = d exp = 0, N/d k=0 MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 6 / 23 Expanding in the Fourier basis We make an orthonormal basis noting that fj , fj = N, j = 0, 1, · · · , N − 1. Remark 9.5 f √j , j = 0, 1, · · · , N − 1, N is an orthonormal basis for RN . We can write any v ∈ RN in the form N−1 1 v= v, fj fj , N j=0 where N−1 2πijk v, fj = vj exp − ˆ = vj . N k=0 MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 7 / 23 The Inverse Fourier transform Remark 9.6 The Fourier transform computes the coeﬃcients of a vector with respect to the Fourier basis. N−1 1 v= ˆ vj fj , N j=0 means N−1 1 2πijk vk = vj exp ˆ , k = 0, 1, · · · , N − 1. N N j=0 So we recover the original coeﬃcients via the inverse transform above. MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 8 / 23 Parseval and Plancherel Inner product is independent of orthonormal basis: Theorem 9.7 (Plancherel’s theorem) 1 v, u = ˆ ˆ v, u , v, u ∈ RN . N Corollary 9.8 (Parseval’s theorem) 1 v =√ v , ˆ v ∈ RN . N MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 9 / 23 Matrix notation We can write the DFT in the form 1 1 1 ··· 1 1 w w2 ··· · · · w N−1 v= 1 ˆ w2 w4 ··· w 2(N−1) v, . . . .. . . . . . . . . . . 1 w N−1 w 2(N−1) w 3(N−1) · · · w (N−1)(N−1) where w = exp − 2πi . N Write v = WN v, where WN (i, j) = w (i−1)(j−1) . ˆ MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 10 / 23 For example For N = 4, w = exp −iπ/4 = −i. Then 1 1 1 1 1 w w2 w3 W4 = 1 w2 w4 w6 1 w3 w6 w9 1 1 1 1 1 −i −1 i = 1 −1 . 1 −1 1 i −1 −i MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 11 / 23 Filtering and convolution Fourier transform gives frequency information. In audio engineering we wish to manipulate the frequencies. Can we arrange: vj → uj = wj vj , ˆ ˆ ˆˆ j = 0, 1, · · · , N − 1. Solved by Deﬁnition 9. The convolution of v, w ∈ RN is N−1 [v ∗ w]j = vk wj−k , j = 0, 1, · · · , N − 1. k=0 Theorem 9. [v ∗ w]j = vj wj , ˆˆ j = 0, 1, · · · , N − 1. MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 12 / 23 Noise reduction Noisy signal has high frequency. Filter the high frequency. DFT v; 1 ˆ ˆ Make w by setting vj = 0, j > M; 2 3 ˆ IDFT w . Example 9. MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 13 / 23 Compression Smooth signal has few low frequencies. Filter the small frequencies. 1 DFT v; 2 Make w by setting vj = 0, if |ˆj | < ; ˆ ˆ v 3 ˆ IDFT w . MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 14 / 23 FFT idea First reorder rows and columns: 1 1 1 1 1 1 1 1 1 −i −1 i → 1 1 −1 −1 . 1 −1 1 −1 1 −1 −i i 1 i −1 −i 1 −1 i −i We need to keep track of changes in inputs and outputs - Binary Bit Reverse order 1 1 1 1 1 1 0 0 1 1 0 0 1 1 −1 −1 1 −1 0 0 0 0 1 1 = , 1 −1 −i i 0 0 1 1 1 −1 0 0 1 −1 i −i 0 0 1 −1 0 0 −i i which is a sparse factorisation. MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 15 / 23 Decimation formula Theorem 9.10 Let N = 2M. Then M−1 2πijk v2j ˆ = (vk + vk+M ) exp − M k=0 M−1 πik 2πijk v2j+1 ˆ = exp (vk − vk+M ) exp − , M M k=0 j = 0, 1, · · · , M − 1. Proof N−1 2πi(2j + 1)k v2j+1 ˆ = vk exp − N k=0 M−1 2πi(2j + 1)k 2πi(2j + 1)k + M = vk exp − + vk+M exp − 2M 2M k=0 M−1 πik 2πijk = exp − (vk − vk+M ) exp − . M M k=0 MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 16 / 23 Complexity Remark 9.11 We have replaced a length N DFT with 2 N/2 DFTs. Number of complex multiplications to compute DFT is O(N 2 ). After decimation we have 2 times (N/2)2 = N 2 /2, so a saving of 1/2. Now we iterate. The cost has been N/2 complex multiplications. MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 17 / 23 The FFT Suppose N = 2m . We can do m = log2 N decimations. Each takes N/2 complex operations. 2 times N/2 to N/4 is 2 times N/4 complex multiplications = N/2 multiplications. Total complexity N log2 N, 2 which brings large saving. MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 18 / 23 Translation operator Deﬁnition 9.9 The translation operator Tk : RN → RN is deﬁned by [Tl v]j = v (j − l), j = 0, · · · , N − 1. Periodic extention of v. MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 19 / 23 Lemma 9.10 2πijl [Tˆv]j = exp − l ˆ vj . N Proof N−1 2πijk [Tˆv]j l = [Tl v]k exp − N k=0 N−1 2πijk = v (k − l) exp − N k=0 N−1 2πij(k + l) = v (k) exp − N k=0 2πijl = exp − ˆ vj . N MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 20 / 23 Why wavelets? Theorem 9. Suppose {Tk v : k = 0, 1, · · · , N − 1} is an orthogonal basis for RN . Then |ˆk | = constant, v ∀k. Proof We have, using Plancherel Tk v, Tj v = Tk v, Tj v = 0, k = j. Using Lemma 9. we see that N−1 2πikl 2πijl Tk v, Tj v = exp − ˆ vl exp ˆ vl N N l=0 N−1 2πi(k − j)l = exp − |vl |2 = 0. N l=0 |ˆl |2 is orthogonal to fj (Fourier basis) for j = 1, 2, · · · , N − 1. v v parallel to f0 = [1, 1, · · · , 1]T . ˆ MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 21 / 23 Uncertainty principle Localisation in space means no localisation in frequency. Next step, use translates of two vectors: RN = span({T2k u} ∪ {T2k v} : k = 0, 1, · · · , N/2 − 1). This is the Wavelet Paradigm. MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 22 / 23 Summary In this lecture we have: 1 Properties of the DFT; 2 FFT algorithm. 3 Smoothing and compression. 4 What is wrong with Fourier series. 5 The wavelet paradigm. MAGIC015 Lecture 9: The Discrete Fourier Transform Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester) November 5, 2009 23 / 23