#21 (11/08/2023)

Discrete Fourier Transforms (DFT)

It is often necessary to obtain Fourier series applied on a set of finite data points. If one uses the conventional Fourier series expansion1, (1) the infinite series must be truncated and (2) numerical integration must be used to compute the Fourier coefficient, cm, thus, introducing two layers of approximation.
The Discrete Fourier Transform (DFT) has advantages over the conventional Fourier series expansion in (1) requiring only N data points, (2) the original function can be reproduced exactly at the selected N points, (3) the coefficients of DFT do not require integration to be computed, (4) the DFT coefficients can be computed by Fast Fourier Transforms (FFT), and (5) When N is sufficiently large, DFT approaches the classical Fourier series.
cft1.jpg dft1.jpg
The table below is the difference and similarity between CFT and DFT.



CFT
f(x) = ∞
∑
k=−∞ 
 Fk ei k x
Fk = 1
2 π
 ⌠
⌡ 		2 π

0 
 f(x) e−i k x d x
DFT
fl = N−1
∑
k=0 
 Fk ei k [(2π)/(N)] l
Fk = 1
N
 N−1
∑
l = 0 
 fl e−i k [(2π)/(N)]l
(1)

Derivation

The derivation of the DFT coefficient, Fk, closely follows that of the CFT coefficient, ck. Note that x in CFT is roughly equivalent to [(2π)/(N)]l in DFT. The formulas for IDFT (inverse Discrete Fourier Transform) and DFT (Discrete Fourier Transform) are expressed as
fl
=
N−1
∑
k=0 
 Fk  ei k [(2π)/(N)] l,
(2)
Fk
=
1
N
 N−1
∑
l = 0 
 fl  e−i k [(2π)/(N)]l.
(3)
Multiply Eq. (2) by e−i m [(2 π)/(N)] l to obtain
fl  e−i m [(2 π)/(N)] l = N−1
∑
k=0 
 Fk  e i k [(2 π)/(N)]l e − i m [(2 π)/(N)]l.
(4)
Take summation of Eq. (4) over l to get
N−1
∑
l = 0 
 fl  e− im [(2π)/(N)] l
=
N−1
∑
l = 0 
 N−1
∑
k=0 
 Fk  e i k [(2 π)/(N)]l  e − i m [(2 π)/(N)]l
=
(interchange the order of summation)
=
N−1
∑
k=0 
 Fk N−1
∑
l = 0 
 ei (k−m) [(2 π)/(N)] l
=
⎛
⎜
⎜
⎝ 		Use N−1
∑
l = 0 
 ei (k−m) [(2π)/(N)]l = ⎧
⎪
⎨
⎪
⎩
N
k−m = 0, ±N, ±2N …
0
k−m ≠ 0, ±N, ±2 N …
 ⎞
⎟
⎟
⎠
=
N Fm.
(5)
Hence,
Fm = 1
N
 N − 1
∑
l = 0 
 fl  e− i  m [(2 π)/(N)] l.
(6)
Note that in Eq. (5), we used an identity

N−1
∑
l = 0 
 ei k [(2π)/(N)] l  e−i m [(2π)/(N)] l = ⎧
⎪
⎨
⎪
⎩
0
k ≠ m
N
k = m, m±N, m±2 N …
(7)
which is a discrete version of
⌠
⌡ 		π

−π 
 ei k x  e−i m x d x = ⎧
⎪
⎨
⎪
⎩
0
k ≠ m
2 π
k = m
(8)
(Proof of Identity (7))
Define
ω ≡ ei (k−m) [(2 π)/(N)].
so that
N−1
∑
l = 0 
 ei k [(2π)/(N)] l e−i m [(2π)/(N)] l
=
N−1
∑
l = 0 
 ωl
=
1 + ω+ ω2 + ω3 + …ωN−1.
(9)
If
k − m = 0, ±N, ±2 N, ±3N …,
ω = e0, ei [(2π)/(N)] N , ei [(2π)/(N)]2 N ... = 1 so that Eq. (9) becomes
N−1
∑
l = 0 
 ei k [(2π)/(N)] l  e−i m [(2π)/(N)] l = 1 + 1 + …+1 = N.
