傅立葉變換用於分析各種濾波器的頻率特性。對於影象，二維離散傅立葉變換(2D Discrete Fourier Transform/DFT)用於尋找頻域。快速傅立葉變換(Fast Fourier Transform/FFT)的快速演算法用於計算DFT。

對於一個正弦訊號，x(t)=Asin(2πft)，我們可以說 f 是訊號的頻率，如果它的頻率域被接受，我們可以看到 f 的峰值。如果訊號被取樣來形成一個離散訊號，我們得到相同的頻率域，但是在[−π,π] or [0,2π]範圍內是週期性的（或者N-point的DFT，[0，N]）。你可以把影象看作一個訊號，從兩個方向進行取樣。所以在X和Y方向上進行傅立葉變換就會得到影象的頻率表示。

更直觀的是，對於正弦訊號，如果振幅在短時間內變化得非常快，你可以說它是一個高頻訊號。如果它變化緩慢，它是一個低頻訊號。你可以把同樣的想法擴充套件到圖片上。影象的振幅在哪裡變化？在邊緣點，或噪音。所以我們可以說，邊和噪聲是影象中的高頻內容。如果振幅沒有很大的變化，那就是低頻分量。

Numpy中的傅立葉變換

首先，我們將看到如何使用Numpy找到傅立葉變換。Numpy有一個FFT包來完成這個工作。np.fft.fft2()為我們提供了一個複雜陣列的頻率轉換。它的第一個引數是輸入影象，它是灰度影象。第二個引數是可選的，它決定了輸出陣列的大小。如果它大於輸入影象的大小，則輸入影象在計算FFT之前填充了0。如果它小於輸入影象，輸入影象將被裁剪。如果沒有引數傳遞，輸出陣列的大小將與輸入相同。

一旦得到了結果，零頻率元件(DC元件)將位於左上角。如果你想把它帶到中心，你需要在兩個方向上通過傳入 N/2 來改變結果。這是由函式np.fft.fftshift()完成的。(它更容易分析)。一旦你找到頻率變換，你就能找到大小譜。

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

img = cv.imread('messi5.jpg', 0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20 * np.log(np.abs(fshift))

plt.subplot(121
 
),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

你可以在中心看到更多的白色區域，表示低頻率的內容更多。

所以你找到了頻率變換。現在你可以在頻域做一些運算，比如高通濾波和重建影象，也就是找到逆DFT。為此，你只需用一個矩形視窗大小的60x60來移除低頻部分。然後使用np.fft.ifftshift()應用反向移動，使DC元件再次出現在左上角。然後使用np.ifft2()函式找到反FFT。再一次結果，將會是一個複數。你可以取它的絕對值。

rows, cols = img.shape
crow, ccol = rows/2, cols/2
fshift[crow-30:crow+30, ccol-30:ccol+30] = 0
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(img_back, cmap = 'gray')
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])

plt.show()

結果表明，高通濾波是一種邊緣檢測操作。這是我們在影象梯度章節中看到的。這也表明，大多數影象資料都存在於低頻區域。

如果您仔細觀察結果，特別是JET顏色的最後一個影象，您可以看到一些工件(我用紅色箭頭標記的一個例項)。它顯示了一些波紋狀的結構，叫做振鈴效應。它是由我們用來遮蔽的矩形視窗造成的。這個掩模被轉換成sinc形狀，導致了這個問題。因此，矩形視窗不用於過濾。更好的選擇是高斯視窗。

OpenCV中的傅立葉變換

OpenCV提供了cv.dft()和cv.idft()函式。它返回與前面相同的結果，但是有兩個通道。第一個通道將會有結果的實部，第二個通道將會有一個虛部。輸入影象首先應該轉換為np.float32。我們將會看到如何去做。

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

img = cv.imread('messi5.jpg', 0)

dft = cv.dft(np.float32(img), flags=cv.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)

magnitude_spectrum = 20 * np.log(cv.magnitude(dft_shift[:,:,0], dft_shift[:,:,1]))

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

你也可以使用cv.cartToPolar()，它可以在一次拍攝中同時返回大小和相位。
現在我們要做的是逆DFT。在前面，我們建立了一個HPF，這次我們將看到如何移除影象中的高頻內容，即我們將LPF應用到影象中。它實際上模糊了影象。為此，我們先建立一個具有高值(1)低頻率的掩碼，即我們通過低頻內容，而在高頻區域則是0。

