[DSP] Section 4: 푸리에 변환

[DSP] Section 4

[REF | AudioSignalProcessingForML]

10. Fourier Transform: The Intuition

Intuition

• Decompose a complex sound into its frequency components
• Compare signal with sinusoids of various frequencies
• For each freq. we get a magnitude and a phase
• High magnitude indicates high similarity between the signal and a sinusoid

From time to frequency domain

• Time -> TF -> Freq.

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

audio_path = "./raw/10_audio/piano_c.wav"
ipd.Audio(audio_path)

signal, sr = librosa.load(audio_path)

# # plot waveform
plt.figure(figsize=(10, 5))
librosa.display.waveplot(signal, sr=sr, alpha=0.5)
plt.show()

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

# plot spectrum
plt.figure(figsize=(10, 5))
plt.plot(frequency[:5000], magnitude[:5000])
plt.xlabel("Freq. (Hz)")
plt.ylabel("Magnitude")
plt.show()

len(signal)

33968

d = 1/sr
d

4.5351473922902495e-05

d_523 = 1/523
d_523

0.0019120458891013384

d_400_samples = 400*d
d_400_samples

0.018140589569160998

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

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

Sine Wave

\(sin(2\pi \times (ft-\varphi))\)

# create sinusoids
f = 523
phase = 0
phase2 = 0.2
phase3 = 0.5

sin = 0.5 * np.sin(2*np.pi * (f*t-phase))
sin2 = 0.5 * np.sin(2*np.pi * (f*t-phase2))
sin3 = 0.5 * np.sin(2*np.pi * (f*t-phase3))

plt.figure(figsize=(10, 5))
plt.plot(t[10000:10400], sin[10000:10400], color="r")
plt.plot(t[10000:10400], sin2[10000:10400], color="y")
plt.plot(t[10000:10400], sin3[10000:10400], color="b")

plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.show()

# compare signal and sinusoids
f = 523
phase = 0
phase2 = 0.55

sin = 0.1*np.sin(2*np.pi*(f*t-phase))
sin2 = 0.1*np.sin(2*np.pi*(f*t-phase2))

plt.figure(figsize=(10, 5))
plt.plot(t[10000:10400], signal[10000:10400])
plt.plot(t[10000:10400], sin[10000:10400], color="r")
plt.plot(t[10000:10400], sin2[10000:10400], color="y")

plt.fill_between(t[10000:10400], sin2[10000:10400]*signal[10000:10400])

plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.show()

Fourier transform: Step by step

• Choose a frequency -> Optimise phase -> Calculate magnitude

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

Fourier transform

\(\varphi_f = argmax_{ \varphi \in [0, 1) } (\int s(t) \cdot \sin(2\pi \cdot (ft-\varphi) ) \cdot dt)\)

\[d_f = max_{\varphi \in [0, 1)} (\int s(t) \cdot sin(2\pi \cdot (ft-\varphi)) \cdot dt)\]

Reconstructing a signal

• Superimpose sinusoids
• Weight them by the relative magnitude
• Use relative phase
• Original signal and FT have same information

Inverse Fourier Transform

11. Complex numbers for audio signal processing

Why bother with complex numbers?

• Fourier transform -> magnitude and phase
• Magnitude is a real number
• … someting with magnitude + phase? -> COMPLEX Number

The genesis of CNs (복소수의 기원)

• $\sqrt -1 = i^2 = -1$
• $c=a+ib$

Polar coordinate representation

• $cos(\gamma) = { a \over

  c

  }$

• $sin(\gamma) = { b \over

  c

  }$

• ${ sin(\gamma) \over cos(\gamma) } = { b \over a }$
• $\gamma = arctan( { b \over a } )$

$a=	c	\cdot cos(\gamma) \ \ \ b=	c	\cdot sin(\gamma)$
$c=a+ib =	c	\cdot ( cos(\gamma) + isin(\gamma) )$

Euler formula

$e^{i\gamma} = cos(\gamma) + isin(\gamma)$

Euler identity

