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))
ax = fig.add_subplot(111)
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)
# with numpy broadcasting
# 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)
# actually it is broadcasting
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]: