5.23.2 Discrete Fourier Transform
Let N be an integer.
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. 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=[x0,x1,...xN−1] then
FN(x)=y is a periodic sequence of period N, defined by:
|(FN,ωN(x))k=yk=||xjωN−k· j, k=0..N−1 |
where ωN is a primitive N-th root of unity.
The discrete Fourier transform may be computed faster than by
computing each yk 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 FN is a bijective transformation
on periodic sequences such that
Inside Xcas the discrete Fourier transform and its inverse
are denote by fft and ifft:
fft(x)= FN(x), ifft(x)= FN−1(x)
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:
|N*FN(x · y)||=||FN(x) * FN(y)|
|FN(x * y)||=||FN(x) · FN(y)
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.
Compute the values of P(x) at
This is just the discrete Fourier transform of c since
Input, for example :
|x=ak=ωN−k=exp(||), k=0..N−1 |
Here the coefficients of P are [0,1,1,0],
N=4 and ω=exp(2iπ/4)=i.
- Compute the values of P(x) at
This is N times the inverse fourier transform of c since
Input, for example :
^2 and w:=i
Hence, the coefficients of P are [0,1,1,0],
N=4 and ω=exp(2iπ/4)=i.
We find of course the same values as above...
- Trigonometric interpolation
Let f be periodic function of period 2π, assume that f(2kπ/N)=fk
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)= || pj exp(ijx), p(xk)=fk|
2 ℜ(||pkexp(ikx))+ℜ(p|| exp(i||)) |
- Fourier series
Let f be a periodic function of period 2π, such that
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 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 cn obtained by the
trapezoidal rule, then
Indeed, since xk=2kπ/N and f(xk)=yk:
if n≥0, cn=yn
- if n<0 cn=yn+N
The coefficients of the trigonometric polynomial that interpolates f
at x=2kπ/N are
- If f is a trigonometric polynomial P of degree m≤ N/2,
the trigonometric polynomial that interpolate f=P is P, the numeric
approximation of the coefficients are in fact exact (cn=cn).
- More generally, we can compute cn−cn.
Suppose that f is equal to its Fourier series, i.e. that :
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
By reversing the two sums, we get
Conclusion: if |n|<N/2, c_n−cn is a sum of cj of large indexes
(at least N/2 in absolute value), hence is small (depending on the
rate of convergence of the Fourier series).
| ||=||...cn−2· N+cn−N+cn+cn+N+cn+2·
After a division by N=8, we get
Hence bk=0 and ak=c−k+ck is equal to 1 if k=1,2 and 0 otherwise.
- Convolution Product
are given by the vector of their coefficients
a=[a0,a1,..an−1] and b=[b0,b1,..bm−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=[a0,a1,..an−1,0..0] and b=[b0,b1,..bm−1,0..0], then:
We compute FN(a), FN(b), then ab=FN−1(FN(a)· FN(b))
using the properties
|NFN(x · y)=FN(x) * FN(y),
FN(x * y)=FN(x) · FN(y) |