FFT in Python

Generating Sinusoids¶

In [1]:

import numpy as np
import matplotlib.pyplot as plt
# Generate a few sinusoids

samplerate = 1000  # 1000 points every second
# 1 second of audio
timeAxis  = np.linspace(0, 1, samplerate)
freq1 = 5
freq2 = 10
freq3 = 25

sineWave1 = 0.5*np.sin(2*np.pi*freq1*timeAxis)
sineWave2 = 0.5*np.sin(2*np.pi*freq2*timeAxis)
sineWave3 = 0.5*np.sin(2*np.pi*freq3*timeAxis)

# Sum of SineWaves
sineWaveSum = sineWave1 + sineWave2 + sineWave3
sineWaveSum = sineWaveSum/np.max(sineWaveSum)

In [2]:

# Create subplots
fig, axes = plt.subplots(4, 1, figsize=(10, 8))  # 4 rows, 1 column

# Plot each sine wave in its respective subplot
axes[0].plot(timeAxis, sineWave1, color='b')
axes[0].set_title(f'Sine Wave 1: {freq1} Hz')
axes[0].set_xlabel('t/s')
axes[0].set_ylabel('Amplitude = x(t)')

axes[1].plot(timeAxis, sineWave2, color='g')
axes[1].set_title(f'Sine Wave 2: {freq2} Hz')
axes[1].set_xlabel('t/s')
axes[1].set_ylabel('Amplitude = x(t)')

axes[2].plot(timeAxis, sineWave3, color='r')
axes[2].set_title(f'Sine Wave 3: {freq3} Hz')
axes[2].set_xlabel('t/s')
axes[2].set_ylabel('Amplitude = x(t)')

axes[3].plot(timeAxis, sineWaveSum, color='k')
axes[3].set_title(f'Sine Wave Sum')
axes[3].set_xlabel('t/s')
axes[3].set_ylabel('Amplitude = x(t)')


# Clean
plt.tight_layout()
plt.show()

No description has been provided for this image

Implementing and Visualizing the FFT¶

FFT (Analysis): $$ X[k] = \sum_{n=0}^{N-1} x[n]e^{-j2\pi kn/N} $$

where,
$$e^{-j2\pi kn/N} = \cos(\frac{2\pi kn}{N}) - j\sin(\frac{2\pi kn}{N})$$

Length of N ~ Length of DFT

https://ringbuffer.org/dsp/Fourier_Transform/discrete-fourier-transform/

Also checkout the DFT lectures from the Audio Signal Processing Coursera course: https://www.coursera.org/learn/audio-signal-processing/lecture/EZRXC/dft-1

In [3]:

from scipy.fft import fft, fftfreq, ifft

#Generate the sinusoid

samplerate = 1000  
fs = samplerate

# 1 second of audio
timeAxis  = np.linspace(0, 1, samplerate)
freq1 = 100

sineWave1 = 0.5*np.sin(2*np.pi*freq1*timeAxis)

# Sum of SineWaves
sineWaveSum = sineWave1 
sineWaveSum = sineWaveSum/np.max(sineWaveSum)

Setup for FFT¶

In [4]:

# Number of samples in the input signal
nX = len(sineWaveSum)
# Find the next power of 2 for zero-padding
nxpadded = 2 ** int(np.ceil(np.log2(nX)))

Compute fft with scipy.fft.fft function¶

In [5]:

# Compute FFT
X = fft(sineWaveSum, nxpadded)

Discarding half of the result because the DFT of real signals is symmetric around 0¶

In [6]:

