# 來源：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 

# 建立 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()