This website works better with desktop in both themes, for mobile devices please change to light theme.

# Fourier Transform#

[1]:

import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import style
from IPython import display
from matplotlib import animation

style.use('ggplot')

# %matplotlib inline


## Euler’s Number(e)#

[2]:

math.e,np.e

[2]:

(2.718281828459045, 2.718281828459045)

[3]:

plt.plot(np.exp(np.arange(-10,30,0.1)))

[3]:

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


### $$e^{\pi i}$$#

[4]:

N = 20
fig = plt.figure(figsize=(10,5))
ax1 = plt.subplot(1,2,1)
ax2 = plt.subplot(1,2,2)

ax1.set_xlim((-2, 2))
ax1.set_ylim((-2, 2))
ax1.set_xlabel("real")
ax1.set_ylabel("imaginary")

ax2.set_xlim((0,N))
ax2.set_ylim((0,N))

txt_title = ax1.set_title('data points : 0')
plot1, = ax1.plot([], [], 'b-', lw=2)
plot2, = ax2.plot([], [], 'b-', lw=2)

def draw_frame(n):
a = np.arange(n+1)
data = np.exp(1j * a)
plot1.set_data(data.real,data.imag)
plot2.set_data(range(n+1),a)
txt_title.set_text(f"data points : {n+1}")
return (plot1,)

anim = animation.FuncAnimation(fig, draw_frame, frames=N, interval=1000, blit=True)

display.HTML(anim.to_html5_video())

[4]:


### $$e^{-2\pi\frac{i}{N}}$$#

[5]:

N = 100
fig = plt.figure(figsize=(10,5))
ax1 = plt.subplot(1,2,1)
ax2 = plt.subplot(1,2,2)

ax1.set_xlim((-2, 2))
ax1.set_ylim((-2, 2))
ax1.set_xlabel("real")
ax1.set_ylabel("imaginary")

ax2.set_xlim((0,N))
ax2.set_ylim((0,N))

txt_title = ax1.set_title('data points : 0')
plot1, = ax1.plot([], [], 'b-', lw=2)
plot2, = ax2.plot([], [], 'b-', lw=2)

def draw_frame(n):
a = np.arange(n+1)
data = np.exp(-2j * np.pi * a/N)
plot1.set_data(data.real,data.imag)
plot2.set_data(range(n+1),a)
txt_title.set_text(f"data points : {n+1}")
return (plot1,)

anim = animation.FuncAnimation(fig, draw_frame, frames=N, interval=100, blit=True)

display.HTML(anim.to_html5_video())

[5]:


## Sine wave#

$$A_{shift} + A \sin(2{\pi}{f}{t} + \phi )$$

A = Amplitude
$$A_{shift}$$ = Amplitude Shift
f = frequency
t = time
$$\phi$$ = phase shift
[6]:

np.linspace(0,10,5) # evenly spaced 5 numbers between 0 to 10

[6]:

array([ 0. ,  2.5,  5. ,  7.5, 10. ])

[7]:

def sine_wave(n_samples, signal_freq, n_cycles = 10, amplitude = 1, amp_shift = 0, phase_shift = 0):
"""Sine wave generation
"""
sampling_freq = (n_samples * signal_freq) / n_cycles
# n_samples = int((sampling_freq / signal_freq) * n_cycles)

t_step = 1 / sampling_freq
time = np.linspace(0, (n_samples-1) * t_step , n_samples)

f_step = sampling_freq / n_samples
freq = np.linspace(0, (n_samples-1) * f_step , n_samples)

wave = amp_shift + amplitude * np.sin( (2 * np.pi * signal_freq * time) + phase_shift)

return wave, time, freq

[8]:

N = 100
wave, time, freq = sine_wave(
n_samples = N,
signal_freq = 10,
n_cycles = 1,
amplitude = 1,
amp_shift = 0,
phase_shift = 0
)

fig,ax = plt.subplots(1,2,figsize=(10,5))

