Your Federal Quarterly Tax Payments are due April 15th Get Help Now >>

MAGIC015 Lecture 9 The Discrete Fourier Transform by qxc16070

VIEWS: 13 PAGES: 23

									         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
                                                   ˆ
    Definition 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 coefficients 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 coefficients 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
    Definition 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




    Definition 9.9
    The translation operator Tk : RN → RN is defined 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

								
To top