圖像的傅里葉變換
傅里葉變換（Fourier Transform）是非常重要的數學分析工具，同時也是一種非常重要的信號處理方法。我記得本科課程電路原理中有提到過，但由于計算過于復雜，好像是超出考試范圍了，所以并沒有深入學習。最近實驗中需要對圖像進行濾波處理，文獻中提到的方法通常是經過傅里葉變換之后對頻域進行過濾，將圖像中的低頻信息與高頻信息區分開來。
理解傅里葉變換對非數學專業的人來說比較難以理解的原因主要有兩方面。首先，由于涉及到比較復雜繁瑣的數學操作，看到下面的這兩個公式，一般人可能當時就蒙了：
$$ \widehat{f}(t) = \int_{-\infty}^{\infty} f(x) e^{-2\pi ixt}dx; $$
$$ f(x) = \int_{-\infty}^{\infty}\widehat f (t) e^{2\pi itx}dt;
$$
另一方面的原因則是由于變換過程比較抽象，很難從直覺上去把握在傅里葉變換過程中到底發生了什么。關于傅里葉變換的科普文章，有一篇是知乎上傳閱較廣的《傅里葉分析之掐死教程》，作者用了盡可能少的數學公式和圖形分解來解決這兩方面的問題。在此之前，國外有個專門科普數學概念的網站（Better Explained）也寫了一篇類似的科普文章，但更徹底的是，全文都不涉及到任何數學公式和推斷，完全用英語來向讀者解釋傅里葉變換過程。第一篇文章以音樂和樂譜為例，第二篇作者用“奶昔”作為例子，我想了半天終于找到一個更通俗的例子：
我們來想象一下，假設這塊面餅的厚度是$N$層面條，每一層都是由一根彎曲成正弦曲線形狀的面條排列而成，有些面條波浪較大，也就是排列較為稀疏，而密有些排列較密集。我們把這$N$層面條擠壓到一起，就得到上圖這一塊雜亂無序、世間獨一無二的面餅。傅里葉變換所做的事就是把上面的過程反過來，我們可以從一塊完整的面餅得到最初的$N$層面條。如果有些人的口味比較特殊，喜歡波浪大又稀疏的面，于是我們就將排列太緊湊（高頻面）剔除之后再重新壓制成一塊新的面餅，這就是我們最終想要的濾波（Filter）的過程。
帶著對面餅的想象，我們來看一種更為抽象、優雅的描述（from Wikipedia）：
一般的波形或者說信號（Signal）都是基于時間尺度上的采樣結果，因此也稱為時域（Time Domain），而上面泡面的例子和我們將要處理的圖像信號則是基于空間尺度上的采樣，但好像并沒有“空域（Space Domain）”這一說，畢竟我們對空間的感知仍然依賴于時間。不過在空間尺度上我們可以更直觀地認為信號是靜止，例如下面這張圖像（灰度圖），其實是由 250x250個像素點組成，每個像素點的灰度值（$[0, 255]$）就是基于像素坐標的空間采樣的結果：
右邊的3D Fourier就是一塊長相奇怪的面餅。
實踐
對傅里葉變換有了大概的了解之后可以先動手嘗試一下，來更加直觀地感受一下（實際上完全可以在不理解的情況下，直接上手）。這里用到的是OpenCV + Numpy，實際上OpenCV和Numpy都提供了快速傅里葉變換（FFT）算法：
import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
img = cv.imread('Joseph_Fourier_250.jpg', 0)
f = np.fft.fft2(img) # 快速傅里葉變換算法得到頻率分布
fshift = np.fft.fftshift(f) # 默認結果中心點位置是在左上角，轉移到中間位置
fimg = np.log(np.abs(fshift)) # fft 結果是復數，求絕對值結果才是振幅
# 展示結果
plt.subplot(121), plt.imshow(img, 'gray'), plt.title('Original Fourier')
plt.subplot(122), plt.imshow(fimg, 'gray'), plt.title('Fourier Fourier')
plt.show()
右邊圖就是頻率分布圖譜，其中越靠近中心的位置頻率越低，越亮（灰度值越高）的位置代表該頻率的信號振幅越大。fft的結果是復數形式，保留了圖像的全部信息，但去絕對值得到的頻譜圖只表現了振幅而沒有體現相位。
回想一下高中時候學過的三角函數：
$$ sin(x) = A(\omega x+\varphi) = A(2\pi fx + \varphi)
$$
一個正弦波是由下面三個參數決定的：
- 角速度（頻率）$\omega = 2\pi f$ ；
- 振幅$A$；
- 相位$\varphi$。
除了上面這個公式之外，還可以用另外一種形式來（唯一地）表示一個正弦波（from BetterExplained）：
即：
$$ cos(x) + i sin(x) \Leftrightarrow a + i b
$$
所以說，fft的復數結果保留了正弦波成分的所有信息，但頻譜圖只展現了頻率和振幅的分布。因此可以根據fft的結果還原原始圖像，但是我們做傅里葉變換的目的并不是為了觀察圖像的頻率分布（至少不是最終目的），更多情況下是為了對頻率進行過濾。過濾的方法一般有三種：低通（Low-pass）、高通（High-pass）、帶通（Band-pass）。所謂低通就是保留圖像中的低頻成分，過濾高頻成分，可以把過濾器想象成一張漁網，根據上文對頻譜圖的解讀，想要低通過濾器，就是將高頻區域的信號全部拉黑，而低頻區域全部保留：
img = cv.imread('Joseph_Fourier_250.jpg', 0)
f = np.fft.fft2(img) # 快速傅里葉變換算法得到頻率分布
fshift = np.fft.fftshift(f) # 默認結果中心點位置是在左上角，轉移到中間位置
lpButterMask = butterWorthLPF(img.shape[:2], 24, 2)
hpButterMask = butterWorthHPF(img.shape[:2], 36, 2)
lpFshift = fshift * lpButterMask
maskedInvf = np.fft.ifft2(np.fft.ifftshift(lpFshift))
lpfImg = np.abs(maskedInvf)
hpFshift = fshift * hpButterMask
maskedInvf = np.fft.ifft2(np.fft.ifftshift(hpFshift))
hpfImg = np.abs(maskedInvf)
plt.subplot(221), plt.imshow(lpButterMask, 'gray'), plt.title('Butterworth LPF')
plt.subplot(222), plt.imshow(lpfImg, 'gray'), plt.title('LPF Image')
plt.subplot(223), plt.imshow(hpButterMask, 'gray'), plt.title('HPF Spectrum')
plt.subplot(224), plt.imshow(hpfImg, 'gray'), plt.title('HPF Image')
plt.show() 
很顯然，濾波器的選擇也是很重要的，這里用到的是 Butterworth 濾波器，有興趣的可以自己實現一下:P
總結
傅里葉變換真是又偉大、又深奧、又方便！至少對于一個微積分已經忘記差不多、三角函數公式都要搜索半天才能回憶起來的“文科生”來說，從頭學習一遍簡直是又有虐腦又有驚喜！這篇文章可能更加偏重于記錄我自己的消化過程，如果想要更加細致、深入深刻以及深入淺出地介紹，請參考下方參考鏈接。
這兩天被虐腦的感覺真是酸爽，打算接續這一篇下去寫一個專題，就叫#BlowYourMind3000#，哈哈哈，下篇預告《貝葉斯》。
參考
-
Schyns, P. G., & Oliva, A. (1994). From blobs to boundary edges: evidence for time- and spatial-scale-dpendent scene recognition.Psychological Science, 5(5), 195-200.