Fast Fourier Transform Implementation

Gino Prasad

12/17/2022

Continuous Fourier Transform

$\text{Let } f \text{ be a periodic function with period } 2 \pi$

$proj_{cos(nx)}(f) = \frac{cos \cdot f}{cos \cdot cos} = \frac{1}{\pi} \int_{-\pi}^{\pi}{f(x)(cos(nx))dx}$

$\text{Let } f \text{ be a periodic function with period } T$

$\omega = \frac{2 \pi}{T}$

$proj_{cos(n\omega x)}(f) = \frac{cos \cdot f}{cos \cdot cos} = \frac{2}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}}{f(x)(cos(n\omega x))dx} = \frac{\omega}{\pi} \int_{-\frac{T}{2}}^{\frac{T}{2}}{f(x)(cos(n\omega x))dx}$

$\text{Let } a_n = proj_{cos(n\omega x)}(f) = \frac{2}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}}{f(x)(cos(n\omega x))dx}$

$\text{Let } b_n = proj_{sin(n\omega x)}(f) = \frac{2}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}}{f(x)(sin(n\omega x))dx}$

If $n = 0$, $proj_1(f) = \frac{1 \cdot f}{1 \cdot 1} = \frac{1}{T} \int_{\frac{-T}{2}}^{\frac{T}{2}}f(x)dx = \frac{a_0}{2}$

$\frac{a_0}{2} + \sum_{n=1}^{\infty} a_n cos(n \omega x) + b_n sin(n \omega x)$

$cos(x) + (i) sin(x) = e^{ix}$

$cos(x) = \frac{e^{i x} + e^{-i x}}{2}$

$sin(x) = \frac{e^{i x} - e^{-i x}}{2 i}$

$ = \frac{a_0}{2} + \sum_{n=1}^{\infty} a_n \frac{e^{i n \omega x} + e^{-i n \omega x}}{2} + b_n \frac{e^{i n \omega x} - e^{-i n \omega x}}{2i}$

$= \frac{a_0}{2} + \sum_{n=1}^{\infty} a_n \frac{e^{i n \omega x} + e^{-i n \omega x}}{2} + b_n \frac{e^{i n \omega x} - e^{-i n \omega x}}{2i}$

$= \frac{a_0}{2} + \sum_{n=1}^{\infty} \frac{a_n - (i) b_n}{2} e^{i n \omega x} + \frac{a_n + (i) b_n}{2} e^{-i n \omega x}$

$= \frac{a_0}{2} + \sum_{n=1}^{\infty} \frac{a_n - (i) b_n}{2} e^{i n \omega x} + \sum_{n=1}^{\infty} \frac{a_n + (i) b_n}{2} e^{-i n \omega x}$

$\text{Let } c_n = \{$

$\text{n < 0: } \frac{a_{-n} + (i)b_{-n}}{2}$

$\text{n = 0: } \frac{a_0}{2}$

$\text{n < 0: } \frac{a_{n} - (i)b_{n}}{2}$

$\}$

$f(x) = \sum_{-\infty}^{\infty} c_n e^{i n \omega x}$

For $n \geq 0$:

$c_n = \frac{2}{T} \int_{\frac{-T}{2}}^{\frac{T}{2}} \frac{(cos(n \omega x) - (i)sin(n \omega x)) f(x)}{2}$

$= \frac{1}{T} \int_{\frac{-T}{2}}^{\frac{T}{2}} (cos(n \omega x) - (i) sin(n \omega x)) f(x)$

$= \frac{1}{T} \int_{\frac{-T}{2}}^{\frac{T}{2}} f(x) e^{- i n \omega x} $

For $n \leq 0$:

$c_n = \frac{2}{T} \int_{\frac{-T}{2}}^{\frac{T}{2}} \frac{(cos(n \omega x) + (i)sin(n \omega x)) f(x)}{2}$

Let n = |n|

$= \frac{1}{T} \int_{\frac{-T}{2}}^{\frac{T}{2}} (cos(n \omega x) - (i) sin(n \omega x)) f(x)$

$= \frac{1}{T} \int_{\frac{-T}{2}}^{\frac{T}{2}} f(x) e^{- i n \omega x} $

Discreet Fourier Transform

$\text{Let } \mathbf{f} \text{ be a vector in } \mathbb{R}^n$

$\text{To maintain orthogonality , the largest period (non-constant) (} \omega_0 \text{) has to be } T = \mathbf{n}$

$\text{Written as an complex exponential, this frequency is}$

$\omega_f = e^{\frac{-2 \pi i}{n}}$

$\hat{f_k} = \sum_{j=0}^{n-1}f_j e^{\frac{-2 \pi}{n} j k}$

$\hat{f_k} = \sum_{j=0}^{n-1}f_j w_f^{j k}$

$\hat{f_k} = |f_k| e^{i \phi}$

$|f_k| \text{ is the amplitude, } \phi \text{ is the phase shift}$

$\mathbf{\hat{f}} = \begin{bmatrix} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_f & \omega_f^2 & \omega_f^3 & \cdots & \omega_f^{n-1} \\ 1 & \omega_f^2 & \omega_f^4 & \omega_f^6 & \cdots & \omega_f^{2(n-1)} \\ 1 & \omega_f^3 & \omega_f^6 & \omega_f^9 & \cdots & \omega_f^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_f^{n-1} & \omega_f^{2(n-1)} & \omega_f^{3(n-1)} & \cdots & \omega_f^{(n-1)^2} \end{bmatrix} \mathbf{f}$

Fast Fourier Transform

import numpy as np
import matplotlib.pyplot as plt

