Fourier Transform¶

Like we have seen earlier, the fourier transformation is a central element in signal processing. It is our best tool to analyse signal in our two favourite dimensions : time and frequency. Let's dig a little to understand this tool 😎

Let's begin by understanding how a complexe signal is made from many pure ton (simple periodic signal mix together). First let's import a sound like in our prevous article and plot his waveform.

In [1]:

import librosa
import librosa.display
import scipy as sp
import IPython.display as ipd
import matplotlib.pyplot as plt
import numpy as np

import librosa import librosa.display import scipy as sp import IPython.display as ipd import matplotlib.pyplot as plt import numpy as np

In [2]:

#load audio file in the player
audio_path = "./files/piano_c5.wav"
ipd.Audio(audio_path)

#load audio file in the player audio_path = "./files/piano_c5.wav" ipd.Audio(audio_path)

Out[2]:

In [3]:

# load audio file
signal, sr = librosa.load(audio_path)

# load audio file signal, sr = librosa.load(audio_path)

In [4]:

# plot waveform
plt.figure(figsize=(18, 8))
librosa.display.waveshow(signal, sr=sr, alpha=0.5)
plt.show()

# plot waveform plt.figure(figsize=(18, 8)) librosa.display.waveshow(signal, sr=sr, alpha=0.5) plt.show()

Sinusoide superposition¶

As we said earlier, complexe signals are simple the sum of simple signal, let's construct a sum signal with three other signals $p_{1,2,3}$ where :

$$\begin{cases} p_1=sin(2\pi ft)\\ p_2=sin(2\pi 2ft)\\ p_3=sin(2\pi 3ft)\\ \end{cases}$$

Let's see the shape of the signal sum 😎

In [17]:

#define frequency and time 
f = 1
t = np.linspace(0, 10, 10000)

sin = np.sin(2*np.pi * (f * t))
sin2 = np.sin(2*np.pi * (2*f * t))
sin3 = np.sin(2*np.pi * (3*f * t))

#sum signals
sum_signal = sin + sin2 + sin3

#plot all the 4 graphs 
plt.figure(figsize=(15, 10))
plt.subplot(4, 1, 1)
plt.plot(t, sum_signal, color="r")
plt.title('Sum of signals (1+2+3)')

plt.subplot(4, 1, 2)
plt.plot(t, sin)
plt.title('Signal 1')

plt.subplot(4, 1, 3)
plt.plot(t, sin2)
plt.title('Signal 2')

plt.subplot(4, 1, 4)
plt.plot(t, sin3)
plt.title('Signal 3')
plt.show()

#define frequency and time f = 1 t = np.linspace(0, 10, 10000) sin = np.sin(2*np.pi * (f * t)) sin2 = np.sin(2*np.pi * (2*f * t)) sin3 = np.sin(2*np.pi * (3*f * t)) #sum signals sum_signal = sin + sin2 + sin3 #plot all the 4 graphs plt.figure(figsize=(15, 10)) plt.subplot(4, 1, 1) plt.plot(t, sum_signal, color="r") plt.title('Sum of signals (1+2+3)') plt.subplot(4, 1, 2) plt.plot(t, sin) plt.title('Signal 1') plt.subplot(4, 1, 3) plt.plot(t, sin2) plt.title('Signal 2') plt.subplot(4, 1, 4) plt.plot(t, sin3) plt.title('Signal 3') plt.show()

Fourier transform using scipy¶

Let's use the fft.fft() scipy function in order to plot the FT of our signal :

In [5]:

# derive spectrum using FT
ft = sp.fft.fft(signal)
magnitude = np.absolute(ft)
frequency = np.linspace(0, sr, len(magnitude))

# derive spectrum using FT ft = sp.fft.fft(signal) magnitude = np.absolute(ft) frequency = np.linspace(0, sr, len(magnitude))

In [6]:

# plot spectrum
plt.figure(figsize=(18, 8))
plt.plot(frequency[:5000], magnitude[:5000]) # magnitude spectrum
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")
plt.show()

# plot spectrum plt.figure(figsize=(18, 8)) plt.plot(frequency[:5000], magnitude[:5000]) # magnitude spectrum plt.xlabel("Frequency (Hz)") plt.ylabel("Magnitude") plt.show()

Zoom in our audio signal¶

If you zoom on the piano note you will see the sinosoid appears like magic 🪄

In [8]:

# zomm in to the waveform
samples = range(len(signal))
t = librosa.samples_to_time(samples, sr=sr)

plt.figure(figsize=(18, 8))
plt.plot(t[10000:10400], signal[10000:10400]) 
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.show()