(10)
On the other hand, if
k − m ≠ 0, ±N, ±2 N, ±3N …,
Eq. (9) becomes
1 + ω+ ω2 + ω3 + …ωN−1
=
1 − ωN
1−ω
=
1− ei (k−m) 2 π
1−ei (k−m) [(2π)/(N)]
=
0     if     k −m ≠ 0, ±N, ±2N …
(11)
Therefore, it was shown that
N−1
∑
l = 0 
 ei (k−m) [(2π)/(N)]l = ⎧
⎪
⎨
⎪
⎩
N
k−m = 0, ±N, ±2N …
0
k−m ≠ 0, ±N, ±2 N …
(12)
Note that only half of Fk's are independent, i.e.

Fk =
FN−k
 
 .
(13)
The power spectrum, Pk, is symmetrical around its center, i.e.

Pk = PN−k.
(14)

Relation between DFT and CFT

The classical Fourier series expansion of f(x) over [0, 2 π] is
f(x)
=
∞
∑
k=−∞ 
 ck ei k x ,
(15)
ck
=
1
2π
 ⌠
⌡ 		2π

0 
 f(x) e− i k x dx.
(16)
If we consider an interval, 0 < x < 2 π, to be sampled at 0, [(2 π)/(N)], [(4 π)/(N)], [(6 π)/(N)] …, denoted as
fl ≡ f⎛
⎝ 		2 π
N
 l⎞
⎠ 		,     ω ≡ e[(2 π)/(N)]i .
(17)
Then, Eq. (3) becomes
Fk
=
1
N
 N−1
∑
l = 0 
 fl  ω− k l
=
1
N
 N−1
∑
l = 0 
 ⎛
⎝ 		∞
∑
k′=−∞ 
 ck′  ωk′l⎞
⎠ 		ω−k l
=
∞
∑
k′=−∞ 
 ck′⎛
⎝ 		1
N
 N−1
∑
l = 0 
 ω(k′−k)l⎞
⎠ 		.
=
⎛
⎜
⎜
⎝ 		Use N−1
∑
l = 0 
 ei (k−k′) [(2π)/(N)]l = ⎧
⎪
⎨
⎪
⎩
N
k − k′ = 0, ±N, ±2N …
0
k − k′ ≠ 0, ±N, ±2 N …
 ⎞
⎟
⎟
⎠
=
∞
∑
k′=−∞ 
 ck′ ×(1   if and only if k′=k, k±N, k±2N, k±3N, …, otherwise 0)
=
ck + ck±N + ck±2 N + ck±3N +…
=
ck + ∞
∑
l=1 
 ck±N l.
(18)
dft_cft.jpg
Equation (18) is the discrete version of the Nyquist sampling theorem. If the data set does not contain components whose frequencies are higher than N/2 or |k| < N/2, DFT is identical to CFS (classical Fourier series).

Convolution Theorem of DFT

The convolution between fl and gl is defined as
fl*gl ≡ 1
N
 N−1
∑
m=0 
 fl−m  gm.
(19)
Note that
fl*gl = gl * fl.
(DFT Convolution Theorem)

fl* gl = 1
N
 N−1
∑
k=0 
 Fk  Gk  ei [(2π)/(N)] k,
or DFT of fl* gl is Fk Gk.

Parseval's theorem of DFT


1
N
 N−1
∑
l = 0 
 fl  
-
g
 
l 
 = N−1
∑
k=0 
 Fk 
-
G
 
k 
 ,
(20)

1
N
 N−1
∑
l = 0 
 ||fl||2 = N−1
∑
k=0 
 || Fk||2.
(21)
Power spectrum in DFT is defined as
Pk ≡ ||Fk||2.
(22)

FFT

The DFT coefficients are given by Eq. (3) and are listed here again for reference

Fk = 1
N
 N−1