ax[0].plot(time,wave,".-")
ax[1].plot(freq,np.abs(np.fft.fft(wave))/N)

plt.show()

[9]:

fig = plt.figure(figsize=(16,8))

wave, time, _ = sine_wave(
n_samples = 50,
signal_freq = 10,
n_cycles = 5
)
ax.plot(time, wave, 'r-*')

wave, time, _  = sine_wave(
n_samples = 40,
signal_freq = 10,
n_cycles = 5
)
ax.plot(time, wave, 'b-*')

wave, time, _  = sine_wave(
n_samples = 30,
signal_freq = 10,
n_cycles = 5
)
ax.plot(time, wave, 'g-*')

wave, time, _  = sine_wave(
n_samples = 10,
signal_freq = 10,
n_cycles = 5
)
ax.plot(time, wave, 'c-*')

plt.show()


## Sine Wave Animation#

[10]:

fig = plt.figure(figsize=(10,5))
ax1 = plt.subplot(1,2,(1,2))

ax1.set_xlim(( 0, 2))
ax1.set_ylim((-2, 2))

txt_title = ax1.set_title('Sine Wave')
plot1, = ax1.plot([], [], 'b.-', lw=2)

def draw_frame(n):
wave, time, _  = sine_wave(
n_samples = 50,
signal_freq = 10,
n_cycles = 5,
phase_shift=n
)
plot1.set_data(time,wave)
return (plot1,)

anim = animation.FuncAnimation(fig, draw_frame, frames=100, interval=100, blit=True)

display.HTML(anim.to_html5_video())

[10]:


## Discrete Fourier Transform#

$$X_{k} = \sum_{n=0}^{N-1} x_n e^{-{i}{2}{\pi}{k}\frac{n}{N}}$$

### Discrete fourier transform with iteration#

[11]:

def freq_plane(x,k,N):
return x * np.exp(-1j *2 * np.pi * k * (np.arange(N)/N))

def fourier_transform_old(x):
N = x.shape[0]
values = []
for k in range(N):
X_k = (freq_plane(x,k,N)).sum()
values.append(X_k)

X = np.array(values)
return X


but this is not optimized way to calculate.

I learnt this new method to use matrix multiplication and broadcasting in machine learning gradient descent calculation

### Optimized discrete fourier transform#

expanding the submission.. and looking for patterns to multiply \begin{align} X_0 &= e^\frac{i2\pi}{N}[ x_0 e^{0.0} + x_1 e^{0.1} + \dots + x_{N-1} e^{0.N-1} ]\\ X_1 &= e^\frac{i2\pi}{N}[ x_0 e^{1.0} + x_1 e^{1.1} + \dots + x_{N-1} e^{1.N-1} ]\\ \dots \\ X_{N-1} &= e^\frac{i2\pi}{N}[ x_0 e^{{N-1}.0} + x_1 e^{{N-1}.1} + ... + x_{N-1} e^{{N-1}{N-1}} ] \end{align}

Transforming in matrix

$$\begin{bmatrix} X_{0} & X_{1} \dots & X_{N-1} \end{bmatrix} = \\ \begin{bmatrix} x_{0} & x_{1} \dots & x_{N-1} \end{bmatrix} * e^{\frac{i2\pi}{N}} \begin{bmatrix} e^{0.0} & e^{0.1} & e^{0.2} & \dots & e^{0.N-1} \\ e^{1.0} & e^{1.1} & e^{1.2} & \dots & e^{1.N-1} \\ \dots \\ e^{N-1.0} & e^{N-1.1} & e^{N-1.2} & \dots & e^{N-1.N-1} \end{bmatrix}$$

Preparing $$k \times n$$

