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',
              delimiter=',', skip_header=1)
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
temp = np.loadtxt(
    "./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
    img = add_noise(img, noise)   # add noise to image

    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