rows, cols = img.shape
crow, ccol = rows/2, cols/2

# 首先建立一個掩碼，中心方為1，其餘為0
mask = np.zeros((rows, cols, 2), np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1

# 應用掩碼，反DFT
fshift = dft_shift * mask
f_ishift = np.fft.ifftshift(fshift)
img_back = cv.idft(f_ishift)
img_back = cv.magnitude(img_back[:,:,0], img_back[:,:,1])

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('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

效能優化的DFT

對於某些陣列的大小來說，DFT的計算效能更好。當陣列大小為2時，它是最快的。這些陣列的大小是2、3和5的乘積，它們的處理效率也相當高。因此，如果您擔心程式碼的效能，您可以在找到DFT之前將陣列的大小修改為任何最佳大小(通過填充零)。對於OpenCV，你必須手動填充0。但是對於Numpy，您指定了FFT計算的新大小，它將自動為您填充0。

那麼如何找到這個最優的大小呢?OpenCV為這個提供了一個函式，cv.getOptimalDFTSize()。它適用於cvdft()和np.fft.fft2()。

In [16]: img = cv.imread('messi5.jpg',0)
In [17]: rows,cols = img.shape
In [18]: print("{} {}".format(rows,cols))
342 548

In [19]: nrows = cv.getOptimalDFTSize(rows)
In [20]: ncols = cv.getOptimalDFTSize(cols)
In [21]: print("{} {}".format(nrows,ncols))
360 576

大小(342,548)被修改為(360，576)。現在讓我們用0(對於OpenCV)來填充它，並找到它們的DFT計算效能。您可以通過建立一個新的大型零陣列並將資料複製到它，或者使用cv.copyMakeBorder()來實現它。

nimg = np.zeros((nrows, ncols))
nimg[:rows, :cols] = img

OR

right = ncols - cols
bottom = nrows - rows
bordertype = cv.BORDER_CONSTANT  # 只是為了避免PDF檔案中的行拆分
nimg = cv.copyMakeBorder(img, 0, bottom, 0, right, bordertype, value=0)

現在我們計算DFT NumPy功能效能比較：

In [22]: %timeit fft1 = np.fft.fft2(img)
10 loops, best of 3: 40.9 ms per loop
In [23]: %timeit fft2 = np.fft.fft2(img,[nrows,ncols])
100 loops, best of 3: 10.4 ms per loop

它顯示了一個4倍的加速。現在我們來試試OpenCV函式。

In [24]: %timeit dft1= cv.dft(np.float32(img),flags=cv.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 13.5 ms per loop
In [27]: %timeit dft2= cv.dft(np.float32(nimg),flags=cv.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 3.11 ms per loop

它也顯示了一個4倍的加速。您還可以看到，OpenCV函式的速度比Numpy函式快了大約3倍

為什麼Laplacian是一個高通濾波器?

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

# 無標度引數的簡單平均濾波器
mean_filter = np.ones((3,3))

# 建立一個高斯濾波器
x = cv.getGaussianKernel(5,10)
gaussian = x*x.T

# 不同的邊緣檢測濾波器
# scharr in x-direction
scharr = np.array([[-3, 0, 3],
                   [-10,0,10],
                   [-3, 0, 3]])

# sobel in x direction
sobel_x= np.array([[-1, 0, 1],
                   [-2, 0, 2],
                   [-1, 0, 1]])

# sobel in y direction
sobel_y= np.array([[-1,-2,-1],
                   [0, 0, 0],
                   [1, 2, 1]])

# laplacian
laplacian=np.array([[0, 1, 0],
                    [1,-4, 1],
                    [0, 1, 0]])

filters = [mean_filter, gaussian, laplacian, sobel_x, sobel_y, scharr]
filter_name = ['mean_filter', 'gaussian','laplacian', 'sobel_x', \
                'sobel_y', 'scharr_x']
fft_filters = [np.fft.fft2(x) for x in filters]
fft_shift = [np.fft.fftshift(y) for y in fft_filters]
mag_spectrum = [np.log(np.abs(z)+1) for z in fft_shift]

for i in xrange(6):
    plt.subplot(2,3,i+1),plt.imshow(mag_spectrum[i],cmap = 'gray')
    plt.title(filter_name[i]), plt.xticks([]), plt.yticks([])

plt.show()