In this post we will talk about Fourier Transform.
Fourier series assumes that the signal is periodic. But what do we do if the isignal is aperiodic? A trick is to extend the fundamental period from T 0 T_{0} T 0 ∞ \infty ∞
Recall the fourier series coefficients for a continous-time signal:
C k = 1 T 0 ∫ − T 0 2 T 0 2 f ( t ) e − i ω 0 k t d t \begin{aligned}
C_{k} &= \frac{1}{T_{0}}\int_{\frac{-T_{0}}{2}}^{\frac{T_{0}}{2}}f(t)e^{-i\omega_{0} kt}\ dt
\end{aligned} C k = T 0 1 ∫ 2 − T 0 2 T 0 f ( t ) e − i ω 0 k t d t The period as a function of frequency ω 0 \omega_{0} ω 0
ω 0 = 2 π T 0 2 T 0 = ω 0 2 π \begin{aligned}
\omega_{0} &= \frac{2\pi}{T_{0}}\\
\frac{2}{T_{0}} &= \frac{\omega_{0}}{2\pi}
\end{aligned} ω 0 T 0 2 = T 0 2 π = 2 π ω 0 And taking the limit to ∞ \infty ∞ ω 0 \omega_{0} ω 0 Δ ω \Delta \omega Δ ω
lim T 0 → ∞ 2 T 0 ≈ Δ ω 2 π \begin{aligned}
\lim_{T_{0} \to \infty}\frac{2}{T_{0}} &\approx \frac{\Delta \omega}{2\pi}
\end{aligned} T 0 → ∞ lim T 0 2 ≈ 2 π Δ ω We re-write C k C_{k} C k
C k = Δ ω 2 π ∫ − T 0 2 T 0 2 f ( t ) e − i ω 0 k t d t \begin{aligned}
C_{k} &= \frac{\Delta \omega}{2\pi}\int_{\frac{-T_{0}}{2}}^{\frac{T_{0}}{2}}f(t)e^{-i\omega_{0} kt}\ dt
\end{aligned} C k = 2 π Δ ω ∫ 2 − T 0 2 T 0 f ( t ) e − i ω 0 k t d t Inserting this into the fourier series:
x ( t ) = ∑ k = − ∞ ∞ { Δ ω 2 π ∫ − T 0 2 T 0 2 f ( t ) e − i ω 0 k t d t } e i k ω 0 t \begin{aligned}
x(t) &= \sum_{k = -\infty}^{\infty}\Big\{ \frac{\Delta \omega}{2\pi}\int_{\frac{-T_{0}}{2}}^{\frac{T_{0}}{2}}f(t)e^{-i\omega_{0} kt}\ dt\Big\}e^{ik\omega_{0}t}
\end{aligned} x ( t ) = k = − ∞ ∑ ∞ { 2 π Δ ω ∫ 2 − T 0 2 T 0 f ( t ) e − i ω 0 k t d t } e i k ω 0 t We re-write the equation as:
x ( t ) = 1 2 π ∑ k = − ∞ ∞ { ∫ − ∞ ∞ f ( t ) e − i ω 0 k t d t } e i k ω t d ω = 1 2 π ∑ k = − ∞ ∞ X ( ω ) e i k ω t d ω \begin{aligned}
x(t) &= \frac{1}{2\pi}\sum_{k = -\infty}^{\infty}\Big\{ \int_{-\infty}^{\infty}f(t)e^{-i\omega_{0} kt}\ dt\Big\}e^{ik\omega t}\ d\omega\\
&= \frac{1}{2\pi}\sum_{k = -\infty}^{\infty}X(\omega)e^{ik\omega t}\ d\omega
\end{aligned} x ( t ) = 2 π 1 k = − ∞ ∑ ∞ { ∫ − ∞ ∞ f ( t ) e − i ω 0 k t d t } e i k ω t d ω = 2 π 1 k = − ∞ ∑ ∞ X ( ω ) e i k ω t d ω Essentially we have broken up the signal into its frequency components as weights and multiplying them by the orthogonal basis function e i k ω t e^{ik\omega t} e i k ω t X ( ω ) X(\omega) X ( ω ) i s i n ( ω ) t i sin(\omega)t i s i n ( ω ) t
X ( ω ) X(\omega) X ( ω )
X ( ω ) = ∫ − ∞ ∞ f ( t ) e − i ω t d t \begin{aligned}
X(\omega) &= \int_{-\infty}^{\infty}f(t)e^{-i\omega t}\ dt
\end{aligned} X ( ω ) = ∫ − ∞ ∞ f ( t ) e − i ω t d t And the following is the Inverse Fourier Transform:
x ( t ) = 1 2 π ∑ k = − ∞ ∞ X ( ω ) e i k ω t d ω \begin{aligned}
x(t) &= \frac{1}{2\pi}\sum_{k = -\infty}^{\infty}X(\omega)e^{ik\omega t}\ d\omega
\end{aligned} x ( t ) = 2 π 1 k = − ∞ ∑ ∞ X ( ω ) e i k ω t d ω We can also write out the Fourier Transform in terms of its magnitude and phase:
X ( ω ) = ∣ X ( ω ) ∣ e i ∠ X ( ω ) \begin{aligned}
X(\omega) &= |X(\omega)|e^{i \angle X(\omega)}
\end{aligned} X ( ω ) = ∣ X ( ω ) ∣ e i ∠ X ( ω ) The whole derivation above is to show that it is possible to do a Fourier Transform on an aperiodic signal.
Dirac defined the delta function as a sum of infinite number of exponentials:
δ ( t ) = 1 2 π ∫ − ∞ ∞ e i ω t d ω \begin{aligned}
\delta(t) &= \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{i\omega t}\ d\omega
\end{aligned} δ ( t ) = 2 π 1 ∫ − ∞ ∞ e i ω t d ω And in the frequency domain:
δ ( t ) = 1 2 π ∫ − ∞ ∞ e i ω t d t \begin{aligned}
\delta(t) &= \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{i\omega t}\ dt
\end{aligned} δ ( t ) = 2 π 1 ∫ − ∞ ∞ e i ω t d t The general version with time/frequency shift:
δ ( t − t 0 ) = 1 2 π ∫ − ∞ ∞ e i ω ( t − t 0 ) d ω δ ( ω − ω 0 ) = 1 2 π ∫ − ∞ ∞ e i ( ω − ω 0 ) t d t \begin{aligned}
\delta(t - t_{0}) &= \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{i\omega (t - t_{0})}\ d\omega\\[5pt]
\delta(\omega - \omega_{0}) &= \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{i(\omega - \omega_{0})t}\ dt
\end{aligned} δ ( t − t 0 ) δ ( ω − ω 0 ) = 2 π 1 ∫ − ∞ ∞ e i ω ( t − t 0 ) d ω = 2 π 1 ∫ − ∞ ∞ e i ( ω − ω 0 ) t d t A periodic signal can be represented in FSC (Fourier Series Coefficients) form:
x ( t ) = ∑ k = − ∞ ∞ C k e i ω 0 t \begin{aligned}
x(t) &= \sum_{k = -\infty}^{\infty}C_{k}e^{i\omega_{0}t}
\end{aligned} x ( t ) = k = − ∞ ∑ ∞ C k e i ω 0 t Taking the fourier transform:
X ( ω ) = F { x ( t ) } = F { ∑ k = − ∞ ∞ C k e i ω 0 t } = ∑ k = − ∞ ∞ C k F { e i ω 0 t } \begin{aligned}
X(\omega) &= \mathcal{F}\{x(t)\}\\
&= \mathcal{F}\Big\{\sum_{k = -\infty}^{\infty}C_{k}e^{i\omega_{0}t}\Big\}\\
&= \sum_{k = -\infty}^{\infty}C_{k}\mathcal{F}\{e^{i\omega_{0}t}\}\\
\end{aligned} X ( ω ) = F { x ( t ) } = F { k = − ∞ ∑ ∞ C k e i ω 0 t } = k = − ∞ ∑ ∞ C k F { e i ω 0 t } It boils down to taking the fourier transform of an exponential signal:
e i ω 0 t = c o s ( ω 0 t ) + i s i n ( ω 0 t ) X ( e i ω 0 t ) = X ( c o s ( ω 0 t ) ) + i X ( s i n ( ω 0 t ) ) \begin{aligned}
e^{i\omega_{0}t} &= \mathrm{cos}(\omega_{0}t) + i \mathrm{sin}(\omega_{0}t)\\
X(e^{i\omega_{0}t}) &= X(\mathrm{cos}(\omega_{0}t)) + i X(\mathrm{sin}(\omega_{0}t))\\
\end{aligned} e i ω 0 t X ( e i ω 0 t ) = c o s ( ω 0 t ) + i s i n ( ω 0 t ) = X ( c o s ( ω 0 t ) ) + i X ( s i n ( ω 0 t ) ) Fourier transform of a cosine:
X ( ω ) = F ( c o s ( ω 0 t ) ) = ∫ − ∞ ∞ 1 2 ( e i ω 0 t + e − i ω 0 t ) e i ω t d t = ∫ − ∞ ∞ 1 2 e i ( ω 0 + ω 0 ) t d t + ∫ − ∞ ∞ 1 2 e i ( ω 0 − ω 0 ) t d t \begin{aligned}
X(\omega) &= \mathcal{F}(\mathrm{cos}(\omega_{0}t))\\
&= \int_{-\infty}^{\infty}\frac{1}{2}(e^{i\omega_{0}t} + e^{-i\omega_{0}t})e^{i\omega t}\ dt\\
&= \int_{-\infty}^{\infty}\frac{1}{2}e^{i(\omega_{0} + \omega_{0})t}\ dt +
\int_{-\infty}^{\infty}\frac{1}{2}e^{i(\omega_{0} - \omega_{0})t}\ dt
\end{aligned} X ( ω ) = F ( c o s ( ω 0 t ) ) = ∫ − ∞ ∞ 2 1 ( e i ω 0 t + e − i ω 0 t ) e i ω t d t = ∫ − ∞ ∞ 2 1 e i ( ω 0 + ω 0 ) t d t + ∫ − ∞ ∞ 2 1 e i ( ω 0 − ω 0 ) t d t As discussed in the Delta Function section, we can see that:
δ ( ω − ω 0 ) = 1 2 π ∫ − ∞ ∞ e i ( ω − ω 0 ) t d t π δ ( ω − ω 0 ) = 1 2 ∫ − ∞ ∞ e i ( ω − ω 0 ) t d t \begin{aligned}
\delta(\omega - \omega_{0}) &= \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{i(\omega - \omega_{0})t}\ dt\\
\pi\delta(\omega - \omega_{0}) &= \frac{1}{2}\int_{-\infty}^{\infty}e^{i(\omega - \omega_{0})t}\ dt\\
\end{aligned} δ ( ω − ω 0 ) π δ ( ω − ω 0 ) = 2 π 1 ∫ − ∞ ∞ e i ( ω − ω 0 ) t d t = 2 1 ∫ − ∞ ∞ e i ( ω − ω 0 ) t d t Hence:
F { c o s ( ω 0 t ) } = ∫ − ∞ ∞ 1 2 e i ( ω + ω 0 ) t d t + ∫ − ∞ ∞ 1 2 e i ( ω − ω 0 ) t d t = π δ ( ω + ω 0 ) + π δ ( ω − ω 0 ) \begin{aligned}
\mathcal{F}\{cos(\omega_{0}t)\} &= \int_{-\infty}^{\infty}\frac{1}{2}e^{i(\omega + \omega_{0})t}\ dt +
\int_{-\infty}^{\infty}\frac{1}{2}e^{i(\omega - \omega_{0})t}\ dt\\
&= \pi\delta(\omega + \omega_{0}) + \pi\delta(\omega - \omega_{0})
\end{aligned} F { c o s ( ω 0 t ) } = ∫ − ∞ ∞ 2 1 e i ( ω + ω 0 ) t d t + ∫ − ∞ ∞ 2 1 e i ( ω − ω 0 ) t d t = π δ ( ω + ω 0 ) + π δ ( ω − ω 0 ) Similarly, for sine:
X ( ω ) = F { s i n ( ω 0 t ) } = i π δ ( ω + ω 0 ) − i π δ ( ω − ω 0 ) \begin{aligned}
X(\omega) &= \mathcal{F}\{\mathrm{sin}(\omega_{0}t)\}\\
&= i\pi\delta(\omega + \omega_{0}) - i\pi\delta(\omega - \omega_{0})
\end{aligned} X ( ω ) = F { s i n ( ω 0 t ) } = i π δ ( ω + ω 0 ) − i π δ ( ω − ω 0 ) Using the above result:
X ( e i ω 0 t ) = X ( c o s ( ω 0 t ) ) + i X ( s i n ( ω 0 t ) ) = π δ ( ω + ω 0 ) + π δ ( ω − ω 0 ) + i 2 ( π δ ( ω + ω 0 ) − i 2 π δ ( ω − ω 0 ) ) = π δ ( ω − ω 0 ) − i 2 π δ ( ω − ω 0 ) + π δ ( ω + ω 0 ) + i 2 π δ ( ω + ω 0 ) = π δ ( ω − ω 0 ) + π δ ( ω − ω 0 ) + π δ ( ω + ω 0 ) − π δ ( ω + ω 0 ) = 2 π δ ( ω − ω 0 ) \begin{aligned}
X(e^{i\omega_{0}t}) &= X(\mathrm{cos}(\omega_{0}t)) + i X(\mathrm{sin}(\omega_{0}t))\\
&= \pi\delta(\omega+\omega_{0}) + \pi\delta(\omega - \omega_{0}) + i^{2}(\pi\delta(\omega + \omega_{0}) - i^{2}\pi\delta(\omega-\omega_{0}))\\
&= \pi\delta(\omega-\omega_{0}) - i^{2}\pi\delta(\omega - \omega_{0}) +
\pi\delta(\omega + \omega_{0}) + i^{2}\pi\delta(\omega + \omega_{0})\\
&= \pi\delta(\omega-\omega_{0}) + \pi\delta(\omega - \omega_{0}) +
\pi\delta(\omega + \omega_{0}) - \pi\delta(\omega + \omega_{0})\\
&= 2\pi\delta(\omega - \omega_{0})
\end{aligned} X ( e i ω 0 t ) = X ( c o s ( ω 0 t ) ) + i X ( s i n ( ω 0 t ) ) = π δ ( ω + ω 0 ) + π δ ( ω − ω 0 ) + i 2 ( π δ ( ω + ω 0 ) − i 2 π δ ( ω − ω 0 ) ) = π δ ( ω − ω 0 ) − i 2 π δ ( ω − ω 0 ) + π δ ( ω + ω 0 ) + i 2 π δ ( ω + ω 0 ) = π δ ( ω − ω 0 ) + π δ ( ω − ω 0 ) + π δ ( ω + ω 0 ) − π δ ( ω + ω 0 ) = 2 π δ ( ω − ω 0 ) Back to the fourier transform of a periodic signal:
X ( ω ) = F { x ( t ) } = F { ∑ k = − ∞ ∞ C k e i ω 0 t } = ∑ k = − ∞ ∞ C k F { e i ω 0 t } = 2 π ∑ k = − ∞ ∞ C k δ ( ω − k ω 0 ) \begin{aligned}
X(\omega) &= \mathcal{F}\{x(t)\}\\
&= \mathcal{F}\Big\{\sum_{k = -\infty}^{\infty}C_{k}e^{i\omega_{0}t}\Big\}\\
&= \sum_{k = -\infty}^{\infty}C_{k}\mathcal{F}\{e^{i\omega_{0}t}\}\\
&= 2\pi \sum_{k = -\infty}^{\infty}C_{k}\delta(\omega - k\omega_{0})
\end{aligned} X ( ω ) = F { x ( t ) } = F { k = − ∞ ∑ ∞ C k e i ω 0 t } = k = − ∞ ∑ ∞ C k F { e i ω 0 t } = 2 π k = − ∞ ∑ ∞ C k δ ( ω − k ω 0 ) This is essentially the same as the discrete fourier series coefficient of the signal but scaled to be 2 π 2\pi 2 π X ( ω ) X(\omega) X ( ω ) 2 π 2\pi 2 π
Up till now, what we have been referring to is the Continuous Fourier Transfor:
X ( ω ) = ∫ − ∞ ∞ f ( t ) e − i ω t d t \begin{aligned}
X(\omega) &= \int_{-\infty}^{\infty}f(t)e^{-i\omega t}\ dt
\end{aligned} X ( ω ) = ∫ − ∞ ∞ f ( t ) e − i ω t d t For DFT:
X ( k ) = ∑ n = 0 N − 1 x ( n ) e − i ω k n ω k = 2 π k N k = { 0 , ⋯ , N − 1 } \begin{aligned}
X(k) &= \sum_{n = 0}^{N - 1}x(n)e^{-i\omega_{k} n}\\[5pt]
\omega_{k} &= 2\pi \frac{k}{N}\\[5pt]
k &= \{0, \cdots, N-1\}
\end{aligned} X ( k ) ω k k = n = 0 ∑ N − 1 x ( n ) e − i ω k n = 2 π N k = { 0 , ⋯ , N − 1 } Instead of t, we replaced it with n n n ω k \omega_{k} ω k N N N 2 π 2\pi 2 π
If you look closely at the DFT equation, you could see that it is essentially finding the correlation of the N N N N N N e − i ω k n e^{-i\omega_{k} n} e − i ω k n
For IDFT (Inverse Discrete FT):
x n = 1 N ∑ k = 0 N − 1 X k e + i ω k n ω k = 2 π k N n = { 0 , ⋯ , N − 1 } \begin{aligned}
x_{n} &= \frac{1}{N}\sum_{k = 0}^{N-1}X_{k}e^{+i \omega_{k}n}\\[5pt]
\omega_{k} &= 2\pi \frac{k}{N}\\[5pt]
n &= \{0, \cdots, N-1\}
\end{aligned} x n ω k n = N 1 k = 0 ∑ N − 1 X k e + i ω k n = 2 π N k = { 0 , ⋯ , N − 1 } So N N N N N N N N N N N N
Note that there are N N N k k k N N N N 2 N^{2} N 2 N N N N l o g ( N ) N \mathrm{log}(N) N l o g ( N )
Let us take a cosine wave with the following info:
f = 1 f s = 16 f n = 16 2 = 8 \begin{aligned}
f &= 1\\
f_{s} &= 16\\
f_{n} &= \frac{16}{2}\\
&= 8
\end{aligned} f f s f n = 1 = 1 6 = 2 1 6 = 8 In other words, the frequency of the cosine curve is 1Hz. Sampling requency f s f_{s} f s f n f_{n} f n
Using R FFT routine to verify the result, we get the following result:
fs <- 8
fn <- fs / 2
f <- 1
omega <- 2 * pi * f
n <- seq ( 0 , 1 , length = fs )
xn <- cos ( omega * n [ 1 : length ( n )])
fft ( xn )
> xn
[ 1 ] 1.0000000
[ 2 ] 0.6234898
[ 3 ] -0.2225209
[ 4 ] -0.9009689
[ 5 ] -0.9009689
[ 6 ] -0.2225209
[ 7 ] 0.6234898
[ 8 ] 1.0000000
> fft ( xn )
[ 1 ] 0.9999999999999990007993+0.0000000i
[ 2 ] 3.8433767734721757669547+1.5919788i
[ 3 ] -0.3019377358048381254640-0.3019377i
[ 4 ] -0.0414390376673379190464-0.1000427i
[ 5 ] 0.0000000000000007771561+0.0000000i
[ 6 ] -0.0414390376673376970018+0.1000427i
[ 7 ] -0.3019377358048382364863+0.3019377i
[ 8 ] 3.8433767734721762110439-1.5919788i We can see index 2 has the highest value. Because we start with n = 0 n = 0 n = 0 k = 1 k = 1 k = 1
index 1 = DC Offset index 2 = 1 Hz index 3 = 2 Hz index 4 = 3 Hz index 5 = 4 Hz index 5 = 3 Hz index 5 = 2 Hz index 5 = 1 Hz \begin{aligned}
\text{index 1} &= \text{DC Offset}\\
\text{index 2} &= \text{1 Hz}\\
\text{index 3} &= \text{2 Hz}\\
\text{index 4} &= \text{3 Hz}\\
\text{index 5} &= \text{4 Hz}\\
\text{index 5} &= \text{3 Hz}\\
\text{index 5} &= \text{2 Hz}\\
\text{index 5} &= \text{1 Hz}
\end{aligned} index 1 index 2 index 3 index 4 index 5 index 5 index 5 index 5 = DC Offset = 1 Hz = 2 Hz = 3 Hz = 4 Hz = 3 Hz = 2 Hz = 1 Hz Anything higher than the Nyquist rate will result in aliasing, so the result is repeated but reflected.
We will do a calculation for index 2:
ω k = 1 = 2 π 1 16 X ( k = 1 ) = ∑ n = 0 N − 1 x ( n ) e − i ω k n = 1.0000000 × e − i ω 1 0 + 0.6234898 × e − i ω 1 1 + − 0.2225209 × e − i ω 1 2 + − 0.9009689 × e − i ω 1 3 + − 0.9009689 × e − i ω 1 4 + − 0.2225209 × e − i ω 1 5 + 0.6234898 × e − i ω 1 6 + 1.0000000 × e − i ω 1 7 \begin{aligned}
\omega_{k = 1} = &2\pi\frac{1}{16}\\
X(k = 1) = &\sum_{n = 0}^{N-1}x(n)e^{-i\omega_{k}n}\\
=\ &1.0000000 \times e^{-i\omega_{1}0} +
0.6234898 \times e^{-i\omega_{1}1} +\\
&-0.2225209 \times e^{-i\omega_{1}2} +
-0.9009689 \times e^{-i\omega_{1}3} +\\
&-0.9009689 \times e^{-i\omega_{1}4} +
-0.2225209 \times e^{-i\omega_{1}5} +\\
&0.6234898 \times e^{-i\omega_{1}6} +
1.0000000 \times e^{-i\omega_{1}7}
\end{aligned} ω k = 1 = X ( k = 1 ) = = 2 π 1 6 1 n = 0 ∑ N − 1 x ( n ) e − i ω k n 1 . 0 0 0 0 0 0 0 × e − i ω 1 0 + 0 . 6 2 3 4 8 9 8 × e − i ω 1 1 + − 0 . 2 2 2 5 2 0 9 × e − i ω 1 2 + − 0 . 9 0 0 9 6 8 9 × e − i ω 1 3 + − 0 . 9 0 0 9 6 8 9 × e − i ω 1 4 + − 0 . 2 2 2 5 2 0 9 × e − i ω 1 5 + 0 . 6 2 3 4 8 9 8 × e − i ω 1 6 + 1 . 0 0 0 0 0 0 0 × e − i ω 1 7 Using R to calculate the value:
omega_k <- 2 * pi * 1 / fs
> 1.0000000 * exp ( -1i * omega_k * 0 ) +
0.6234898 * exp ( -1i * omega_k * 1 ) +
-0.2225209 * exp ( -1i * omega_k * 2 ) +
-0.9009689 * exp ( -1i * omega_k * 3 ) +
-0.9009689 * exp ( -1i * omega_k * 4 ) +
-0.2225209 * exp ( -1i * omega_k * 5 ) +
0.6234898 * exp ( -1i * omega_k * 6 ) +
1.0000000 * exp ( -1i * omega_k * 7 )
[ 1 ] 3.843377+1.591979i We can see it agrees with the value we got from FFT.
FFT is a recursive divide and conquer algorithm that reduces the time complexity of Fourier Transform from O ( N 2 ) O(N^{2}) O ( N 2 ) O ( N l o g ( N ) ) O\Big(N\mathrm{log}(N)\Big) O ( N l o g ( N ) ) X ( ω ) X(\omega) X ( ω )
Firstly, let us do the Fourier Transform the long way:
X ( k ) = ∑ n = 0 N − 1 x ( n ) e − i ω k n k = { 0 , ⋯ , N − 1 } \begin{aligned}
X(k) &= \sum_{n = 0}^{N - 1}x(n)e^{-i\omega_{k} n}\\
k &= \{0, \cdots, N - 1 \}
\end{aligned} X ( k ) k = n = 0 ∑ N − 1 x ( n ) e − i ω k n = { 0 , ⋯ , N − 1 } res <- c ()
N <- length ( xn )
for ( k in 0 : ( N - 1 )) {
res2 <- 0
for ( j in 0 : ( N - 1 )) {
res2 <- res2 + xn [ j + 1 ] * exp ( -1i * 2 * pi * k * j / N )
}
res <- c ( res , res2 )
}
> res
[ 1 ] 1.000000e+00+0.000000e+00i 3.843377e+00+1.591979e+00i
[ 3 ] -3.019377e-01-3.019377e-01i -4.143904e-02-1.000427e-01i
[ 5 ] 8.881784e-16-5.040594e-16i -4.143904e-02+1.000427e-01i
[ 7 ] -3.019377e-01+3.019377e-01i 3.843377e+00-1.591979e+00i Next we divide the series into 2 parts. Even 2 m 2m 2 m 2 m + 1 2m+1 2 m + 1
X ( k ) = ∑ n = 0 N − 1 x ( n ) e − i ω k n = ∑ m = 0 N 2 − 1 x ( 2 m ) e − i ω k 2 m + ∑ m = 0 N 2 − 1 x ( 2 m + 1 ) e − i ω k ( 2 m + 1 ) = ∑ m = 0 N 2 − 1 x ( 2 m ) e − i ω k 2 m + e − i ω k ∑ m = 0 N 2 − 1 x ( 2 m + 1 ) e − i ω k 2 m = ∑ m = 0 N 2 − 1 x ( 2 m ) e − i 2 π k N 2 m + e − i 2 π k N ∑ m = 0 N 2 − 1 x ( 2 m + 1 ) e − i 2 π k N 2 m = ∑ m = 0 N 2 − 1 x ( 2 m ) e − i 2 π k m N 2 + e − i 2 π k N ∑ m = 0 N 2 − 1 x ( 2 m + 1 ) e − i 2 π k m N 2 \begin{aligned}
X(k) &= \sum_{n = 0}^{N - 1}x(n)e^{-i\omega_{k} n}\\
&= \sum_{m = 0}^{\frac{N}{2} - 1}x(2m)e^{-i\omega_{k} 2m} +
\sum_{m = 0}^{\frac{N}{2} - 1}x(2m + 1)e^{-i\omega_{k} (2m + 1)}\\
&= \sum_{m = 0}^{\frac{N}{2} - 1}x(2m)e^{-i\omega_{k} 2m} +
e^{-i\omega_{k}}\sum_{m = 0}^{\frac{N}{2} - 1}x(2m + 1)e^{-i\omega_{k} 2m}\\
&= \sum_{m = 0}^{\frac{N}{2} - 1}x(2m)e^{-i 2\pi \frac{k}{N} 2m} +
e^{-i 2\pi \frac{k}{N}}\sum_{m = 0}^{\frac{N}{2} - 1}x(2m + 1)e^{-i 2\pi \frac{k}{N} 2m}\\
&= \sum_{m = 0}^{\frac{N}{2} - 1}x(2m)e^{-i 2\pi \frac{km}{\frac{N}{2}}} +
e^{-i 2\pi \frac{k}{N}}\sum_{m = 0}^{\frac{N}{2} - 1}x(2m + 1)e^{-i 2\pi \frac{km}{\frac{N}{2}}}
\end{aligned} X ( k ) = n = 0 ∑ N − 1 x ( n ) e − i ω k n = m = 0 ∑ 2 N − 1 x ( 2 m ) e − i ω k 2 m + m = 0 ∑ 2 N − 1 x ( 2 m + 1 ) e − i ω k ( 2 m + 1 ) = m = 0 ∑ 2 N − 1 x ( 2 m ) e − i ω k 2 m + e − i ω k m = 0 ∑ 2 N − 1 x ( 2 m + 1 ) e − i ω k 2 m = m = 0 ∑ 2 N − 1 x ( 2 m ) e − i 2 π N k 2 m + e − i 2 π N k m = 0 ∑ 2 N − 1 x ( 2 m + 1 ) e − i 2 π N k 2 m = m = 0 ∑ 2 N − 1 x ( 2 m ) e − i 2 π 2 N k m + e − i 2 π N k m = 0 ∑ 2 N − 1 x ( 2 m + 1 ) e − i 2 π 2 N k m In each step, we break the equation into two parts and each part is half the size of the parent. We can then recursively reduce the series by half and end up with only 1 element and:
x ( 0 ) e − i ω k × 0 = x ( 0 ) \begin{aligned}
x(0)e^{-i\omega_{k}\times 0} &= x(0)
\end{aligned} x ( 0 ) e − i ω k × 0 = x ( 0 ) Note that the factor e i ω k e^{i\omega_{k}} e i ω k m m m N − 1 N - 1 N − 1
Below is the algorithm implemented in R and using the earlier cosine example:
fft <- function ( x , k = 1 ) {
n <- length ( x )
if ( n == 1 ) {
return ( x )
}
even <- fft3 ( x [ seq ( 1 , n , 2 )], k )
odd <- fft3 ( x [ seq ( 2 , n , 2 )], k )
factor <- exp ( -1i * 2 * pi / n * k )
result <- even + odd * factor
result
}
> for ( k in 1 : length ( xn )) {
+ print ( fft ( xn , k ))
+ }
[ 1 ] 3.843377+1.591979i
[ 1 ] -0.3019377-0.3019377i
[ 1 ] -0.041439-0.1000427i
[ 1 ] 7.771561e-16-5.040594e-16i
[ 1 ] -0.041439+0.1000427i
[ 1 ] -0.3019377+0.3019377i
[ 1 ] 3.843377-1.591979i
[ 1 ] 1+0i Let us do all k k k
fft <- function ( x ) {
n <- length ( x )
if ( n == 1 ) {
return ( x )
}
even <- fft ( x [ seq ( 1 , n , 2 )])
odd <- fft ( x [ seq ( 2 , n , 2 )])
factor1 <- exp ( -1i * 2 * pi / n * seq ( 0 , n / 2 - 1 ))
factor2 <- exp ( -1i * 2 * pi / n * seq ( n / 2 , n - 1 ))
result <- c ( even + factor1 * odd , even + factor2 * odd )
result
}
> fft ( xn )
[ 1 ] 1.000000e+00+0.000000e+00i 3.843377e+00+1.591979e+00i
[ 3 ] -3.019377e-01-3.019377e-01i -4.143904e-02-1.000427e-01i
[ 5 ] 7.771561e-16-6.123234e-17i -4.143904e-02+1.000427e-01i
[ 7 ] -3.019377e-01+3.019377e-01i 3.843377e+00-1.591979e+00i Note we have 2 factors in the function.
Because we have broken down each even and odd by half, they each would have length N 2 \frac{N}{2} 2 N ± N 2 \pm\frac{N}{2} ± 2 N
X ( k + N 2 ) = ∑ m = 0 N 2 − 1 x ( 2 m ) e − i ω k + N 2 2 m + e − i ω k + N 2 ∑ m = 0 N 2 − 1 x ( 2 m + 1 ) e − i ω k + N 2 2 m = ∑ m = 0 N 2 − 1 x ( 2 m ) e − i ω k 2 m + e − i ω k + N 2 ∑ m = 0 N 2 − 1 x ( 2 m + 1 ) e − i ω k 2 m \begin{aligned}
X\Big(k + \frac{N}{2}\Big) &= \sum_{m = 0}^{\frac{N}{2} - 1}x(2m)e^{-i\omega_{k + \frac{N}{2}} 2m} +
e^{-i\omega_{k + \frac{N}{2}}} \sum_{m = 0}^{\frac{N}{2} - 1}x(2m + 1)e^{-i\omega_{k + \frac{N}{2}} 2m}\\
&= \sum_{m = 0}^{\frac{N}{2} - 1}x(2m)e^{-i\omega_{k} 2m} +
e^{-i\omega_{k + \frac{N}{2}}} \sum_{m = 0}^{\frac{N}{2} - 1}x(2m + 1)e^{-i\omega_{k} 2m}
\end{aligned} X ( k + 2 N ) = m = 0 ∑ 2 N − 1 x ( 2 m ) e − i ω k + 2 N 2 m + e − i ω k + 2 N m = 0 ∑ 2 N − 1 x ( 2 m + 1 ) e − i ω k + 2 N 2 m = m = 0 ∑ 2 N − 1 x ( 2 m ) e − i ω k 2 m + e − i ω k + 2 N m = 0 ∑ 2 N − 1 x ( 2 m + 1 ) e − i ω k 2 m Since:
e − i ω k + N 2 = e − i ω k − i ω N 2 = e − i ω k − i 2 π N 2 N = e − i ω k e − i π = − 1 × e − i ω k = − e − i ω k \begin{aligned}
e^{-i\omega_{k + \frac{N}{2}}} &= e^{-i\omega_{k} - i\omega_{\frac{N}{2}}}\\
&= e^{-i\omega_{k} - i 2\pi \frac{N}{2N}}\\
&= e^{-i\omega_{k}}e^{- i \pi}\\
&= -1 \times e^{-i\omega_{k}}\\
&= -e^{-i\omega_{k}}
\end{aligned} e − i ω k + 2 N = e − i ω k − i ω 2 N = e − i ω k − i 2 π 2 N N = e − i ω k e − i π = − 1 × e − i ω k = − e − i ω k X ( k + N 2 ) = ∑ m = 0 N 2 − 1 x ( 2 m ) e − i ω k 2 m + e − i ω k + N 2 ∑ m = 0 N 2 − 1 x ( 2 m + 1 ) e − i ω k 2 m = ∑ m = 0 N 2 − 1 x ( 2 m ) e − i ω k 2 m − e − i ω k ∑ m = 0 N 2 − 1 x ( 2 m + 1 ) e − i ω k 2 m \begin{aligned}
X\Big(k + \frac{N}{2}\Big) &= \sum_{m = 0}^{\frac{N}{2} - 1}x(2m)e^{-i\omega_{k} 2m} +
e^{-i\omega_{k + \frac{N}{2}}} \sum_{m = 0}^{\frac{N}{2} - 1}x(2m + 1)e^{-i\omega_{k} 2m}\\
&= \sum_{m = 0}^{\frac{N}{2} - 1}x(2m)e^{-i\omega_{k} 2m} -
e^{-i\omega_{k}} \sum_{m = 0}^{\frac{N}{2} - 1}x(2m + 1)e^{-i\omega_{k} 2m}
\end{aligned} X ( k + 2 N ) = m = 0 ∑ 2 N − 1 x ( 2 m ) e − i ω k 2 m + e − i ω k + 2 N m = 0 ∑ 2 N − 1 x ( 2 m + 1 ) e − i ω k 2 m = m = 0 ∑ 2 N − 1 x ( 2 m ) e − i ω k 2 m − e − i ω k m = 0 ∑ 2 N − 1 x ( 2 m + 1 ) e − i ω k 2 m In summary:
X ( k ) = ∑ m = 0 N 2 − 1 x ( 2 m ) e − i ω k 2 m + e − i ω k ∑ m = 0 N 2 − 1 x ( 2 m + 1 ) e − i ω k 2 m X ( k + N 2 ) = ∑ m = 0 N 2 − 1 x ( 2 m ) e − i ω k 2 m − e − i ω k ∑ m = 0 N 2 − 1 x ( 2 m + 1 ) e − i ω k 2 m \begin{aligned}
X(k) &= \sum_{m = 0}^{\frac{N}{2} - 1}x(2m)e^{-i\omega_{k} 2m} +
e^{-i\omega_{k}} \sum_{m = 0}^{\frac{N}{2} - 1}x(2m + 1)e^{-i\omega_{k} 2m}\\
X\Big(k + \frac{N}{2}\Big) &= \sum_{m = 0}^{\frac{N}{2} - 1}x(2m)e^{-i\omega_{k} 2m} -
e^{-i\omega_{k}} \sum_{m = 0}^{\frac{N}{2} - 1}x(2m + 1)e^{-i\omega_{k} 2m}
\end{aligned} X ( k ) X ( k + 2 N ) = m = 0 ∑ 2 N − 1 x ( 2 m ) e − i ω k 2 m + e − i ω k m = 0 ∑ 2 N − 1 x ( 2 m + 1 ) e − i ω k 2 m = m = 0 ∑ 2 N − 1 x ( 2 m ) e − i ω k 2 m − e − i ω k m = 0 ∑ 2 N − 1 x ( 2 m + 1 ) e − i ω k 2 m We can simplify the above code using the above equation as follows:
fft <- function ( x ) {
n <- length ( x )
if ( n == 1 ) {
return ( x )
}
even <- fft2 ( x [ seq ( 1 , n , 2 )])
odd <- fft2 ( x [ seq ( 2 , n , 2 )])
factor <- exp ( -1i * 2 * pi / n * seq ( 0 , n / 2 - 1 ))
result <- c ( even + factor * odd , even - factor * odd )
result
}
> fft ( xn )
[ 1 ] 1.000000e+00+0.000000e+00i 3.843377e+00+1.591979e+00i
[ 3 ] -3.019377e-01-3.019377e-01i -4.143904e-02-1.000427e-01i
[ 5 ] 7.771561e-16+0.000000e+00i -4.143904e-02+1.000427e-01i
[ 7 ] -3.019377e-01+3.019377e-01i 3.843377e+00-1.591979e+00i Note that the above algorithm assumes the number of elements is a power of 2. More work (like padding or other ways) are needed to support non-power of 2.
To recap:
x n = 1 N ∑ k = 0 N − 1 X k e + i ω k n ω k = 2 π k N n = { 0 , ⋯ , N − 1 } \begin{aligned}
x_{n} &= \frac{1}{N}\sum_{k = 0}^{N-1}X_{k}e^{+i \omega_{k}n}\\[5pt]
\omega_{k} &= 2\pi \frac{k}{N}\\[5pt]
n &= \{0, \cdots, N-1\}
\end{aligned} x n ω k n = N 1 k = 0 ∑ N − 1 X k e + i ω k n = 2 π N k = { 0 , ⋯ , N − 1 } Using the same example as before:
> fft ( xn )
[ 1 ] 0.9999999999999990007993+0.0000000i 3.8433767734721757669547+1.5919788i
[ 3 ] -0.3019377358048381254640-0.3019377i -0.0414390376673379190464-0.1000427i
[ 5 ] 0.0000000000000007771561+0.0000000i -0.0414390376673376970018+0.1000427i
[ 7 ] -0.3019377358048382364863+0.3019377i 3.8433767734721762110439-1.5919788i Firstly, let us do the long calculation that has O ( n 2 ) O(n^{2}) O ( n 2 )
fs <- 8
fn <- fs / 2
f <- 1
omega <- 2 * pi * f
n <- seq ( 0 , 1 , length = fs )
xn <- cos ( omega * n [ 1 : length ( n )])
freqs <- fft ( xn )
N <- length ( xn )
xn2 <- rep ( 0 , N )
for ( n in 1 : N ) {
for ( k in 1 : N ) {
xn2 [ n ] <- xn2 [ n ] + freqs [ k ] * exp ( 1i * 2 * pi * ( k - 1 ) / N * ( n - 1 ))
}
}
> Re ( xn2 / N )
[ 1 ] 1.0000000 0.6234898 -0.2225209 -0.9009689 -0.9009689 -0.2225209 0.6234898 1.0000000
> xn
[ 1 ] 1.0000000 0.6234898 -0.2225209 -0.9009689 -0.9009689 -0.2225209 0.6234898 1.0000000 We can see that the result match with each other. Next we employ inverse FFT. FFT and IFFT have similar form, so the computation for FFT can be easily generalized to IFFT:
gsignal :: ifft ( freqs )
[ 1 ] 1.0000000 0.6234898 -0.2225209 -0.9009689 -0.9009689 -0.2225209 0.6234898 1.0000000 Let us see the representation from both fourier series and fourier transform on x ( t ) x(t) x ( t )
x ( t ) = ∑ k = − ∞ ∞ C k e i k ω 0 t x ( t ) = 1 2 π ∫ − ∞ ∞ X ( ω ) e i k ω t d ω x n = 1 N ∑ k = 0 N − 1 X k e + i ω k n \begin{aligned}
x(t) &= \sum_{k = -\infty}^{\infty}C_{k}e^{ik\omega_{0}t}\\
x(t) &= \frac{1}{2\pi}\int_{-\infty}^{\infty}X(\omega)e^{ik\omega t}\ d\omega\\
x_{n} &= \frac{1}{N}\sum_{k = 0}^{N-1}X_{k}e^{+i \omega_{k}n}
\end{aligned} x ( t ) x ( t ) x n = k = − ∞ ∑ ∞ C k e i k ω 0 t = 2 π 1 ∫ − ∞ ∞ X ( ω ) e i k ω t d ω = N 1 k = 0 ∑ N − 1 X k e + i ω k n The main reason why we do a Fourier tranform rather than FSC is because Fourier Transform can be used for both aperiodic and periodic signals.
Fourier Series
The Intuitive Guide to Fourier Analysis and Spectral Estimation: with Matlab
Iain Explains Signals, Systems, and Digital Comms (YouTube)
Python Numerical Methods
An Introduction and Analysis of the Fast Fourier Transform (Thao Nguyen)