∑
l = 0 
 fl  ω−kl    k = 0, 1, 2, …, (N−1),     ω ≡ e2 πi/N.
(23)
For each Fk, provided that fl's and ω−k's are already stored in computer memory, the number of multiplications required is N and the number of required additions is N−1 so the total number of multiplications is N2 and the total number of additions is N (N−1). If N=1000, the total number of multiplications is one million and the total number of additions is 999,000, thus, requiring close to 2 million operations.
FFT (Fast Fourier Transform) is an algorithm to efficiently compute the DFT coefficients which can reduce the number of operations from N2 (multiplications) to N log2 N. This reduction is dramatic as shown in the table below:
N N2 N log2 N Ratio
10 100 33 3
100 10000 664 15
1000 1,000,000 9,900 100
10000 100 million 0.13 million 750
1 million 1 trillion 20 million 50,100
Given f0, f1, f2, …fN−1, define

g(x) ≡ f0 + f1 x + f2 x2 + …+ fN−1 xN−1.
(24)
By comparing Eq. (24) with Eq. (23), it can be readily seen that the coefficient, Fk, can be computed by

Fk = 1
N
 g(ω−k).
(25)
To illustrate the algorithm of FFT, let N=8, then,

g(x)
=
f0 + f1 x + f2 x2 + f3 x3 + f4 x4 + f5 x5 + f6 x6 + f7 x7
=
(f0 + f2 x2 + f4 x4 + f6 x6) +x ( f1 + f3 x2 + f5 x4+ f7 x6)
=
p(x2) + x q(x2),
(26)
where

⎧
⎪
⎨
⎪
⎩
p(x) ≡ f0 + f2 x + f4 x2 + f6 x3
q(x) ≡ f1 + f3 x + f5 x2 + f7 x3.
(27)
Thus, to compute all the coefficients of Fk is equivalent to computing both p(x2) and q(x2) at x=1, ω, ω2, ω3, ω4, ω5, ω6, ω7, multiplying x by q(x2) and add them together. It should be noted that for

x=1, ω, ω2, ω3, ω4, ω5, ω6, ω7,
(28)
the corresponding x2's are (note ω8=1)

x2
=
1, ω2, ω4, ω6, ω8, ω10, ω12, ω14
=
1, ω2, ω4, ω6, 1, ω2, ω4, ω6.
(29)
It is seen that the second half of x2 are the repetition of the first half, thus, it is NOT necessary to compute p(x2) and q(x2) when x belongs to the second half. Moreover, note that the set of {1, ω2, ω4, ω6} is the roots to

x4=1,
(30)
so computing p(x2) and q(x2) is equivalent to computing p(x) and q(x) where x is the root for x4=1 (a half of the order 8).
Therefore, if the number of multiplications to compute FFT of the N-th order is denoted as M(N), it must satisfy

M(N) = 2 M( N
2
 ) + N.
(31)
where the factor "2" comes from the computation of p(x2) and q(x2), the argument, N/2, of M comes from the equation of x4 = 1 (x(N/2)=1) and the last term of N is for the multiplication between x and q(x2) for N times.
When N=1, g(x)=f0, thus,
M(0)=1
(32)
From Eqs.(31-32), the solution to Eq. (31) is

M(N ) = N log2 N.2
(33)
Similarly, the number of additions can be computed as

A (N) = N log2 N.
(34)

Example 1

Consider a data set of (N=16)

