# Discrete Fourier Transform¶

In [27]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib

In [28]:
nu = 4096
t = linspace(0, 2*pi, nu)
f = 3 + 4*sin(t)

# f is 3*cos(0) + -4*(-sin(2*pi*t/tau))
# zero, 1st harmonic cos,        zero, 1st harmonic sin

# exp(-2*pi*i*k*t/tau) = cos(2*pi*i*k*t/tau) - i sin(2*pi*i*k*t/tau)

f2 = np.real(fft.fft([3, 0], nu)) + np.imag(fft.fft([0, -4], nu))
#f2 = np.conj(np.real(fft.ifft([1, 0], nu)) + np.imag(fft.ifft([1, 0], nu)))*nu

plot(t, f)
plot(t, f2)
grid()

In [31]:
nu = 4
t = linspace(0, 2*pi, nu)
F = fft.fft(np.eye(nu))
print(np.round(F, 1))

[[ 1.+0.j  1.+0.j  1.+0.j  1.+0.j]
[ 1.+0.j  0.-1.j -1.+0.j  0.+1.j]
[ 1.+0.j -1.+0.j  1.+0.j -1.+0.j]
[ 1.+0.j  0.+1.j -1.+0.j  0.-1.j]]

In [32]:
data = np.array([1, 2, -2, 1])

In [33]:
ak = inv(F).dot(data)
ak

Out[33]:
array([ 0.5 +0.j  ,  0.75+0.25j, -1.  +0.j  ,  0.75-0.25j])
In [41]:
from numpy.fft import fft, ifft

f = fft([1, 0], n=4096)
plot(real(f))

Out[41]:
[<matplotlib.lines.Line2D at 0x7fb2018ef460>]
In [43]:
nu = 4
lam = exp(-1j*2*pi/nu)
z = lam**range(nu)
plt.plot(real(z), imag(z), 'x')
np.round(z, 2)
grid()
axis('equal');

