%pylab inline
%pylab is deprecated, use %matplotlib inline and import the required libraries. Populating the interactive namespace from numpy and matplotlib
ak = ifft(f)
freq = fftfreq(nu, 1)
plot(freq[0:10], ak[0:10])
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) Input In [2], in <cell line: 1>() ----> 1 ak = ifft(f) 2 freq = fftfreq(nu, 1) 3 plot(freq[0:10], ak[0:10]) File <__array_function__ internals>:180, in ifft(*args, **kwargs) File ~/anaconda3/envs/aae301/lib/python3.9/site-packages/numpy/fft/_pocketfft.py:314, in ifft(a, n, axis, norm) 312 a = asarray(a) 313 if n is None: --> 314 n = a.shape[axis] 315 inv_norm = _get_backward_norm(n, norm) 316 output = _raw_fft(a, n, axis, False, False, inv_norm) IndexError: tuple index out of range
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()
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]]
data = np.array([0, 1, 0, -1])
ak = inv(F).dot(data)
ak
array([0.+0.j , 0.+0.5j, 0.+0.j , 0.-0.5j])
from numpy.fft import fft, ifft
f = fft([1, 0], n=4096)
plot(real(f))
[<matplotlib.lines.Line2D at 0x7f5f4a9e6730>]
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');
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
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()
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')
Text(0.5, 0, 'freq cycles/year')
np.abs(ak[30])
2.4850696434232598
k_max = np.argmax(ak[1:])
f_max = k_max/T0
p_max = 1/f_max
f_max, p_max
(0.09740259740259741, 10.266666666666666)
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;
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)
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
np.linalg.norm(fft([1, 4, 3, 1], 4096) - np.conj(ifft([1, 4, 3, 1], 4096))*4096)
0.0
fft([1, 2, 3], 4096) - np.conj(ifft([1, 2, 3], 4096))*4096
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)$
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')
[<matplotlib.lines.Line2D at 0x7f5f48d4e9d0>]
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')
[<matplotlib.lines.Line2D at 0x7ff0cf31c970>]
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')
[<matplotlib.lines.Line2D at 0x7ff0cf3a12e0>]
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()
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
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
nu = 128
ak = np.zeros((nu, nu))
ak[0, 0] = 0
ak[1, 0] = 1
ak[0, 1] = 0
ak[1, 1] = 0
ak[3, 4] = 1
ak[1, 5] = 1
ak[29, 1] = 1
ak[10, 40] = 1
res = real(rfft2(ak))
figure(figsize=(10, 10))
imshow(res, cmap='Greys')
axis('tight');
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|)$')
analyze_image('./data/bigfoot.png', noise=100000, cx=30, cy=30)
52.94940125100608 191.76169220701556 -177.67835123574312 189.59257964927286
analyze_image('./data/nessy.jpg', noise=0, cx=50, cy=50)
-0.00021607081899269867 174.08057129632274 -142.27140806315464 163.62149821958812
analyze_image('./data/lizard_noisy.png', noise=0, cx=20, cy=20)
66.8865988453692 197.5550051528936 -193.4932269813338 187.59790402644256