NumPy 中的傅里葉分析
# 來源：NumPy Essentials ch6
繪圖函數
import matplotlib.pyplot as plt
import numpy as np
def show(ori_func, ft, sampling_period = 5):
n = len(ori_func)
interval = sampling_period / n
# 繪制原始函數
plt.subplot(2, 1, 1)
plt.plot(np.arange(0, sampling_period, interval), ori_func, 'black')
plt.xlabel('Time'), plt.ylabel('Amplitude')
# 繪制變換后的函數
plt.subplot(2,1,2)
frequency = np.arange(n / 2) / (n * interval)
nfft = abs(ft[range(int(n / 2))] / n )
plt.plot(frequency, nfft, 'red')
plt.xlabel('Freq (Hz)'), plt.ylabel('Amp. Spectrum')
plt.show()
信號處理
# 生成頻率為 1（角速度為 2 * pi）的正弦波
time = np.arange(0, 5, .005)
x = np.sin(2 * np.pi * 1 * time)
y = np.fft.fft(x)
show(x, y)
# 將其與頻率為 20 和 60 的波疊加起來
x2 = np.sin(2 * np.pi * 20 * time)
x3 = np.sin(2 * np.pi * 60 * time)
x += x2 + x3
y = np.fft.fft(x)
show(x, y)
# 生成方波，振幅是 1崖咨，頻率為 10Hz
# 我們的間隔是 0.05s非凌，每秒有 200 個點
# 所以需要每隔 20 個點設為 1
x = np.zeros(len(time))
x[::20] = 1
y = np.fft.fft(x)
show(x, y)
# 生成脈沖波
x = np.zeros(len(time))
x[380:400] = np.arange(0, 1, .05)
x[400:420] = np.arange(1, 0, -.05)
y = np.fft.fft(x)
show(x, y)
# 生成隨機數
x = np.random.random(100)
y = np.fft.fft(x)
show(x, y)
原理
x = np.random.random(500)
n = len(x)
m = np.arange(n)
k = m.reshape((n, 1))
M = np.exp(-2j * np.pi * k * m / n)
y = np.dot(M, x)
np.allclose(y, np.fft.fft(x))
# True
'''
%timeit np.dot(np.exp(-2j * np.pi * np.arange(n).reshape((n, 1)) * np.arange(n) / n), x)
10 loops, best of 3: 18.5 ms per loop
%timeit np.fft.fft(x)
100000 loops, best of 3: 10.9 μs per loop
'''
# 傅里葉逆變換
M2 = np.exp(2j * np.pi * k * m / n)
x2 = np.dot(y, M2) / n
np.allclose(x, x2)
# True
np.allclose(x, np.fft.ifft(y))
# True
# 創(chuàng)建 10 個 0~9 隨機整數的信號
a = np.random.randint(10, size = 10)
a
# array([7, 4, 9, 9, 6, 9, 2, 6, 8, 3])
a.mean()
# 6.2999999999999998
# 進行傅里葉變換
A = np.fft.fft(a)
A
'''
array([ 63.00000000 +0.00000000e+00j,
-2.19098301 -6.74315233e+00j,
-5.25328890 +4.02874005e+00j,
-3.30901699 -2.40414157e+00j,
13.75328890 -1.38757276e-01j,
1.00000000 -2.44249065e-15j,
13.75328890 +1.38757276e-01j,
-3.30901699 +2.40414157e+00j,
-5.25328890 -4.02874005e+00j,
-2.19098301 +6.74315233e+00j])
'''
A[0] / 10
# (6.2999999999999998+0j)
A[int(10 / 2)]
# (1-2.4424906541753444e-15j)
# A[0] 是 0 頻率的項
# A[1:n/2] 是正頻率項
# A[n/2 + 1: n] 是負頻率項
# 如果我們要把 0 頻率項調整到中間
# 可以調用 fft.fftshift
np.fft.fftshift(A)
'''
array([ 1.00000000 -2.44249065e-15j,
13.75328890 +1.38757276e-01j,
-3.30901699 +2.40414157e+00j,
-5.25328890 -4.02874005e+00j,
-2.19098301 +6.74315233e+00j,
63.00000000 +0.00000000e+00j,
-2.19098301 -6.74315233e+00j,
-5.25328890 +4.02874005e+00j,
-3.30901699 -2.40414157e+00j,
13.75328890 -1.38757276e-01j])
'''
# fft2 用于二維抄谐，fftn 用于多維
x = np.random.random(24)
x.shape = 2,12
y2 = np.fft.fft2(x)
x.shape = 1,2,12
y3 = np.fft.fftn(x, axes = (1, 2))
np.allclose(y2, y3)
# True
應用
from matplotlib import image
# 將上面的圖片保存為 scientist.png
# 并讀入
img = image.imread('./scientist.png')
# 將圖片轉換為灰度圖
# 每個像素是 0.21R + 0.72G + 0.07B
gray_img = np.dot(img[:,:,:3], [.21, .72, .07])
gray_img.shape
# (317L, 661L)
plt.imshow(gray_img, cmap = plt.get_cmap('gray'))
# <matplotlib.image.AxesImage at 0xa6165c0>
plt.show()
# fft2 是二維數組的傅里葉變換
# 將空域轉換為頻域
fft = np.fft.fft2(gray_img)
amp_spectrum = np.abs(fft)
plt.imshow(np.log(amp_spectrum))
# <matplotlib.image.AxesImage at 0xcdeff60>
plt.show()
fft_shift = np.fft.fftshift(fft)
plt.imshow(np.log(np.abs(fft_shift)))
# <matplotlib.image.AxesImage at 0xd201dd8>
plt.show()
# 放大圖像
# 我們向 fft_shift 插入零頻率嫌松，將其尺寸擴大兩倍
m, n = fft_shift.shape
b = np.zeros((int(m / 2), n))
c = np.zeros((2 * m - 1, int(n / 2)))
fft_shift = np.concatenate((b, fft_shift, b), axis = 0)
fft_shift = np.concatenate((c, fft_shift, c), axis = 1)
# 然后再轉換回去
ifft = np.fft.ifft2(np.fft.ifftshift(fft_shift))
ifft.shape
# (633L, 1321L)
ifft = np.real(ifft)
plt.imshow(ifft, cmap = plt.get_cmap('gray'))
# <matplotlib.image.AxesImage at 0xf9a0f98>
plt.show()
