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:
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:
for k=0..N−1.
The Discrete Fourier Transform FN is bijective with inverse
i.e.:
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.
-
fft takes one argument:
x, a list or sequence regarded as one period of a periodic sequence.
- fft(x) returns FN(x), the discrete Fourier
transform of x.
If x has length which is a power of 2, then FN(x) is
computed with the Fast Fourier Transform.
The ifft command computes the inverse discrete Fourier transform.
-
ifft takes one argument:
x, a list or sequence regarded as one period of a periodic sequence.
- ifft(x) returns FN−1(x), the inverse
discrete Fourier transform of x.
If x has length which is a power of 2, then FN−1(x) is
computed with the Fast Fourier Transform.
Examples.
-
Input:
fft(0,1,1,0)
Output:
| ⎡
⎣ | 2.0,−1.0−i,0.0,−1.0+i | ⎤
⎦ |
- Input:
ifft([2,-1-i,0,-1+i])
Output:
The properties of the Discrete Fourier Transform
Definitions.
Let x and y be two periodic sequences of period N.
-
The Hadamard product (notation ·) is defined by:
- the convolution product (notation *) is defined by:
Properties.
FN(x · y) | = | ⎛
⎜
⎜
⎝ | | ⎞
⎟
⎟
⎠ | FN(x) * FN(y) |
|
FN(x * y) | = | FN(x) · FN(y)
|
|
Applications
-
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).
-
Compute the values of P(x) at
x=ak=ωN−k=exp( | | ), k=0..N−1
|
This is just the discrete Fourier transform of c since
Example.
Find the values of P(x+x2) at x=1,i,−1,−i.
Input:
P(x):=x+x^2
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.0,−1.0−i,0.0,−1.0+i | ⎤
⎦ |
Hence:
-
P(1)=2,
- P(−i)=P(ω−1)=−1−i,
- P(−1)=P(ω−2)=0,
- P(i)=P(ω−3)=−1+i.
- Compute the values of P(x) at
This is N times the inverse fourier transform of c since
Example.
Use this method to find the values of P(x+x2) at x=1,i,−1,−i.
Input:
P(x):=x+x^2
Again, 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.0,−1.0+i,0.0,−1.0−i | ⎤
⎦ |
Hence:
-
P(1)=2,
- P(i)=P(ω1)=−1+i,
- P(−1)=P(ω2)=0,
- P(−i)=P(ω3)=−1−i.
You find of course the same values as above.
- 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
Replacing xk by its value in p(x) we get:
In other words, (fk) is the inverse DFT of (pk), hence
If the function f is real, p−k=pk, hence depending
whether N is even or odd:
p(x)=p0+ 2 ℜ( | | pkexp(ikx))
+ℜ(p | | exp(i | | ))
|
if N is even and
if N is odd. - 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:
Hence we want a numeric approximation of
The numeric value of the integral ∫02πf(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
Indeed, since xk=2kπ/N and f(xk)=yk:
Hence:
[c0,..c | | ,c | | ,..cN−1]=
| | FN([y0,y1… y(N−1)])
|
since
-
if n≥0, cn=yn
- if n<0 cn=yn+N
- ωN=exp(2iπ/N), so ωNn=ωNn+N
Properties.
-
The coefficients of the trigonometric polynomial that interpolates f
at x=2kπ/N are
- If f is a trigonometric polynomial of degree m≤ N/2,
then
the trigonometric polynomial that interpolates f is f itself, the numeric
approximation of the coefficients are in fact exact (cn=cn).
- More generally, you can compute cn−cn.
Suppose that f is equal to its Fourier series, i.e. that:
f(t)= | | cmexp(2iπ mt),
| | |cm|<∞
|
Then:
f(xk)=f( | | )=yk= | | cmωNkm,
c_n= | | | ykωN−kn
|
Replace yk by its value in c_n:
If m≠ n (mod N ), ωNm−n is an N-th root of unity different
from 1, hence:
Therefore, if m−n is a multiple of N (m=n+l· N) then
∑k=0N−1ωNk(m−n)=N, otherwise
∑k=0N−1ωNk(m−n)=0.
By reversing the two sums, you get
c_n | = | |
| = | |
| = | … cn−2· N+cn−N+cn+cn+N+cn+2·
N+…
|
|
Conclusion: if |n|<N/2, then c_n−cn is a sum of cj
with large indices
(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)
|
Output:
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=c−k+ck equals 1 for k=1,2 and 0 otherwise.
- 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:
If you know FN(a) and FN(b), then a∗ b=FN−1(FN(a)·
FN(b)), since