#19 (11/11/2024)

Note

Fourier coefficients (integrals) are pseudo-symmetric with respect to m or ω as

C_−m

C_m

(1)

F_−ω

F(ω)

(2)

Proof.

Just use the definition.

Power spectrum is symmetrical, i.e.

P(−ω) = P(ω).

(3)

Proof.

Just use the definition.

Nyquist-Shannon Sampling theorem

A Fourier transformable function, f(x), can be completely represented by its discrete values at x = …−3, 2, 1, 0, 1, 2, … by

f(x)

…+ f(−1)

sin π(x + 1)

π(x + 1)

+f(0)

sin πx

πx

+f(1)

sin π(x − 1)

π(x − 1)

+ …

∞
∑
k=−∞

f(k)

sin π(x−k)

π(x−k)

(4)

if f(x) does not contain higher frequency components than π, i.e.

F(ω) = 0 |ω| > π.

Proof As the bandwidth of f(x) is restricted between −π < ω < π, the inverse Fourier transform of f(x) is

f(x) =

2 π

⌠
⌡

∞

−∞

F(ω) e^{i ωx} dω =

2 π

⌠
⌡

π

−π

F(ω) e^{i ωx} dω,

(5)

where

F(ω) =

⌠
⌡

∞

−∞

f(x) e^−iωx dx.

As F(ω) exists only in the interval [−π, π], it can be expressed by the Fourier series as

F(ω) =

∞
∑
k=−∞

C_k e^{i k ω},

(6)

where

C_k

2 π

⌠
⌡

π

−π

F(ω) e^{−i k ω} dω

2 π

⌠
⌡

π

−π

F(ω) e^{i ω(−k)} dω

f(−k),

(7)

where we used Eq. (5). By substituting Eq. (7) into Eq.(6), one gets

F(ω) =

∞
∑
k=−∞

f(−k) e^{i k ω} =

∞
∑
k=−∞

f(k) e^{− i k ω}.

(8)

Finally, using Eq. (8), f(x) is expressed as

f(x)

2 π

⌠
⌡

π

−π

⎛
⎝

∞
∑
k=−∞

f(k) e^{−i k ω}

⎞
⎠

e^{i ωx} d ω

2 π

∞
∑
k=−∞

⎛
⎝

⌠
⌡

π

−π

e^{i ω(x − k)} dω

⎞
⎠

f(k)

2π

∞
∑
k=−∞

⎛
⎝

2 sin π(x − k)

(x − k)

⎞
⎠

f(k)

∞
∑
k=−∞

f(k)

sin π(x − k)

π(x − k)

(9)

This conclusion can be generalized as

f(x) =

∞
∑
k=−∞

⎛
⎝

⎞
⎠

sin

⎛
⎝

2 πW

⎛
⎝

x−

⎞
⎠

2 πW

⎛
⎝

x−

⎞
⎠

(10)

F(ω) = 0, |ω| ≥ W.

(11)

Derivation

Sampling Rate

Period

Definition

(−π, π)

(angular frequency)

2π

(−

)

(frequency)

(−W, W)

(frequency)

2 W

(−2 πW, 2 πW)

(angular frequency)

(12)

The sampling theorem indicates that if the bandwidth of f(x) is limited to [−W, W], f(x) can be completely reconstructed by sampling the value of f(x) with the interval of τ = 2 W.

Examples: Human ears can hear frequencies up to 22 kHz. Therefore, the CD sample rate is 44.1 kHz. Copper phone lines pass frequencies up to 4 kHz, hence, phone companies sample at 8 kHz.

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 expansion², (1) the infinite series must be truncated and (2) numerical integration must be used to compute the Fourier coefficient, c_m, 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.

The table below is the difference and similarity between CFT and DFT.

CFT

f(x) =

∞
∑
k=−∞

F_k e^{i k x}

F_k =

2 π

⌠
⌡

2 π

0

f(x) e^{−i k x}d x

DFT

f_l =

N−1
∑
k=0

F_k e^{i k [(2π)/(N)] l}

F_k =

N−1
∑
l = 0

f_l e^{−i k [(2π)/(N)]l}

(13)

Derivation

The derivation of the DFT coefficient, F_k, closely follows that of the CFT coefficient, c_k. 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

f_l

N−1
∑
k=0

F_k e^{i k [(2π)/(N)] l},

(14)

F_k

N−1
∑
l = 0

f_l e^{−i k [(2π)/(N)]l}.

(15)

Multiply Eq. (14) by e^{−i m [(2 π)/(N)] l} to obtain

f_l e^{−i m [(2 π)/(N)] l} =

N−1
∑
k=0

F_k e^{i k [(2 π)/(N)]l} e^{− i m [(2 π)/(N)]l}.

(16)

Take summation of Eq. (16) over l to get

N−1
∑
l = 0

f_l e^{− im [(2π)/(N)] l}

N−1
∑
l = 0