# Compute Magnitude Spectrum
Xabs = np.abs(X[:nxpadded // 2])
# Compute phase spectrum
Xphase = np.angle(X[:nxpadded // 2])
# Unwrap the phase
unwrapphase = np.unwrap(Xphase)

Create an array for the frequency bins evenly spaced (using the fftfreq in scipy)¶

In [7]:

frequencies = fftfreq(nxpadded, 1/fs)[:nxpadded // 2]
#Another way -> using formula for bin spacing df = fs/N
frequencies_manual = np.linspace(0, fs/2, nxpadded//2, endpoint=True)

Plotting¶

In [8]:

# Create subplots
fig, axes = plt.subplots(3, 1, figsize=(10, 8))  # 3 rows, 1 column


# Plot each sine wave in its respective subplot
axes[0].plot(frequencies,Xabs)
axes[0].set_title("Magnitude Spectrum")
axes[0].set_xlabel("f/Hz")
axes[0].set_ylabel("|H(f)|")

# Plot each sine wave in its respective subplot
axes[1].plot(frequencies, Xphase)
axes[1].set_title("Phase Spectrum")
axes[1].set_xlabel("f/Hz")
axes[1].set_ylabel("Phase (radians)")

# Plot each sine wave in its respective subplot
axes[2].plot(frequencies, unwrapphase)
axes[2].set_title("Unwrapped Phase")
axes[2].set_xlabel("f/Hz")
axes[2].set_ylabel("Phase (radians)")


plt.tight_layout()
plt.show()

IFFT from Magnitude and Phase¶

FFT (Analysis): $$ X[k] = \sum_{n=0}^{N-1} x[n]e^{-j2\pi kn/N} $$

Inverse FFT (Synthesis) with proper scaling: $$ x[n] = \frac{1}{N}\sum_{k=0}^{N-1} X[k]e^{j2\pi kn/N} $$

Energy conservation (Parseval's theorem): $$\sum_{n=0}^{N-1} |x[n]|^2 = \frac{1}{N}\sum_{k=0}^{N-1} |X[k]|^2$$

Note the $\frac{1}{N}$ scaling factor in IFFT - this is crucial for maintaining energy between time and frequency domains.

Generating the Signal¶

In [9]:

from scipy.fft import fft, fftfreq, ifft

#Generate the sinusoid

samplerate = 1000  

# 1 second of audio
timeAxis  = np.linspace(0, 1, samplerate)
freq1 = 100

sineWave1 = 0.5*np.sin(2*np.pi*freq1*timeAxis)

# Sum of SineWaves
sineWaveSum = sineWave1 
#sineWaveSum = sineWaveSum/np.max(sineWaveSum)

Compute the FFT¶

In [10]:

# Number of samples in the input signal
nX = len(sineWaveSum)
# Find the next power of 2 for zero-padding
nxpadded = 2 ** int(np.ceil(np.log2(nX)))

# FFT
X = fft(sineWaveSum, nxpadded)

# Magnitude and phase calculations
Xabs = np.abs(X[:nxpadded // 2])
Xphase = np.angle(X[:nxpadded // 2])
unwrapphase = np.unwrap(Xphase)

Use the symmetry of DFT and reconstruct X from Xabs and Xphase¶

In [11]:

# Reconstruct
halfspectrum = Xabs * np.exp(1j * unwrapphase) # Eulers formula

# Setup the axis for the full spectrum
fullspectrum = np.zeros(nxpadded, dtype=complex)

fullspectrum[0] = X[0]  # DC for the 0th bin position

# Hermitian symmetry
fullspectrum[1:nxpadded//2] = halfspectrum[1:]
fullspectrum[nxpadded//2] = X[nxpadded//2]  # Nyquist
# Same half but go backwards and take the complex conjugate
fullspectrum[nxpadded//2+1:] = np.conj(halfspectrum[1:][::-1])

In [12]:

# Reconstruct
halfspectrum = Xabs * np.exp(1j * unwrapphase)
fullspectrum = np.zeros(nxpadded, dtype=complex)
fullspectrum[0] = X[0]  # DC
fullspectrum[1:nxpadded//2] = halfspectrum[1:]
fullspectrum[nxpadded//2] = X[nxpadded//2]  # Nyquist
fullspectrum[nxpadded//2+1:] = np.conj(halfspectrum[1:][::-1])

# IFFT with 1/N scaling, ifft already handles for it
Signal = ifft(fullspectrum, nxpadded)
Signal = Signal.real[:samplerate]
plt.plot(Signal)

Out[12]:

[<matplotlib.lines.Line2D at 0x782af0070ce0>]

In [13]:

# Create subplots
fig, axes = plt.subplots(2, 1, figsize=(10, 8))  # 4 rows, 1 column

# Plot each sine wave in its respective subplot
axes[0].plot(timeAxis, sineWaveSum, color='b')
axes[0].set_title(f'Sine Wave : {100} Hz')
axes[0].set_xlabel('t/s')
axes[0].set_ylabel('Amplitude  = x(t)')

axes[1].plot(timeAxis, Signal, color='g')
axes[1].set_title('Reconstructed WaveForm')
axes[1].set_xlabel('t/s')
axes[1].set_ylabel('Amplitude = x(t)')

# Clean
plt.tight_layout()
plt.show()

FFT of common sequences¶

Impulse

In [14]:

samplerate = 1000  

# 1 second
timeAxis  = np.linspace(0, 1, samplerate)
imp = np.zeros(samplerate)
imp[0] = 1

#Compute FFT

# Number of samples in the input signal
nX = len(imp)
# Find the next power of 2 for zero-padding
nxpadded = 2 ** int(np.ceil(np.log2(nX)))

# FFT
X = fft(imp, nxpadded)

# Magnitude and phase calculations
Xabs = np.abs(X[:nxpadded // 2])
Xphase = np.angle(X[:nxpadded // 2])
unwrapphase = np.unwrap(Xphase)
frequencies = fftfreq(nxpadded, 1/fs)[:nxpadded // 2]

In [15]:

# Create subplots
fig, axes = plt.subplots(4, 1, figsize=(10, 8))  # 3 rows, 1 column

# Plot each sine wave in its respective subplot
axes[0].plot(timeAxis,imp)
axes[0].set_title("Impulse")
axes[0].set_xlabel("n")
axes[0].set_ylabel("Amplitude = x[n]")


# Plot each sine wave in its respective subplot
axes[1].plot(frequencies,Xabs)
axes[1].set_title("Magnitude Spectrum")
axes[1].set_xlabel("f/Hz")
axes[1].set_ylabel("|H(f)|")

# Plot each sine wave in its respective subplot
axes[2].plot(frequencies, Xphase)
axes[2].set_title("Phase Spectrum")
axes[2].set_xlabel("f/Hz")
axes[2].set_ylabel("Phase (radians)")

# Plot each sine wave in its respective subplot
axes[3].plot(frequencies, unwrapphase)
axes[3].set_title("Unwrapped Phase")
axes[3].set_xlabel("f/Hz")
axes[3].set_ylabel("Phase (radians)")


plt.tight_layout()
plt.show()

BoxCar

In [16]:

samplerate = 1000  

# 1 second
timeAxis  = np.linspace(0, 1, samplerate)

width = 200  # Width of boxcar
boxcar = np.zeros(samplerate)
boxcar[200: 200+width] = 1  # Set middle portion to 1


#Compute FFT

# Number of samples in the input signal
nX = len(boxcar)
# Find the next power of 2 for zero-padding
nxpadded = 2 ** int(np.ceil(np.log2(nX)))

# FFT
X = fft(boxcar, nxpadded)

# Magnitude and phase calculations
Xabs = np.abs(X[:nxpadded // 2])
Xphase = np.angle(X[:nxpadded // 2])
unwrapphase = np.unwrap(Xphase)
frequencies = fftfreq(nxpadded, 1/fs)[:nxpadded // 2]

In [17]:

# Create subplots
fig, axes = plt.subplots(4, 1, figsize=(10, 8))  # 3 rows, 1 column

# Plot each sine wave in its respective subplot
axes[0].plot(timeAxis,boxcar)
axes[0].set_title("Boxcar")
axes[0].set_xlabel("n")
axes[0].set_ylabel("Amplitude = x[n]")


# Plot each sine wave in its respective subplot
axes[1].plot(frequencies,Xabs)
axes[1].set_title("Magnitude Spectrum")
axes[1].set_xlabel("f/Hz")
axes[1].set_ylabel("|H(f)|")

# Plot each sine wave in its respective subplot
axes[2].plot(frequencies, Xphase)
axes[2].set_title("Phase Spectrum")
axes[2].set_xlabel("f/Hz")
axes[2].set_ylabel("Phase (radians)")

# Plot each sine wave in its respective subplot
axes[3].plot(frequencies, unwrapphase)
axes[3].set_title("Unwrapped Phase")
axes[3].set_xlabel("f/Hz")
axes[3].set_ylabel("Phase (radians)")


plt.tight_layout()
plt.show()