快速傅立葉變換的應用非常廣。用FFT(快速傅立葉變換)能將時域的數字訊號轉換為頻域訊號。轉換為頻域訊號之後我們可以很方便地分析出訊號的頻率成分，在頻域上進行處理，最終還可以將處理完畢的頻域訊號通過IFFT(逆變換)轉換為時域訊號，實現許多在時域無法完成的訊號處理演算法。

　　在python的scipy類庫中提供了快速傅立葉變換包fftpack。下面的例子中低頻的原始訊號加入了高頻噪聲，我們可以通過快速傅立葉變換來分析其頻譜，去掉噪聲再重構訊號。

import numpy as np
from scipy import fftpack  # The scipy.fftpack module allows to compute fast Fourier transforms
 

import pylab as pl

time_step = 0.02   # 訊號取樣間隔Δt
period = 5.        # 正弦訊號週期T
time_vec = np.arange(0, 20, time_step) # 生成取樣點,取樣次數N=20

# 模擬帶高頻噪聲的訊號
sig = np.sin(2 * np.pi / period * time_vec) + 0.5 * np.random.randn(time_vec.size)

'''
fftfreq函式返回離散頻率序列f, 頻域間隔:Δf=1/(N*Δt)
scipy.fftpack.fftfreq(n, d=1.0)
f = [0, 1, ...,   n/2-1,     -n/2, ..., -1] / (d*n)   if n is even
f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n)   if n is odd
 
'''
sample_freq = fftpack.fftfreq(sig.size, d=time_step) 

sig_fft = fftpack.fft(sig) # compute the fast Fourier transform

pidxs = np.where(sample_freq > 0) # only the positive part of the spectrum needs to be used
freqs, power = sample_freq[pidxs], np.abs(sig_fft)[pidxs]
freq = freqs[power.argmax()]

 
# 去除高頻噪聲
sig_fft[np.abs(sample_freq) > freq] = 0

# 用傅立葉反變換，重構去除噪聲的訊號
main_sig = fftpack.ifft(sig_fft)

pl.figure()
pl.plot(freqs, power)
pl.xlabel('Frequency [Hz]')
pl.ylabel('plower')
axes = pl.axes([0.3, 0.3, 0.5, 0.5])
pl.title('Peak frequency')
pl.plot(freqs[:8], power[:8])
pl.setp(axes, yticks=[])

pl.figure()
pl.plot(time_vec, sig)
pl.plot(time_vec, main_sig, linewidth=3)
pl.xlabel('Time [s]')
pl.ylabel('Amplitude')
pl.show()

　　下圖可以看出在原始訊號頻率0.2Hz處，幅值最大，其它高於0.2Hz的都為噪聲訊號。用FFT方法消除噪聲就是對含噪聲訊號的頻譜進行處理，將噪聲所在頻段的$X(k)$值全部置零後，再對處理後的$X(k)$進行離散傅立葉逆變換(IFFT)可得原始訊號的近似結果。這種方法要求噪聲與訊號的頻譜不在同一頻段，否則將很難處理。

　　下圖中的粗線代表通過fft及ifft從帶噪聲的訊號中提取的有用訊號

• 二維離散傅立葉變換

　　一個影象為M×N的函式f(x,y)的離散傅立葉變換F(u,v)為:

$$F(u,v)=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi (ux/M+vy/N)}$$

　　(u,v)=(0,0)位置的傅立葉變換為$F(0,0)=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)=\bar{f}(x,y)$，即f(x,y)的均值，原點(0,0)的傅立葉變換是影象的平均灰度。F(0,0)稱為頻率譜的直流分量，其它F(u,v)稱為交流分量。傅立葉變換中出現的變數u和v通常稱為頻率變數，空間頻率可以理解為等相位線在x,y座標軸投影的截距的倒數。

　　相應的空間頻率分別為：$u=\frac{1}{X},\quad v=\frac{1}{Y}$，對影象訊號而言，空間頻率是指單位長度內亮度作週期性變化的次數。

　　為了更好理解影象的傅立葉變換用下面幾幅圖來進行說明：

　　下圖右邊是原始影象，左邊是通過傅立葉變換得到的影象（沒有移頻，即原點在影象的左上角）。由於原始影象的灰度只在Y軸上發生週期性變化，因此對應的傅立葉影象中只在左側邊緣出現一個亮點，代表原始影象Y軸上灰度變化的頻率。如果右側影象水平黑白條紋越密集（變化越快），則左側影象上的亮點越遠離原點向下；若右側水平黑白條紋越稀疏（變化越慢），則左側影象上的亮點越接近原點。

　　下面這幅圖右側的原始影象在X軸和Y軸上都有灰度的變化，因此對應的傅立葉影象上的亮點不再處於座標軸上（邊緣上）：

　　如果原始影象中存在兩種頻率，則顯然其傅立葉變換的影象中會出現兩個對應的離散點：

• 影象傅立葉變換的理解　

　　把影象想象成沿著兩個方向採集的訊號。所以對影象同時進行X方向和Y方向的傅立葉變換，我們就會得到這幅影象的頻域表示（頻譜圖），也就是影象梯度的分佈圖，當然頻譜圖上的各點與影象上各點並不存在一一對應的關係。二維傅立葉變換的座標軸意義是頻率，越靠近原點，頻率越低，對應於影象中畫素值變化速度比較慢的部分，一般是物體的主體、背景等。越遠離原點，頻率越高，對應於影象中畫素值變化速度快的那部分，例如物體的邊界。對一般影象來說靠近原點周圍比較亮，遠離原點的地方比較暗，也就是影象裡低頻部分分量多，高頻分量少。一張圖片中顯然是邊緣部分少，而大部分都是顏色相近/灰度相近的區域。

　　不同頻率資訊在影象結構中有不同的作用。影象的主要成分是低頻資訊，它形成了影象的基本灰度等級，對影象結構的決定作用較小；中頻資訊決定了影象的基本結構，形成了影象的主要邊緣結構；高頻資訊形成了影象的邊緣和細節，是在中頻資訊上對影象內容的進一步強化。傅立葉變換提供另外一個角度來觀察影象，可以將影象從灰度分佈轉化到頻率分佈上來觀察影象的特徵。

• 利用傅立葉變換進行濾波

　　Numpy的FFT包可以幫助我們實現快速傅立葉變換。函式np.fft.fft2()可以對訊號進行頻率轉換，輸出結果是一個複雜的陣列。函式的第一個引數是輸入影象，要求是灰度格式。第二個引數是可選的，決定輸出陣列的大小。預設輸出陣列的大小和輸入影象大小一樣；如果輸出結果比輸入影象大，輸入影象就需要在進行FFT前補0；如果輸出結果比輸入影象小的話，輸入影象就會被切割。

import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread('input.png',0)

# fft2() computes the n-dimensional discrete Fourier Transform
f = np.fft.fft2(img)

# Shift the zero-frequency component to the center of the spectrum
fshift = np.fft.fftshift(f)

magnitude_spectrum = 20*np.log(np.abs(fshift))

rows, cols = img.shape
crow, ccol = rows/2 , cols/2
fshift[crow-40:crow+40, ccol-40:ccol+40] = 0  # remove the low frequencies by masking with a rectangular window of size 80x80
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)

plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back, cmap = 'gray')
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.show()

　　上圖的結果顯示高通濾波其實是一種邊界檢測操作。

　　OpenCV中相應的函式是cv2.dft()和cv2.idft()，和前面輸出的結果一樣，但是是雙通道的。第一個通道是結果的實數部分，第二個通道是結果的虛數部分，輸入影象要首先轉換成np.float32格式。OpenCV中的函式cv2.dft() 和cv2.idft()要比Numpy快，但是Numpy函式更加使用者友好。

　　當陣列的大小為某些值時DFT的效能會更好。當陣列的大小是2的冪時DFT效率最高；當陣列的大小是 2，3，5的倍數時效率也會很高。所以如果你想提高程式碼的執行效率時，你可以修改輸入影象的大小（補0）。對於OpenCV必須自己手動補0，但是Numpy只需要指定FFT運算的大小，它會自動補 0。那我們怎樣確定最佳大小呢？ OpenCV提供了一個函式:cv2.getOptimalDFTSize()。它可以同時被cv2.dft()和np.fft.fft2()使用。

　　讓我們使用 IPython中的魔法命令%timeit 來測試一下：

　　陣列的大小從(186,295)變成了(192,300)，現在我們為它補 0，然後看看效能有沒有提升。

　　現在我們看看Numpy的表現：

　　速度提高了3倍。我們再看看OpenCV的表現：

　　我們也會發現OpenCV的速度比Numpy還要快2倍。

　　在前面的部分我們實現了一個HPF（高通濾波），現在我們來實現LPF（低通濾波）將高頻部分去除。其實就是對影象進行模糊操作。首先我們需要構建一個掩模，與低頻區域對應的地方設定為 1, 與高頻區域對應的地方設定為 0。

import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread('input.png',0)

rows, cols = img.shape
nrows = cv2.getOptimalDFTSize(rows) # Returns the optimal DFT size for a given vector size
ncols = cv2.getOptimalDFTSize(cols)
nimg = np.zeros((nrows,ncols))
nimg[:rows,:cols] = img

dft = cv2.dft(np.float32(nimg), flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)

crow, ccol = nrows/2 , ncols/2

