Previous Up Next

6.26.3  Discrete Fourier Transform and the Fast Fourier Transform

For any integer N, the Discrete Fourier Transform (DFT) is a transformation FN defined on the set of periodic sequences of period N; it depends on a choice of a primitive N-th root of unity ωN. For sequences with complex coefficients, we take:

ωN=e
i π
N
 

If x is a periodic sequence of period N, defined by the vector x=[x0,x1,… xN−1] then FN(x)=y is a periodic sequence of period N, defined by:

(FN(x))k=yk=
N−1
j=0
xjωNk· j,

for k=0..N−1.

The Discrete Fourier Transform FN is bijective with inverse

 FN−1 =
1
N
 
FNN
     on  ℂ

i.e.:

(FNN−1(x))k=
1
N
N−1
j=0
xjωNk· j 

The Fast Fourier Transform (FFT) is an efficient way to compute the discrete Fourier transform; faster than computing each term individually. Xcas implements the FFT algorithm to compute the discrete Fourier transform when the period of the sequence is a power of 2.

The fft command computes the discrete Fourier transform.

The ifft command computes the inverse discrete Fourier transform.


Examples.

The properties of the Discrete Fourier Transform

Definitions. Let x and y be two periodic sequences of period N.


Properties.

FN(x · y)=



1
N



FN(x) * FN(y)
FN(x * y)=FN(x) · FN(y)

Applications

  1. Value of a polynomial
    Define a polynomial P(x)=∑j=0N−1cjxj by the vector of its coefficients c:=[c0,c1,..cN−1], where zeroes may be added so that N is a power of 2 (so the Fast Fourier Transform can be used).
  2. Trigonometric interpolation
    Let f be periodic function of period 2π and let fk=f(2kπ/N) for k=0..(N−1). Find a trigonometric polynomial p that interpolates f at xk=2kπ/N, that is find pj, j=0..N−1 such that
    p(x)= 
    N−1
    j=0
     pj xj,    p(xk)=fk
    Replacing xk by its value in p(x) we get:
    N−1
    j=0
     pj exp(i
    j2kπ
    N
    ) = fk
    In other words, (fk) is the inverse DFT of (pk), hence
    (pk)= 
    1
    N
     FN(  (fk)  ) 
    If the function f is real, pk=pk, hence depending whether N is even or odd:
    p(x)=p0+ 2 ℜ(
    N
    2
    −1
    k=0
    pkexp(ikx)) +ℜ(p
     
    N
    2
     exp(i
    Nx
    2
    )) 
    if N is even and
    p(x)=p0+ 2 ℜ(
    N−1
    2
    k=0
    pkexp(ikx))
    if N is odd.
  3. Fourier series
    Let f be a periodic function of period 2π and let yk= f(xk) where xk=2kπ/N for k=0..N−1. Suppose that the Fourier series of f converges to f (this will be the case if for example f is continuous). If N is large, a good approximation of f will be given by:
     
    N
    2
     ≤ n<
    N
    2
     cn exp(inx
    Hence we want a numeric approximation of
    cn=
    1
     


    0
    f(t)exp(−int)dt 
    The numeric value of the integral ∫0f(t)exp(−int)dt can be computed by the trapezoidal rule (note that the Romberg algorithm would not work here because the Euler Mac Laurin development has its coefficients equal to zero, since the integrated function is periodic, hence all its derivatives have the same value at 0 and at 2π). If c_n is the numeric value of cn obtained by the trapezoidal rule, then
    c_n=
    1
    N
    N−1
    k=0
    ykexp(−2i
    nkπ
    N
    ),     −
    N
    2
     ≤ n<
    N
    2
     
    Indeed, since xk=2kπ/N and f(xk)=yk:
                            f(xk)exp(−inxk)=
    ykexp(−2i
    nkπ
    N
    ), 
    f(0)exp(0)=f(2π)exp(−2i
    nNπ
    N
    )
    =y0=yN
    Hence:
    [c0,..c
     
    N
    2
    −1
    ,c
     
    N
    2
    +1
    ,..cN−1]=
    1
    N
    FN([y0,y1… y(N−1)]) 
    since


    Properties.


    Example.
    Input:

    f(t):=cos(t)+cos(2*t)
    x:=f(2*k*pi/8)$(k=0..7)

    Output:

    2,
    2
    2
    ,−1,−
    2
    2
    ,0,−
    2
    2
    ,−1,
    2
    2

    Input:

    fft(x)

    Output:


    0.0,4.0,4.0,0.0,0.0,0.0,4.0,4.0

    Dividing by N=8, you get

    c0=0,c1=0.5,c2=0.5,c3=0.0,
    c−4=0.0,c−3=0.0,c−2=0.5,=c−1=0.5

    Hence bk=0 and ak=ck+ck equals 1 for k=1,2 and 0 otherwise.

  4. Convolution Product
    If P(x)=∑j=0n−1ajxj and Q(x)=∑j=0m−1bjxj are given by the vectors of their coefficients a=[a0,a1,..an−1] and b=[b0,b1,..bm−1], you can compute the product of these two polynomials using the DFT. The product of polynomials is the convolution product of the periodic sequence of their coefficients if the period is greater or equal to (n+m). Therefore we complete a (resp. b) with m+p (resp. n+p) zeros, where p is chosen such that N=n+m+p is a power of 2. If a=[a0,a1,..an−1,0..0] and b=[b0,b1,..bm−1,0..0], then:
    P(x)Q(x)=
    n+m−1
    j=0
    (a*b)jxj 
    If you know FN(a) and FN(b), then ab=FN−1(FN(aFN(b)), since
    FN(x * y)=FN(x) · FN(y

Previous Up Next