In [45]:
def my_ift(x):
nu = len(x)
lam = exp(-1j*2*pi/nu)
F = np.zeros((nu, nu), dtype=complex)
ak_vect = np.hstack([np.arange(0, nu//2 + 1), -np.arange(nu//2 - (nu%2 == 0), 0, -1)])
j_vect = np.arange(nu)
for k in ak_vect:
for j in j_vect:
F[j, k] = lam**(k*j)
return np.conj(F).dot(x)/nu

In [46]:
sunspot = np.genfromtxt('./data/sunspots.csv',
plt.plot(sunspot[:, 0], sunspot[:, 1])
xlabel('t, year');
ylabel('intensity')
title('sunspot intensity')
grid()

In [47]:
T0 = sunspot[-1, 0] - sunspot[0, 0]
ak = my_ift(sunspot[:, 1])
k = arange(-50, 50)
bar(k, abs(ak[k])**2, width=0.8)
grid()
ylabel('$|a_k|^2$')
xlabel('freq cycles/year')

Out[47]:
Text(0.5, 0, 'freq cycles/year')
In [48]:
np.abs(ak[30])

Out[48]:
2.48506964342326
In [49]:
k_max = np.argmax(ak[1:])
f_max = k_max/T0
p_max = 1/f_max
f_max, p_max

Out[49]:
(0.09740259740259741, 10.266666666666666)
In [51]:
from datetime import datetime
"./data/global_temp.txt",
skiprows=22, dtype=float)
plt.plot(temp[:, 0], temp[:, 1])
plt.xlabel('year')
plt.ylabel('temp')
plt.grid()
plt.title('temperature')
temp;

In [52]:
T0 = temp[-1, 0] - temp[0, 0]
nu = len(temp[:, 0])
ak = my_ift(temp[:, 1])
k = np.arange(-20, 21)

a0_skip = 2
k_max = np.argmax(np.abs(ak[a0_skip:nu//2])**2) + a0_skip
freq_max = k_max/T0
print('dominant frequency', freq_max, 'cycles/year')
print('dominant period', 1/freq_max, 'years/cycle')

plt.bar(k/T0, np.abs(ak[k])**2, width=0.8/T0);
plt.ylabel('$|a_k|^2$')
plt.xlabel('cycles/year')
plt.grid()
plt.plot([-freq_max, freq_max],
np.abs(ak[k_max])**2*np.array([1, 1]), '*', markersize=10, color='y')
print(freq_max, ak[k_max])
plt.title('power spectrum')

dominant frequency 0.014925373134328358 cycles/year
dominant period 67.0 years/cycle
0.014925373134328358 (0.008510537691405189-0.1303208928383332j)

Out[52]:
Text(0.5, 1.0, 'power spectrum')

Here we want to ignore the $a_0$ and $a_1$ term since they are a constant and the first harmonic which is not very interesting. The $a_1$ power is high, but this is talking aboug the global trend of temperature increase

# Reconstructing Sinusoids from FFT/IFFT¶

In [62]:
np.linalg.norm(fft([1, 4, 3, 1], 4096) - np.conj(ifft([1, 4, 3, 1], 4096))*4096)

Out[62]:
0.0
In [63]:
fft([1, 2, 3], 4096) - np.conj(ifft([1, 2, 3], 4096))*4096

Out[63]:
array([0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j])

$\cos(\theta) - i \sin(\theta)$

$f = 3 + 10 \cos(3*t)$

In [91]:
nu = 4096
tau = 1
w0 = 2*pi/tau
t = linspace(0, tau, nu)
p = np.real(fft([3, 0, 0, 10], 4096))
plot(t, p, label='approx')
plot(t, 3 + 10*cos(3*t*w0), label='orig')

Out[91]:
[<matplotlib.lines.Line2D at 0x7fb1f0ab1a90>]
In [92]:
nu = 4096
tau = 1
w0 = 2*pi/tau
t = linspace(0, tau, nu)
p = np.real(ifft([3, 0, 0, 10], 4096))*4096
plot(t, p, label='approx')
plot(t, 3 + 10*cos(3*t*w0), label='orig')

Out[92]:
[<matplotlib.lines.Line2D at 0x7fb1f0e65400>]
In [95]:
nu = 4096
tau = 1
w0 = 2*pi/tau
t = linspace(0, tau, nu)
p = np.real(ifft([3], 4096))*4096 + np.imag(ifft([0, 0, 0, 10], 4096))*4096
plot(t, p, label='approx')
plot(t, 3 + 10*sin(3*t*w0), label='orig')

Out[95]:
[<matplotlib.lines.Line2D at 0x7fb1f072cca0>]
In [96]:
nu = 4096
t = np.linspace(0, 2*pi, nu)

plot(t, nu*imag(ifft([0, -1], nu)), linewidth=10, label='ifft');
plot(t, imag(fft([0, 1], nu)), linewidth=3, color='y', label='fft');
plt.legend()
plt.grid()


# Nyquist Sampling Theorem/ Aliasing¶

In [97]:
t = arange(0, 45*pi, pi/100)
T = 2*pi
dt = 1.1*T
t0 = 0.1

ts = arange(t0, 45*pi, dt)

figure(figsize=(15, 5))
plot(t, cos(t))
plot(ts, cos(ts), 'o-')

print('initial sample time: {:0.2g} seconds'.format(t0))
print('sampling interval: {:0.2g} seconds'.format(dt))

initial sample time: 0.1 seconds
sampling interval: 6.9 seconds


# Image Processing¶

In [98]:
def add_noise(data, n):
data2 = np.copy(data)
for i in range(n):
j = np.random.randint(0, data2.shape[0])
k = np.random.randint(0, data2.shape[1])
data2[j, k] = int(255*np.random.rand())
return data2

In [146]:
nu = 128
ak = np.zeros((nu, nu))
ak[0, 0] = 0
ak[1, 0] = 0
ak[0, 1] = 0
ak[1, 1] = 0
ak[3, 4] = 0
ak[1, 5] = 0
ak[29, 1] = 1
ak[10, 40] = 1
res = real(rfft2(ak))
figure(figsize=(10, 10))
imshow(res, cmap='Greys')
axis('tight');

In [122]:
import cv2
import numpy as np
from matplotlib import pyplot as plt

def analyze_image(file, noise=0, cx=30, cy=30):
img = cv2.imread(file, 0)  # read in image as greyscale, 0 means greyscale

f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)  # move a0 to the center of the image to make filtering easier

fshift_orig = np.copy(fshift)  # make a copy so we can modify fshift
mag_spectrum = 20*np.log10(np.abs(fshift_orig))  # magnitude in dB

rows, cols = img.shape
crow,ccol = rows//2 , cols//2

# gaussian kernel
colrow = np.array(np.meshgrid(np.arange(img.shape[1]), np.arange(img.shape[0])))
z = 1 - np.array(np.exp(-((colrow[1] - crow)**2/(2*cy**2) + (colrow[0] - ccol)**2/(2*cx**2))))

fshift_high = np.copy(fshift)

fshift_high = np.multiply(fshift_high, z)  # high pass filter using gaussian
#fshift_high[crow-cy:crow+cy, ccol-cx:ccol+cx] = 0  # high pass filter, set small ak's to zero

mag_spectrum_high = 20*np.log10(np.abs(fshift_high) + 1e-20)
f_ishift = np.fft.ifftshift(fshift_high)  # shift back s o in right form
img_high = np.fft.ifft2(f_ishift)

img_high = np.real(img_high)  # remove numerical complex numbers

img_low = img - img_high
mag_spectrum_low = 20*np.log10(np.abs(fshift_orig - fshift_high) + 1e-20)

vmin = 0
vmax = 255

print(np.min(img_low), np.max(img_low))
print(np.min(img_high), np.max(img_high))

plt.figure(figsize=(15, 4))
plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image')
plt.axis('tight')

plt.subplot(122),plt.imshow(mag_spectrum, cmap = 'gray', vmin=vmin, vmax=vmax)
plt.colorbar()
plt.axis('tight')
plt.title('Magnitude Spectrum, dB, $20\log(|a_k|)$')
plt.show()

plt.figure(figsize=(15, 4))
plt.subplot(121),plt.imshow(img_high, cmap = 'gray', vmin=-vmax, vmax=vmax)
plt.title('Image after High Pass Filter')
plt.axis('tight')

plt.subplot(122), plt.imshow(mag_spectrum_high, cmap = 'gray', vmin=vmin, vmax=vmax)
plt.colorbar()
plt.axis('tight')
plt.title('Magnitude Spectrum, dB, $20\log(|a_k|)$')

plt.figure(figsize=(15, 4))
plt.subplot(121),plt.imshow(img_low, cmap = 'gray', vmin=vmin, vmax=vmax)
plt.title('Image after Low Pass Filter')
plt.axis('tight')

plt.subplot(122), plt.imshow(mag_spectrum_low, cmap = 'gray', vmin=vmin, vmax=vmax)
plt.colorbar()
plt.axis('tight')
plt.title('Magnitude Spectrum, dB, $20\log(|a_k|)$')

In [124]:
analyze_image('./data/bigfoot.png', noise=100000, cx=30, cy=30)

50.48705388824312 199.23914838771958
-181.63053222806593 189.77568719281993