N−1
∑
k=0

F_k e^{i k [(2 π)/(N)]l} e^{− i m [(2 π)/(N)]l}

(interchange the order of summation)

N−1
∑
k=0

F_k

N−1
∑
l = 0

e^{i (k−m) [(2 π)/(N)] l}

⎛
⎜
⎜
⎝

Use

N−1
∑
l = 0

e^{i (k−m) [(2π)/(N)]l} =

⎧
⎪
⎨
⎪
⎩

k−m = 0, ±N, ±2N …

k−m ≠ 0, ±N, ±2 N …

⎞
⎟
⎟
⎠

N F_m.

(17)

Hence,

F_m =

N − 1
∑
l = 0

f_l e^{− i m [(2 π)/(N)] l}.

(18)

Note that in Eq. (17), we used an identity

N−1
∑
l = 0

e^{i k [(2π)/(N)] l} e^{−i m [(2π)/(N)] l} =

⎧
⎪
⎨
⎪
⎩

k ≠ m

k = m, m±N, m±2 N …

(19)

which is a discrete version of

⌠
⌡

π

−π

e^{i k x} e^{−i m x} d x =

⎧
⎪
⎨
⎪
⎩

k ≠ m

2 π

k = m

(20)

(Proof of Identity (19))

Define

ω ≡ e^{i (k−m) [(2 π)/(N)]}.

so that

N−1
∑
l = 0

e^{i k [(2π)/(N)] l} e^{−i m [(2π)/(N)] l}

N−1
∑
l = 0

ω^l

1 + ω+ ω² + ω³ + …ω^N−1.

(21)

k − m = 0, ±N, ±2 N, ±3N …,

ω = e⁰, e^{i [(2π)/(N)] N}, e^{i [(2π)/(N)]2 N} ... = 1 so that Eq. (21) becomes

N−1
∑
l = 0

e^{i k [(2π)/(N)] l} e^{−i m [(2π)/(N)] l} = 1 + 1 + …+1 = N.

(22)

On the other hand, if

k − m ≠ 0, ±N, ±2 N, ±3N …,

Eq. (21) becomes

1 + ω+ ω² + ω³ + …ω^N−1

1 − ω^N

1−ω

1− e^{i (k−m) 2 π}

1−e^{i (k−m) [(2π)/(N)]}

0 if k −m ≠ 0, ±N, ±2N …

(23)

Therefore, it was shown that

N−1
∑
l = 0

e^{i (k−m) [(2π)/(N)]l} =

⎧
⎪
⎨
⎪
⎩

k−m = 0, ±N, ±2N …

k−m ≠ 0, ±N, ±2 N …

(24)

Note that only half of F_k's are independent, i.e.

F_k =

F_N−k

(25)

The power spectrum, P_k, is symmetrical around its center, i.e.

P_k = P_N−k.

(26)

Relation between DFT and CFT

The classical Fourier series expansion of f(x) over [0, 2 π] is

f(x)

∞
∑
k=−∞

c_k e^{i k x} ,

(27)

c_k

2π

⌠
⌡

2π

0

f(x) e^{− i k x} dx.

(28)

If we consider an interval, 0 < x < 2 π, to be sampled at 0, [(2 π)/(N)], [(4 π)/(N)], [(6 π)/(N)] …, denoted as

f_l ≡ f

⎛
⎝

2 π

⎞
⎠

, ω ≡ e^{[(2 π)/(N)]i}.

(29)

Then, Eq. (15) becomes

F_k

N−1
∑
l = 0

f_l ω^{− k l}

N−1
∑
l = 0

⎛
⎝

∞
∑
k′=−∞

c_k′ ω^k′l

⎞
⎠

ω^{−k l}

∞
∑
k′=−∞

c_k′

⎛
⎝

N−1
∑
l = 0

ω^(k′−k)l

⎞
⎠

⎛
⎜
⎜
⎝

Use

N−1
∑
l = 0

e^{i (k−k′) [(2π)/(N)]l} =

⎧
⎪
⎨
⎪
⎩

k − k′ = 0, ±N, ±2N …

k − k′ ≠ 0, ±N, ±2 N …

⎞
⎟
⎟
⎠

∞
∑
k′=−∞

c_k′ ×(1 if and only if k′=k, k±N, k±2N, k±3N, …, otherwise 0)

c_k + c_k±N + c_{k±2 N} + c_k±3N +…

c_k +

∞
∑
l=1

c_{k±N l}.

(30)

Equation (30) 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 f_l and g_l is defined as

f_l * g_l ≡

N−1
∑
m=0

f_l−m g_m.

(31)

Note that

f_l * g_l = g_l * f_l.

(DFT Convolution Theorem)

f_l * g_l =

N−1
∑
k=0

F_k G_k e^{i [(2π)/(N)] k},

or DFT of f_l * g_l is F_k G_k.

Parseval's theorem of DFT