def val_fft(f):
    domain_matrix = np.array([[(-2 * np.pi * j * k) / len(f) for j in range(len(f))] for k in range(len(f))])
    cos_matrix = np.cos(domain_matrix)
    sin_matrix = np.sin(domain_matrix)
    
    cos, sin = np.dot(cos_matrix, f), np.dot(sin_matrix, f)

    return cos, sin

def fft(f, exponent_sign=-1):
    f = f.astype(np.float64)
    zero = np.zeros_like(f)[:1]
    if len(f) == 1:
        return f, zero
    assert not len(f) & 1
    
    evens = f[::2]
    odds = f[1::2]
    
    fft_sum_cos, fft_sum_sin = fft(evens, exponent_sign=exponent_sign)
    fft_odds_cos, fft_odds_sin = fft(odds, exponent_sign=exponent_sign)
        
        
    series_domain = np.array([(exponent_sign * 2 * np.pi * n) / len(f) for n in range(len(f) // 2)])
    cos_series = np.cos(series_domain)
    sin_series = np.sin(series_domain)
    
    fft_sum_cos += (cos_series * fft_odds_cos) - (sin_series * fft_odds_sin) 
    fft_sum_sin += (cos_series * fft_odds_sin) + (sin_series * fft_odds_cos) 
    
    diff = np.array([(evens - odds).sum()])    
    fft_sum_cos = np.concatenate((fft_sum_cos, diff, fft_sum_cos[:0:-1]))
    fft_sum_sin = np.concatenate((fft_sum_sin, zero, -fft_sum_sin[:0:-1]))
    
    return fft_sum_cos, fft_sum_sin

def get_amplitude_spectrum(f_hat_cos, f_hat_sin):
    return np.concatenate(([0], (len(f_hat_cos) / np.arange(1, len(f_hat_cos))))), np.linalg.norm(np.stack((f_hat_cos, f_hat_sin), axis=1), axis=1)

def predict_frequencies(f, threshold=0.2, decimals=1, min_freq_threshold=0.01):
    fft_cos, fft_sin = fft(f)
    periods, amplitude_spectrum = get_amplitude_spectrum(fft_cos, fft_sin)
    
    print(f"Using Amplitude Threshold {len(f) * threshold}")
    candidates = list(set([round(freq, decimals) for freq in periods[amplitude_spectrum > (len(f) * threshold)]]))
    candidates.sort()
    return [candidate for candidate in candidates if candidate > (len(f) * min_freq_threshold)]

def simulate_wave(n, period, offset=None, noise_amplitude=0.05):
    if offset is None:
        offset = np.random.uniform(0, 2 * np.pi)
    return np.cos(((2 * np.pi / period) * np.arange(0, n)) + offset) + np.random.normal(scale=noise_amplitude, size=n)

Inverse Fast Fourier Transform

$f_n = proj_1 (\frac{1}{N} \sum_{k=0}^{N-1} \hat{f}_k e^{2 \pi i n k / N})$

def val_ifft(f_hat_cos, f_hat_sin, rec=False):
    # Encountered floating point rounding errors
    assert len(f_hat_cos) == len(f_hat_sin)
    N = len(f_hat_cos)
    reconstruction_domain_matrix = np.array([[(2 * np.pi * j * k) / N for j in range(N)] for k in range(N)])
    ifft_cos = (np.dot(np.cos(reconstruction_domain_matrix), f_hat_cos) - np.dot(np.sin(reconstruction_domain_matrix), f_hat_sin))
    ifft_sin = (np.dot(np.cos(reconstruction_domain_matrix), f_hat_sin) + np.dot(np.sin(reconstruction_domain_matrix), f_hat_cos))
    
    if not rec:
        return ifft_cos / N
    return ifft_cos, ifft_sin

def ifft(f_hat_cos, f_hat_sin):
    assert len(f_hat_cos) == len(f_hat_sin)
    real_cos, real_sin = fft(f_hat_cos, exponent_sign=1)
    imag_cos, imag_sin = fft(f_hat_sin, exponent_sign=1)
    imag_cos, imag_sin = -imag_sin, imag_cos
    return (real_cos + imag_cos) / len(f_hat_cos), (real_sin + imag_sin) / len(f_hat_sin)

n = 2 ** 12

wave_periods = [50, 70, 100]
wave_periods.sort()

waves = np.array([simulate_wave(n, period) for period in wave_periods])
waveform = np.sum(waves, axis=0)

x_axis = list(range(n))
for period, wave in zip(wave_periods, waves):
    plt.plot(x_axis, wave, label=f'wave period {period}')
plt.plot(x_axis, waveform, label='combined waveform')
plt.legend()

<matplotlib.legend.Legend at 0x7f9e50d07940>

fft_cos, fft_sin = fft(waveform)
reconstruction, reconstruction_imaginary = ifft(fft_cos, fft_sin)
periods, amplitude_spectrum = get_amplitude_spectrum(fft_cos, fft_sin)
frequencies = predict_frequencies(waveform)
print("FINISHED COMPUTING")

plt.title('Original Waveform')
plt.plot(x_axis, waveform)
plt.figure()

plt.title('Reconstruction')
plt.plot(x_axis, reconstruction)
plt.figure()

plt.title('Amplitude Spectrum')
plt.bar(x=periods, height=amplitude_spectrum, width=10)

for wave_period in wave_periods:
    plt.axvline(wave_period, color='red', linewidth=0.5, linestyle='dashed', label=f'Wave Period {wave_period}')
plt.legend()
plt.show()

print(f"\n\n\tDetected Waves with Periods: {frequencies}")
print(f"\n\n\tActual Periods: {wave_periods}")

Using Amplitude Threshold 819.2
FINISHED COMPUTING

	Detected Waves with Periods: [50.0, 69.4, 70.6, 99.9]


	Actual Periods: [50, 70, 100]