# zomm in to the waveform samples = range(len(signal)) t = librosa.samples_to_time(samples, sr=sr) plt.figure(figsize=(18, 8)) plt.plot(t[10000:10400], signal[10000:10400]) plt.xlabel("Time (s)") plt.ylabel("Amplitude") plt.show()

Now let's return to our simple signal above (our 3 signals) and let's see how their fourier transformation looks like 🤓

In [20]:

# Function to plot the time domain and frequency domain representation of signals
def plot_signals_and_fourier_transforms():
    # Define frequency and time
    f = 1
    t = np.linspace(0, 10, 10000)
    
    # Define signals
    sin = np.sin(2*np.pi * (f * t))
    sin2 = np.sin(2*np.pi * (2*f * t))
    sin3 = np.sin(2*np.pi * (3*f * t))
    
    # Sum of signals
    sum_signal = sin + sin2 + sin3
    
    # Plotting time domain representation of signals
    plt.figure(figsize=(15, 20))
    
    signals = [sum_signal, sin, sin2, sin3]
    titles = ['Sum of signals (1+2+3)', 'Signal 1', 'Signal 2', 'Signal 3']
    
    # Plotting frequency domain representation of signals using Fourier Transform
    for i, signal in enumerate(signals):
        # Fourier Transform
        fft = np.fft.fft(signal)
        freq = np.fft.fftfreq(t.shape[-1], d=(t[1]-t[0]))
        
        plt.subplot(8, 1, i+5)
        plt.plot(freq, np.abs(fft))
        plt.title(f'Fourier Transform of {titles[i]}')
        plt.xlim(0, 5)  # Limit frequency axis for better visualization
    
    plt.tight_layout()
    plt.show()

plot_signals_and_fourier_transforms()

# Function to plot the time domain and frequency domain representation of signals def plot_signals_and_fourier_transforms(): # Define frequency and time f = 1 t = np.linspace(0, 10, 10000) # Define signals sin = np.sin(2*np.pi * (f * t)) sin2 = np.sin(2*np.pi * (2*f * t)) sin3 = np.sin(2*np.pi * (3*f * t)) # Sum of signals sum_signal = sin + sin2 + sin3 # Plotting time domain representation of signals plt.figure(figsize=(15, 20)) signals = [sum_signal, sin, sin2, sin3] titles = ['Sum of signals (1+2+3)', 'Signal 1', 'Signal 2', 'Signal 3'] # Plotting frequency domain representation of signals using Fourier Transform for i, signal in enumerate(signals): # Fourier Transform fft = np.fft.fft(signal) freq = np.fft.fftfreq(t.shape[-1], d=(t[1]-t[0])) plt.subplot(8, 1, i+5) plt.plot(freq, np.abs(fft)) plt.title(f'Fourier Transform of {titles[i]}') plt.xlim(0, 5) # Limit frequency axis for better visualization plt.tight_layout() plt.show() plot_signals_and_fourier_transforms()

Complex number in Fourier Transform¶

Complex numbers, consisting of a real part and an imaginary part, represent both magnitude and phase. In signal processing, this is crucial for understanding not just the amplitude of a frequency component (magnitude) but also its timing or phase relationship with other components.

The Fourier transform decomposes a signal into its constituent frequencies, and using complex numbers allows each component's phase information to be retained alongside its amplitude.

Euler's formula states that $e^{iθ}=cos(θ)+isin(θ)$, where $i$ is the imaginary unit. This formula is the cornerstone of connecting complex exponential functions with trigonometric functions, which are fundamental in describing oscillatory phenomena like waves.

Complex numbers are not just a mathematical convenience. They are a fundamental aspect of how we understand, analyze, and manipulate signals. Their use in Fourier transforms and signal processing allows for a more nuanced and comprehensive analysis of signals, capturing both their amplitude and phase information, which is crucial for a wide range of audio processing applications.

This animation is from the excellent youtube video of 3Blue1Brown we will try to implement one of the fancy transformation you've seen on this video. Let's begin by writing some helpers functions create_signal and plot_signal.

In [21]:

import cmath

import cmath

In [30]:

def create_signal(frequency, time):
    sin = np.sin(2 * np.pi * (frequency * time))
    sin2 = np.sin(2 * np.pi * (2 * frequency * time))
    sin3 = np.sin(2 * np.pi * (3 * frequency * time))
    return sin + sin2 + sin3

def plot_signal(signal, time):
    plt.figure(figsize=(16, 9))
    plt.plot(signal, time)
    plt.xlabel("Time")
    plt.ylabel("Intensity")
    plt.show()