N−1
∑
l = 0

f_l

N−1
∑
k=0

F_k

(32)

N−1
∑
l = 0

||f_l||² =

N−1
∑
k=0

|| F_k||².

(33)

Power spectrum in DFT is defined as

P_k ≡ ||F_k||².

(34)

FFT

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

F_k =

N−1
∑
l = 0

f_l ω^−kl k = 0, 1, 2, …, (N−1), ω ≡ e^{2 πi/N}.

(35)

For each F_k, provided that f_l'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 N² 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 N² (multiplications) to N log₂ N. This reduction is dramatic as shown in the table below:

N	N²	N log₂ 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 f₀, f₁, f₂, …f_N−1, define

g(x) ≡ f₀ + f₁ x + f₂ x² + …+ f_N−1 x^N−1.

(36)

By comparing Eq. (36) with Eq. (35), it can be readily seen that the coefficient, F_k, can be computed by

F_k =

g(ω^−k).

(37)

To illustrate the algorithm of FFT, let N=8, then,

g(x)

f₀ + f₁ x + f₂ x² + f₃ x³ + f₄ x⁴ + f₅ x⁵ + f₆ x⁶ + f₇ x⁷

(f₀ + f₂ x² + f₄ x⁴ + f₆ x⁶) +x ( f₁ + f₃ x² + f₅ x⁴+ f₇ x⁶)

p(x²) + x q(x²),

(38)

where

⎧
⎪
⎨
⎪
⎩

p(x) ≡ f₀ + f₂ x + f₄ x² + f₆ x³

q(x) ≡ f₁ + f₃ x + f₅ x² + f₇ x³.

(39)

Thus, to compute all the coefficients of F_k is equivalent to computing both p(x²) and q(x²) at x=1, ω, ω², ω³, ω⁴, ω⁵, ω⁶, ω⁷, multiplying x by q(x²) and add them together. It should be noted that for

x=1, ω, ω², ω³, ω⁴, ω⁵, ω⁶, ω⁷,

(40)

the corresponding x²'s are (note ω⁸=1)

x²

1, ω², ω⁴, ω⁶, ω⁸, ω¹⁰, ω¹², ω¹⁴

1, ω², ω⁴, ω⁶, 1, ω², ω⁴, ω⁶.

(41)

It is seen that the second half of x² are the repetition of the first half, thus, it is NOT necessary to compute p(x²) and q(x²) when x belongs to the second half. Moreover, note that the set of {1, ω², ω⁴, ω⁶} is the roots to

x⁴=1,

(42)

so computing p(x²) and q(x²) is equivalent to computing p(x) and q(x) where x is the root for x⁴=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.

(43)

where the factor "2" comes from the computation of p(x²) and q(x²), the argument, N/2, of M comes from the equation of x⁴ = 1 (x^(N/2)=1) and the last term of N is for the multiplication between x and q(x²) for N times.

When N=1, g(x)=f₀, thus,

M(0)=1

(44)

From Eqs.(43-44), the solution to Eq. (43) is

M(N ) = N log₂ N.³

(45)

Similarly, the number of additions can be computed as

A (N) = N log₂ N.

(46)

Example 1

Consider a data set of (N=16)

f₁ = {0, 1, 2, 3, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }

(47)

which looks like

The power spectrum, |C_k|², of f₁ is

Now consider

f₂={0, 0, 0, 0, 0, 0, 1, 2, 3, 2, 1, 0, 0, 0, 0, 0}

(48)

The power spectrum, |C_k|², of f₂ is

Example 2 (2-D)

Consider two different images

The power spectra of the images are shown as

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

Our definition

f_l

=

N−1
∑
k=0
F_k e^{i k [(2π)/(N)]l}
(49)

F_k

=

1
N
N−1
∑
k=0
f_l e^{−i k [(2π)/(N)]l}
(50)
MATLAB/Octave

f_l

=

1
N
N−1
∑
k=0
F_k e^{i k [(2π)/(N)]l}
(51)

F_k

=

N−1
∑
k=0
f_l e^{−i k [(2π)/(N)]l}
(52)
Mathematica

f_l

=

1
√N
N−1
∑
k=0
F_k e^{−i k [(2π)/(N)]l}
(53)

F_k

=

1
√N
N−1
∑
k=0
f_l e^{i k [(2π)/(N)]l}
(54)

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

Footnotes:

ω = 2 πf

where ω is an angular frequency (unit: rad/time) and f is a (regular) frequency (unit: Hz).

f(x)=

∞
∑
m=−∞

c_m e^{i mx}, c_m =

2π

⌠
⌡

2π

0

f(x)e^−imx dx.

³ A more aggressive algorithm yields

M(N) = 2 M(N/2) + N/2−2,

which results in

M(N) =

log₂ N −

N + 2.

File translated from T_EX by T_TH, version 4.03.
On 15 Nov 2024, 12:36.