5.23.2  Discrete Fourier Transform

Let N be an integer. The Discrete Fourier Transform (DFT) is a transformation F_N defined on the set of periodic sequences of period N, it depends on a choice of a primitive N-th root of unity ω_N. If the DFT is defined on sequences with complex coefficients, we take:

If x is a periodic sequence of period N, defined by the vector x=[x₀,x₁,...x_N−1] then F_N(x)=y is a periodic sequence of period N, defined by:

where ω_N is a primitive N-th root of unity. The discrete Fourier transform may be computed faster than by computing each y_k individually, by the Fast Fourier Transform (FFT). Xcas implements the FFT algorithm to compute the discrete Fourier transform only if N is a power of 2.

The properties of the Discrete Fourier Transform

The Discrete Fourier Transform F_N is a bijective transformation on periodic sequences such that

i.e. :

Inside Xcas the discrete Fourier transform and its inverse are denote by fft and ifft:

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

• The Hadamard product (notation ·) is defined by:

  (x · y)k = xk yk

• the convolution product (notation *) is defined by:

  (x * y)k=

  N−1 ∑ j=0

  xjyk−j

Properties :

Applications

1. Value of a polynomial
   Define a polynomial P(x)=∑_j=0^N−1c_jx^j by the vector of its coefficients c:=[c₀,c₁,..c_N−1], where zeroes may be added so that N is a power of 2.
   • Compute the values of P(x) at

     x=ak=ωN−k=exp(

     −2ikπ N

     ),    k=0..N−1

     This is just the discrete Fourier transform of c since

     P(ak)=

     N−1 ∑ j=0

     cj(ωN−k)j=FN(c)k

     Input, for example :
     P(x):=x+x^2; w:=i
     Here the coefficients of P are [0,1,1,0], N=4 and ω=exp(2iπ/4)=i.
     Input :
     fft([0,1,1,0])
     Output :
     [2,-1-i,0,-1+i]
     hence
     • P(1)=2,
     • P(-i)=P(w^-1)=-1-i,
     • P(-1)=P(w^-2)=0,
     • P(i)=P(w^-3)=-1+i.
   • Compute the values of P(x) at

     x=bk=ωNk=exp(

     2ikπ N

     ),    k=0..N−1

     This is N times the inverse fourier transform of c since

     P(ak)=

     N−1 ∑ j=0

     cj(ωNk)j=NFN−1(c)k

     Input, for example :
     P(x):=x+x^2 and w:=i
     Hence, the coefficients of P are [0,1,1,0], N=4 and ω=exp(2iπ/4)=i.
     Input :
     4*ifft([0,1,1,0])
     Output :
     [2,-1+i,0,-1-i]
     hence :
     • P(1)=2,
     • P(i)=P(w^1)=-1+i,
     • P(-1)=P(w^2)=0,
     • P(-i)=P(w^3)=-1-i.
     We find of course the same values as above...
2. Trigonometric interpolation
   Let f be periodic function of period 2π, assume that f(2kπ/N)=f_k for k=0..(N−1). Find a trigonometric polynomial p that interpolates f at x_k=2kπ/N, that is find p_j, j=0..N−1 such that

   p(x)=

   N−1 ∑ j=0

   pj exp(ijx),    p(xk)=fk

   Replacing x_k by its value in p(x) we get:

   N−1 ∑ j=0

   pj exp(i

   j2kπ N

   ) = fk

   In other words, (f_k) is the inverse DFT of (p_k), hence

   (pk)=

   1 N

   FN(  (fk)  )

   If the function f is real, p_−k=p_k, 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 ))  p(x) = p0+ 2 ℜ( N−1 2 ∑ k=0 pkexp(ikx))

3. Fourier series
   Let f be a periodic function of period 2π, such that

   f(xk)=yk,    xk=

   2kπ N

   , 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 2π

   ∫

   2π 0

   f(t)exp(−int)dt

   The numeric value of the integral ∫₀^2πf(t)exp(−int)dt may 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 c_n obtained by the trapezoidal rule, then

   c_n=

   1 2π

   2π N

   N−1 ∑ k=0

   ykexp(−2i

   nkπ N

   ),     −

   N 2

   ≤ n<

   N 2

   Indeed, since x_k=2kπ/N and f(x_k)=y_k:

   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
   • if n≥0, c_n=y_n
   • if n<0 c_n=y_n+N
   • ω_N=exp(2iπ/N), then ω_Nⁿ=ω_N^n+N

   Properties

   • The coefficients of the trigonometric polynomial that interpolates f at x=2kπ/N are

     pn=cn,    −

     N 2

     ≤ n<

     N 2

   • If f is a trigonometric polynomial P of degree m≤ N/2, then

     f(t)=P(t)=

     m−1 ∑ k=−m

     ckexp(2ikπ t)

     the trigonometric polynomial that interpolate f=P is P, the numeric approximation of the coefficients are in fact exact (c_n=c_n).
   • More generally, we can compute c_n−c_n.
     Suppose that f is equal to its Fourier series, i.e. that :
     

     f(t)=

     +∞ ∑ m=−∞

     cmexp(2iπ mt),

     +∞ ∑ m=−∞

     |cm|<∞

     Then :

     f(xk)=f(

     2kπ N

     )=yk=

     +∞ ∑ m=−∞

     cmωNkm,    c_n=

     1 N

     N−1 ∑ k=0

     ykωN−kn

     Replace y_k by its value in c_n:

     c_n=

     1 N

     N−1 ∑ k=0

     +∞ ∑ m=−∞

     cmωNkmωN−kn

     If m≠ n (mod N ), ω_N^m−n is an N-th root of unity different from 1, hence:

     ωN(m−n)N=1,

     N−1 ∑ k=0

     ωN(m−n)k=0

     Therefore, if m−n is a multiple of N (m=n+l· N) then ∑_k=0^N−1ω_N^k(m−n)=N, otherwise ∑_k=0^N−1ω_N^k(m−n)=0. By reversing the two sums, we get

     c_n = 1 N +∞ ∑ m=−∞ cm N−1 ∑ k=0 ωNk(m−n)    = +∞ ∑ l=−∞ c(n+l· N)    = ...cn−2· N+cn−N+cn+cn+N+cn+2· N+.....

     Conclusion: if |n|<N/2, c_n−c_n is a sum of c_j of large indexes (at least N/2 in absolute value), hence is small (depending on the rate of convergence of the Fourier series).

   Example Input :
   

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

   Then :

   x={2,sqrt(2)/2,-1,(-sqrt(2)/2,0,(-sqrt(2))/2,-1,sqrt(2)/2}
   fft(x)=[0.0,4.0,4.0,0.0,0.0,0.0,4.0,4.0]

   After a division by N=8, we get

   c₀=0,c₁=4.0/8,c₂=4.0/8,c₃=0.0,
   c₋₄=0.0,c₋₃=0.0,c₋₂=4.0/8,=c₋₁=4.0/8

   Hence b_k=0 and a_k=c_−k+c_k is equal to 1 if k=1,2 and 0 otherwise.

4. Convolution Product
   If P(x)=∑_j=0ⁿ⁻¹a_jx^j and Q(x)=∑_j=0^m−1b_jx^j are given by the vector of their coefficients a=[a₀,a₁,..a_n−1] and b=[b₀,b₁,..b_m−1], we may 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=[a₀,a₁,..a_n−1,0..0] and b=[b₀,b₁,..b_m−1,0..0], then:

   P(x)Q(x)=

   n+m−1 ∑ j=0

   (a*b)jxj

   We compute F_N(a), F_N(b), then ab=F_N⁻¹(F_N(a)· F_N(b)) using the properties

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