$e^{i\pi} + 1 = 0$

Polar coordinates

$c=|c| \cdot e^{i\gamma}$

• $

  c

  $: Scales distance from origin

• $e^{i\gamma}$: Direction of a number in the complex plane

12. Defining the Fourier transform using complex numbers

The intuition

• Use magnitude and phase as polar coordinates
• Encode both coefficients in a single complex number

Complex Fourier transform coefficients

\(c_f = { d_f \over \sqrt 2 } \cdot e^{-i2\pi\varphi_f}\)

• $d_f$: magnitude
• $\gamma = 2\pi\varphi_f$

Continuous audio signal

$g(t) \ \ \ g: R \rightarrow R $

Complex Fourier transform

$\hat g(f) = c_f$
$\hat g: R \rightarrow C$

\[\hat g(f) = \int g(t) \cdot e^{-i2\pi ft} dt\]

import cmath
import matplotlib.pyplot as plt
import numpy as np

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

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_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 create_pure_tone(freq, time):
    angle = -2*np.pi*freq*time
    return np.cos(angle) + 1j * np.sin(angle)

def plot_fourier_transform(pure_tone_freq, 
                           signal_freq, 
                           time, 
                           plot_centre_of_gravity=False,
                           plot_sum=False):
    
    # create sinusoid and signal
    pure_tone = create_pure_tone(pure_tone_freq, time)
    signal = create_signal(signal_freq, 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=(8, 5))
    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_signal(signal, time):
    plt.figure(figsize=(8, 5))
    plt.plot(signal, time)
    plt.xlabel("Time")
    plt.ylabel("Intensity")
    plt.show()

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

time = np.linspace(0, 1, 1000)
plot_fourier_transform(pure_tone_freq=1,
                       signal_freq=1,
                       time=time
                       )

time = np.linspace(0, 1, 1000)
plot_fourier_transform(pure_tone_freq=1.1,
                       signal_freq=1,
                       time=time
                       )

time = np.linspace(0, 2, 1000)
plot_fourier_transform(pure_tone_freq=1.1,
                       signal_freq=1,
                       time=time
                       )

time = np.linspace(0, 10, 1000)
plot_fourier_transform(pure_tone_freq=1.1,
                       signal_freq=1,
                       time=time
                       )

time = np.linspace(0, 10, 1000)
plot_fourier_transform(pure_tone_freq=2,
                       signal_freq=1,
                       time=time
                       )

time = np.linspace(0, 10, 1000)
plot_fourier_transform(pure_tone_freq=1,
                       signal_freq=1,
                       time=time,
                       plot_centre_of_gravity=True
                       )

time = np.linspace(0, 10, 1000)
plot_fourier_transform(pure_tone_freq=1,
                       signal_freq=1,
                       time=time,
                       plot_centre_of_gravity=True,
                       plot_sum=True
                       )

Complex Fourier transform

\(\hat g(f) = \int g(t) \cdot e^{-i\pi ft} dt = \int g(t) \cdot cos(-2\pi ft)dt + i \int g(t) \cdot sin(-2\pi ft)dt\) (Real part + Imaginart part)

Magnitude and phase

$c_f = { d_f \over \sqrt 2 } \cdot e^{-i2\pi \varphi_f}$
$d_f = \sqrt2 \cdot |\hat g(g)|$
$\varphi_f = - { \gamma f \over 2\pi}$

Inverse Fourier transform (IFT)

Forier representation

$g(t) = \int c_f \cdot e^{i2\pi ft} df$

• $e^{i2\pi ft}$: Pure tone of freq. f
• $c_f \cdot e^{i2\pi ft}$: Weight pure tone with magnitude and add phase
• $\int c_f \cdot e^{i2\pi ft} df$: Add up all (weighted) sinusoids

A Fourier roundtrip

$\hat g(f) = \int g(t) \cdot e^{-i2\pi ft} dt \ \ (FT)$
$g(t) = \int c_f \cdot e^{i2\pi ft} df \ \ (IFT)$