# create a mask first, center square is 1, remaining all zeros
mask = np.zeros((nrows, ncols, 2), np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1

# apply mask and inverse DFT
fshift = dft_shift * mask

f_ishift = np.fft.ifftshift(fshift)

# Computes the inverse Discrete Fourier Transform of a 1D or 2D array
img_back = cv2.idft(f_ishift)  

# You can also use cv2.cartToPolar() which returns both magnitude and phase
img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1]).reshape((nrows,ncols)) 
img_back = img_back[:rows,:cols]

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_back, cmap = 'gray')
plt.title('Image after LPF'), plt.xticks([]), plt.yticks([])
plt.show()

　　下面的一個例子利用二維傅立葉變換來消除影象上存在的週期性噪聲。週期噪聲一般產生於影象採集過程中的電氣或電機的干擾，表現為影象中週期性的衝擊，如下圖：

　　通過觀察傅立葉變換後的頻譜圖，我們就可以看出影象的能量分佈，如果頻譜圖中暗的點數更多，那麼實際影象是比較柔和的（因為各點與鄰域差異都不大，梯度相對較小）；反之，如果頻譜圖中亮的點數多，那麼實際影象一定是尖銳的，邊界分明且邊界兩邊畫素差異較大的。對頻譜移頻到影象中心以後，可以看出影象的頻率分佈是以中心為圓心，對稱分佈的。將頻譜移頻到中心除了可以清晰地看出影象頻率分佈以外，還有一個好處，它可以分離出有周期性規律的干擾訊號。比如正弦干擾，一幅帶有正弦干擾，移頻到中心的頻譜圖上可以看出除了中心以外還存在以某一點為中心，對稱分佈的亮點集合，這個集合就是干擾噪音產生的，這時可以很直觀的通過在該位置放置帶阻濾波器消除干擾。

　　若原點設在中心，其頻譜能量集中分佈在中心附近；若原點設在左上角，那麼影象訊號能量將集中在四個角上。這是由二維傅立葉變換本身性質決定的。變換之後的影象在原點平移之前四角是低頻，最亮，平移之後中間部分是低頻，最亮，亮度大說明低頻的能量大。

　　這張登月照片上面有著很有規律的噪點，那麼其FFT頻譜圖上面就會有非常規則的點。這些點就是噪聲在頻域空間的對應。

　　我們需要保留頻譜影象中有用的部分，去掉高頻噪聲。下面的Python程式碼中設定了一個比例，保留這些有用資訊，將高頻區域的頻率設為0，消除噪聲影響。注意程式中沒有進行移頻操作。

import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread('moonlanding.png',0)

rows, cols = img.shape
nrows = cv2.getOptimalDFTSize(rows) # Returns the optimal DFT size for a given vector size
ncols = cv2.getOptimalDFTSize(cols)

f = np.fft.fft2(img, [nrows,ncols])

magnitude_spectrum_1 = 20 * np.log(1 + np.abs(f))

# Define the fraction of coefficients (in each direction) we keep
keep_fraction = 0.1

# Set to zero all rows with indices between row*keep_fraction and # row*(1-keep_fraction):
f[nrows*keep_fraction:nrows*(1-keep_fraction)] = 0

# Similarly with the columns:
f[:, ncols*keep_fraction:ncols*(1-keep_fraction)] = 0

magnitude_spectrum_2 = 20 * np.log(1 + np.abs(f))

img_back = np.abs(np.fft.ifft2(f))
img_back = img_back[:rows,:cols]

plt.subplot(221),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(222),plt.imshow(magnitude_spectrum_1, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.subplot(223),plt.imshow(img_back, cmap = 'gray')
plt.title('Reconstructed Image'), plt.xticks([]), plt.yticks([])
plt.subplot(224),plt.imshow(magnitude_spectrum_2, cmap = 'gray')
plt.title('Filtered Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

　　在Mathematica軟體中可以很方便的用函式ImagePeriodogram[image]來顯示影象的頻譜(shows the squared magnitude of the discrete Fourier transform (power spectrum) of image)。在該函式中使用可選引數Alignment可以設定零頻率的位置，預設值是Alignment→Center，即影象中心頻率最低。如果設定為Alignment→{Left,Top}，那麼零頻分量將被放置在影象左上角。

參考：

http://old.sebug.net/paper/books/scipydoc/frequency_process.html

Scipy : high-level scientific computing

SciPy Lecture Notes 中文版

影象傅立葉變換

傅立葉變換在影象處理中的應用

知乎-傅立葉變換有哪些具體的應用？

How to use 2D Fourier analysis to clean the noise in an image

影象的傅立葉變換

An Intuitive Explanation of Fourier Theory

What does frequency domain denote in case of images?