f1 = {0, 1, 2, 3, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
(35)
which looks like
fft_2.jpg
The power spectrum, |Ck|2, of f1 is
fft_3.jpg
Now consider

f2={0, 0, 0, 0, 0, 0, 1, 2, 3, 2, 1, 0, 0, 0, 0, 0}
(36)
fft_4.jpg
The power spectrum, |Ck|2, of f2 is
fft_5.jpg

Example 2 (2-D)

Consider two different images
fft2d_1.jpg,fft2d_2.jpg
The power spectra of the images are shown as
fft2d_3.jpg,fft2d_4.jpg
They are identical each other even though the spatial representations are different. It can be indeed shown that the correlation between the two graphs of the power spectra is 100%.

Example 3

Using the fft (Fast Fourier Transform) and ifft (Inverse Fourier Transform) functions available in MATLAB/ OCTAVE, remove the noise from the data above and restore the original signal.
You can copy and paste the following line to enter the data into MATLAB/OCTAVE. 
y=[2.75428, 5.71513, 4.74625, 4.0357, 3.58133, 4.29361, 3.17913,3.96698, 4.947, 2.50101, 5.18169, 2.42088, 4.19325, 5.36498, 3.02793,5.3784, 2.79834, 5.13506, 5.43307, 5.70375, 4.02817, 3.45414,3.15312, 3.77185, 3.1477, 3.61882, 4.29079, 1.62228, 5.45282,5.21009, 2.99312, 3.5315, 2.37493, 4.56933, 3.66105, 2.94791,4.00497, 5.01201, 2.42359, 3.34126, 2.95821, 1.60697, 2.69763,3.32045, 4.58761, 3.78392, 1.14806, 1.12364, 2.98555, 1.68054,2.34178, 0.954301, 2.95354, 1.85866, 0.70389, 2.74476, 1.86708,2.5442, 2.26406, 0.984498, 3.01632, 2.65318, 0.928562, 2.69881,1.08165, 2.03826, 3.19213, 4.30941, 1.39572, 3.1275, 2.88962,4.00469, 3.20351, 2.21581, 1.29336, 3.77378, 0.952478, 1.04008,1.25416, 1.67683, 3.71773, 3.11416, 3.59578, 1.28689, 1.28631,1.03774, 3.23732, 3.15302, 4.76595, 3.5586, 4.6036, 3.40278, 3.93159,0.996068, 2.28392, 3.97436, 1.31199, 3.37306, 1.59347, 4.81479,4.98611, 2.97305, 4.99383, 3.80786, 1.95455, 4.56209, 2.11904,3.26038, 1.42664, 4.30224, 3.67952, 4.92539, 1.49919, 1.60278,3.95577, 2.42314, 2.48873, 1.54822, 2.63348, 3.43014, 2.17739,3.19454, 2.07749, 3.67022, 2.26272, 1.30863, 2.18551, 4.97789,1.43625, 1.88595, 1.21594, 2.87569, 1.17509, 2.75502, 2.71201,3.12858, 4.85521, 2.33088, 3.93246, 1.87746, 3.53247, 1.94075,2.44728, 3.5839, 2.47802, 3.85354, 1.45955, 0.983903, 1.264, 3.57013,0.273807, 0.978381, 0.770786, 2.59596, -0.0634878, 2.94635, 0.40255,0.609216, 1.95247, 0.503754, 0.189706, 2.87551, 2.57051, 3.12867,1.11012, 1.33686, 0.473327, -0.158264, 2.87489, 1.67097, -0.85459,-1.06732, 1.62827, 2.06016, -1.22761, 1.79616, 0.640145, -0.812508,-1.50028, -1.54575, -0.217228, 2.06445, -0.0252311, 1.31984,-1.09519, -1.55624, 0.602962, 1.33464, 0.883919, 0.142988, -0.886204,0.426426, 0.893315, 1.30856, -1.24124, 0.239628, 1.96819, 1.91064,-1.39036, 1.02973, -0.120615, 1.24169, -1.40214, 1.03507, -1.46984,1.58711, 1.01167, 2.08498, 1.97875, 1.48291, -1.26303, -0.930141,-0.59183, -0.350317, 1.93194, 0.953808, 0.829422, 0.662458, 1.50385,-0.943274, 1.21817, 1.12564, 1.28737, 0.4419, -0.222186, 2.34941,1.18368, 1.9319, -0.194216, 3.35556, 2.80241, 2.51639, 0.537432,0.625744, 2.86291, 2.28998, 4.0202, 3.9159, 3.92103, 0.377672,2.28483, 2.40107, 1.61928, 0.578327, 4.38002, 0.645356, 3.7589,3.62083, 4.1443, 1.89589, 2.23305, 1.4046, 4.11185, 2.37104, 3.31836,2.83245];
This variable, y, contains 256 data points.
Adjust the default directory, i.e. 


cd c:/tmp
Method 1 


cd c:/tmp
% CD to the "c:/tmp" directory
%
y=[...............];
%
plot(y);
f1=fft(y); % Fast Fourier Transform of y, i.e. frequency domain expression of y
plot(abs(f1)); % Since f1 is a set of complex numbers, we take only the absolute values.
%
% Examine the graph of f1 (in frequency) and determine the cut-off values.
%
% Filtering out data with small amplitude
%
for i=1:256
 if abs(f1(i)) < 50
   f1(i)=0;
 end;
end;
%
% (alternative way)
% f1(abs(f1)<50)=0;
%
% Inverse FFT on f1
%
y2=ifft(f1);
%
% y2 is real but there is some dust so chop them off.
%
y2=real(y2);
plot(y2);
Method 2 


cd c:/tmp
% CD to the "c:/tmp" directory
%
y=[...............];
%
plot(y);
f1=fft(y);
plot(abs(f1));
%
% define filter
%
fil=zeros(1, 256);
fil(1, 1:10)=1;
fil(1,250:256)=1;
%
f2=fil.*f1;
%
plot(abs(f2));
y2=ifft(f2);
plot(real(y2));

Example 4

Wave file 
cd c:\tmp
x=audioread('cricket.wav');
[n, m]=size(x)

y=fft(x);
y2=abs(y);
plot(y2);
semilogx(y2);
% around 10^5
%
fil=zeros(n,m);
w1=40000;w2=100000;
fil(w1:w2, :)=1;
fil(n-w2:n-w1, :)=1;
%
y3=fil.*y;
x2=ifft(y3);
%
fs=44100;
%
audiowrite('cricket2.wav', real(x2), fs);
Note: Use wavread and wavwrite in Octave instead of audioread and audiowrite. Use wavwrite(ifft(y), 44100, 'newwavefile.wav');
Filtered wave file

Definition of DFT

  1. Our definition
    fl
    =
    N−1
    ∑
    k=0 
    		 Fk ei k [(2π)/(N)]l
    		(37)
    Fk
    =
    1
    N
    		 N−1
    ∑
    k=0 
    		 fl e−i k [(2π)/(N)]l
    		(38)
  2. MATLAB/Octave

    fl
    =
    1
    N
    		 N−1
    ∑
    k=0 
    		 Fk ei k [(2π)/(N)]l
    		(39)
    Fk
    =
    N−1
    ∑
    k=0 
    		 fl e−i k [(2π)/(N)]l
    		(40)
  3. Mathematica

    fl
    =
    1
    √N
    		 N−1
    ∑
    k=0 
    		 Fk e−i k [(2π)/(N)]l
    		(41)
    Fk
    =
    1
    √N
    		 N−1
    ∑
    k=0 
    		 fl ei k [(2π)/(N)]l
    		(42)
octave.exe:1> x=1:8
x =

   1   2   3   4   5   6   7   8

octave.exe:2> fft(x)
ans =

 Columns 1 through 4:

   36.0000 +  0.0000i   -4.0000 +  9.6569i   -4.0000 +  4.0000i   -4.0000 +  1.6569i

 Columns 5 through 8:

   -4.0000 +  0.0000i   -4.0000 -  1.6569i   -4.0000 -  4.0000i   -4.0000 -  9.6569i
fourier_mathematica.jpg

Footnotes:

1
f(x)= ∞
∑
m=−∞ 
 cm ei mx,     cm = 1
2π
 ⌠
⌡ 		2π

0 
 f(x)e−imx dx.
2 A more aggressive algorithm yields
M(N) = 2 M(N/2) + N/2−2,
which results in
M(N) = N
2
 log2 N − 3
2
 N + 2.


File translated from TEX by TTH, version 4.03.
On 07 Nov 2023, 21:17.