def create_signal(frequency, time): sin = np.sin(2 * np.pi * (frequency * time)) sin2 = np.sin(2 * np.pi * (2 * frequency * time)) sin3 = np.sin(2 * np.pi * (3 * frequency * time)) return sin + sin2 + sin3 def plot_signal(signal, time): plt.figure(figsize=(16, 9)) plt.plot(signal, time) plt.xlabel("Time") plt.ylabel("Intensity") plt.show()

In [23]:

def calculate_centre_of_gravity(mult_signal):
    x_centre = np.mean([x.real for x in mult_signal])
    y_centre = np.mean([x.imag for x in mult_signal])
    return x_centre, y_centre

def calculate_centre_of_gravity(mult_signal): x_centre = np.mean([x.real for x in mult_signal]) y_centre = np.mean([x.imag for x in mult_signal]) return x_centre, y_centre

In [24]:

def calculate_sum(mult_signal):
    x_sum = np.sum([x.real for x in mult_signal])
    y_sum = np.sum([x.imag for x in mult_signal])
    return x_sum, y_sum

def calculate_sum(mult_signal): x_sum = np.sum([x.real for x in mult_signal]) y_sum = np.sum([x.imag for x in mult_signal]) return x_sum, y_sum

In [25]:

def create_pure_tone(frequency, time):
    angle = -2 * np.pi * frequency * time
    return np.cos(angle) + 1j * np.sin(angle)

def create_pure_tone(frequency, time): angle = -2 * np.pi * frequency * time return np.cos(angle) + 1j * np.sin(angle)

In [26]:

def plot_fourier_transform(pure_tone_frequency, 
                           signal_frequency, 
                           time, 
                           plot_centre_of_gravity=False,
                           plot_sum=False):
    
    # create sinusoid and signal
    pure_tone = create_pure_tone(pure_tone_frequency, time)
    signal = create_signal(signal_frequency, time)
    
    # multiply pure tone and signal
    mult_signal = pure_tone * signal
    
    X = [x.real for x in mult_signal]
    Y = [x.imag for x in mult_signal]

    plt.figure(figsize=(15, 10))
    plt.plot(X, Y, 'o')

    # calculate and plot centre of gravity
    if plot_centre_of_gravity:
        centre_of_gravity = calculate_centre_of_gravity(mult_signal)
        plt.plot([centre_of_gravity[0]], [centre_of_gravity[1]], marker='o', markersize=10, color="red")


    # calculate and plot sum 
    if plot_sum:
        integral = calculate_sum(mult_signal)
        plt.plot([integral[0]], [integral[1]], marker='o', markersize=10, color="green")

    
    # set origin axes
    ax = plt.gca()
    ax.grid(True)
    ax.spines['left'].set_position('zero')
    ax.spines['right'].set_color('none')
    ax.spines['bottom'].set_position('zero')
    ax.spines['top'].set_color('none')

    if not plot_sum:
        plt.xlim(-3, 3)
        plt.ylim(-3, 3)

    plt.show()

def plot_fourier_transform(pure_tone_frequency, signal_frequency, time, plot_centre_of_gravity=False, plot_sum=False): # create sinusoid and signal pure_tone = create_pure_tone(pure_tone_frequency, time) signal = create_signal(signal_frequency, time) # multiply pure tone and signal mult_signal = pure_tone * signal X = [x.real for x in mult_signal] Y = [x.imag for x in mult_signal] plt.figure(figsize=(15, 10)) plt.plot(X, Y, 'o') # calculate and plot centre of gravity if plot_centre_of_gravity: centre_of_gravity = calculate_centre_of_gravity(mult_signal) plt.plot([centre_of_gravity[0]], [centre_of_gravity[1]], marker='o', markersize=10, color="red") # calculate and plot sum if plot_sum: integral = calculate_sum(mult_signal) plt.plot([integral[0]], [integral[1]], marker='o', markersize=10, color="green") # set origin axes ax = plt.gca() ax.grid(True) ax.spines['left'].set_position('zero') ax.spines['right'].set_color('none') ax.spines['bottom'].set_position('zero') ax.spines['top'].set_color('none') if not plot_sum: plt.xlim(-3, 3) plt.ylim(-3, 3) plt.show()

In [31]:

time = np.linspace(0, 10, 10000)
signal = create_signal(frequency=1, time=time)
plot_signal(time, signal)

time = np.linspace(0, 10, 10000) signal = create_signal(frequency=1, time=time) plot_signal(time, signal)

In [32]:

time = np.linspace(0, 1, 10000)
plot_fourier_transform(pure_tone_frequency=1.1,
                       signal_frequency=1,
                       time=time,
                       plot_centre_of_gravity=False,
                       plot_sum=False)

time = np.linspace(0, 1, 10000) plot_fourier_transform(pure_tone_frequency=1.1, signal_frequency=1, time=time, plot_centre_of_gravity=False, plot_sum=False)