\begin{align} k &= \begin{bmatrix} k_0 \\ k_1 \\ k_2 \\ \dots \\ k_N-1 \end{bmatrix} \\ n &= \begin{bmatrix} n_0 & n_1 & n_2 & \dots & n_N-1 \end{bmatrix}\\ \text{Broadcasting }\\ \\ k \times n &= \begin{bmatrix} k_0.n_0 & k_0.n_1 & \dots & k_0.n_{N-1} \\ k_1.n_0 & k_1.n_1 & \dots & k_1.n_{N-1} \\ \dots \\ k_{N-1}.n_0 & k_{N-1}.n_1 & \dots & k_{N-1}.n_{N-1} \\ \end{bmatrix} \end{align}

[12]:

def fourier_transform(x):
"""Discrete Fourier Transform (DFT).

Notes:
This is faster and optimized approach.
using numpy matrix multiplication here.
"""

N = x.shape[0]
n = np.arange(N) # an array of size N-1
k = n.reshape(-1,1) # freq idxs with (N-1,1)

# k (N-1,1) is broadcasted over n (N-1,) and result is
m_value = np.exp(-2 * np.pi * k * (n/N) * 1j)
# (N-1,N-1)

X = x @ m_value
# (N-1,) * (N-1,N-1)
return X


### Winding plane#

here we are using an iterative approach to do the transform. because we want each plane for ever frequency to visualize

[13]:

N = 100
wave, time, freqs  = sine_wave(
n_samples = N,
signal_freq = 50,
n_cycles = 5
)

[14]:

plane = freq_plane(wave,0,N)
plt.plot(plane.real,plane.imag,'.-')
plt.show()

[15]:

plane = freq_plane(wave,1,N)
plt.plot(plane.real,plane.imag,'.-')
plt.show()

[16]:

plane = freq_plane(wave,5,N)
plt.plot(plane.real,plane.imag,'.-')
plt.show()

[17]:

N = 200
wave, time, freqs  = sine_wave(
n_samples = N,
signal_freq = 50,
n_cycles = 5
)

fig = plt.figure(figsize=(10,10))
ax1 = plt.subplot(2,2,1)
ax2 = plt.subplot(2,2,2)
ax3 = plt.subplot(2,2,3)
ax4 = plt.subplot(2,2,4)

ax1.set_xlim((-1.5, 1.5))
ax1.set_ylim((-1.5, 1.5))

ax2.set_xlim((-10, freqs[-1]))
ax2.set_ylim((-1e-12, 1e-12))

ax3.set_xlim((-10, freqs[-1]))
ax3.set_ylim((-200, 200))

ax4.set_xlim((-10, freqs[-1]))
ax4.set_ylim((-200, 200))

txt_title = fig.suptitle('This is a somewhat long figure title', fontsize=16)
ax1.set_title('Winding plane/ Complex plane')
ax2.set_title('Real Values')
ax3.set_title('Imaginary values')
ax4.set_title('Magnitude')
plot1, = ax1.plot([], [], 'g.-', lw=2, label='plane')
plot2, = ax1.plot([], [], 'k*', label='com')
plot3, = ax2.plot([], [])
plot4, = ax3.plot([], [])
plot5, = ax4.plot([], [])

ax1.legend()

real_com = []
imag_com = []
coms = []
i_com = []

def draw_frame(n):
plane = freq_plane(wave,n,N)
com = plane.sum()
real_com.append(com.real)
imag_com.append(com.imag)
coms.append(abs(com)/n)
i_com.append(freqs[n])

plot1.set_data(plane.real,plane.imag)
plot2.set_data(com.real,com.imag)
plot3.set_data(i_com,real_com)
plot4.set_data(i_com,imag_com)
plot5.set_data(i_com,coms)
txt_title.set_text(f"frequency : {freqs[n]}")
return (plot1,plot2,plot3,plot4,plot5)

anim = animation.FuncAnimation(fig, draw_frame, frames=N, interval=1000, blit=True)

display.HTML(anim.to_html5_video())

<ipython-input-17-358f6cd05818>:51: RuntimeWarning: divide by zero encountered in double_scalars
coms.append(abs(com)/n)
<ipython-input-17-358f6cd05818>:51: RuntimeWarning: divide by zero encountered in double_scalars
coms.append(abs(com)/n)

[17]: