在之前的文章中，我們介紹了傅立葉變換的本質和很多基本性質，現在，該聊聊代碼實現的問題了。

為了方便起見，本文采用的編程語言是 Python3，矩陣處理用 numpy，圖像處理則使用 OpenCV3。

離散傅立葉變換

首先，回憶一下離散傅立葉變換的公式：
\[ \begin{eqnarray} F(u, v)&=&\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x, y)e^{-j2\pi ux/M}e^{-j2\pi vy/N} \&=&\frac{1}{MN}\sum_{y=0}^{N-1}\lbrace \sum_{x=0}^{M-1}f(x, y)e^{-j2\pi ux/M}\rbrace e^{-j2\pi vy/N} \end{eqnarray} \]

根據這種可分性，我們可以將二維的計算分為兩個一維進行。

現在，考慮如何計算 \(F(u,y)=\sum_{x=0}^{M-1}f(x, y)e^{-j2\pi ux/M}\)，這個式子中的 \(y\) 可以當作是常數，所以這其實是關於 \(x\) 的一維運算。

根據這個式子，可以得到：
\[ F(0,y)=\sum_{x=0}^{M-1}f(x,y)e^{-j2\pi 0 x} \F(1,y)=\sum_{x=0}^{M-1}f(x,y)e^{-j2\pi x/M} \\dots \F(M-1,y)=\sum_{x=0}^{M-1}f(x,y)e^{-j2\pi (M-1)x/M} \]

當然，由於圖片是二維的，所以 \(f(x,y)\) 對應的向量實際上應該是：
\[ \begin{bmatrix} f(0,y_1) & f(0,y_2) & \dots & f(0,y_N) \f(1,y_1) & f(1,y_2) & \dots & f(1,y_N) \\dots \f(M-1, y_1) & f(M-1,y_2) & \dots & f(M-1,y_N) \end{bmatrix} \]
同理，得到的 \(F(u,y)\) 也是一個二維矩陣。

現在，我們還是先考慮怎麽實現這個一維的計算。

首先，需要先把 \(W_M\) 這個矩陣表示出來。註意到，這個矩陣實際上可以由 \(\begin{bmatrix} W_M^0 \\ W_M^1 \\ \vdots \\ W_M^{M-1} \end{bmatrix}\) \(\times\) \(\begin{bmatrix} W_M^0 & W_M^1 & \dots & W_M^{M-1} \end{bmatrix}\) 得到。借助 numpy 強大的矩陣處理能力，可以很方便的計算出這個矩陣。示例如下：

def dftmtx(M):
    n = np.asmatrix(np.arange(M))
    return np.exp((-2j * np.pi / M) * n.transpose() * n)

np.asmatrix 是把 M 維的向量變成 1 \(\times\) M 的矩陣的格式，因為只有矩陣才有 transpose() 操作。np.exp 會把 \(exp\) 函數作用到矩陣的每個元素中。

得到這個矩陣後，最關鍵的一步其實就做完了，我們可以用這個矩陣計算出 \(F(u,y)\)：

# input表示輸入圖像，M是圖像的高
M = input.shape[0]
F = dftmtx(M) * input

得到 \(F(u,y)\) 後，剩下的是要對 y 這一維進行同樣的操作：\(F(u,v)=\frac{1}{MN}\sum_{y=0}^{N-1}F(u,y) e^{-j2\pi vy/N}\)。同樣地，我們需要計算一個 \(W_N\) 的矩陣。幸運的是，這個矩陣的計算方法和之前的 \(W_M\) 一模一樣，這樣一來，我們已經可以得到完整的計算方法了：

# 傅立葉變換函數
def dft2d(input):
    M, N = input.shape[0], input.shape[1]
    return dftmtx(M) * input * dftmtx(N) / (M * N)

接下來我們把頻譜圖打印出來。傅立葉頻譜圖是實部和虛部的平方和，需要註意的是，由於數值顯示的問題，我們需要將頻譜圖用 log 函數增強後，再標定到 [0, 255] 之間才能看清。代碼如下：

# 將像素值標定到[0，255]之間
def scale_intensity(image):
    min = image.min()
    max = image.max()
    image = (image - min) / (max - min) * 255.0
    return image

# 計算頻譜圖
def spectrogram(image):
    dft = dft2d(image)
    spec = np.sqrt(np.power(np.real(dft), 2) + np.power(np.imag(dft), 2))
    spec = np.log(0.5 + spec) * 10
    spec = scale_intensity(spec)
    return spec
 
image = cv2.imread("your_file.png", cv2.IMREAD_GRAYSCALE)
spec = spectrogram(image)
cv2.imwrite("spec.png", spec)

結果展示：

左圖是原圖，右圖是頻譜圖。如果仔細看的話，可以發現頻譜圖四個角有一些白色的點。這是因為圖片中低頻成分居多，而頻譜圖四個角代表的就是低頻分量（至於為什麽四個角是低頻，我也沒搞懂）。

實踐中，人們習慣於把低頻都聚集到圖片中心，這樣方便後續的操作。根據平移性質：
\[ F(u-\frac{M}{2},v-\frac{N}{2}) =f(x,y)(-1)^{x+y} \]
要把頻譜圖的低頻部分平移到中心，需要將整個頻譜圖平移 \((M/2, N/2)\) 個單位，也就是需要對原圖乘以 \((-1)^{x+y}\)。代碼如下：

def shift_image(image):
    M, N = image.shape[0], image.shape[1]
    for x in range(M):
        for y in range(N):
            image[x, y] *= np.power(-1, x + y)
    return image

  
image = cv2.imread("your_file.png", cv2.IMREAD_GRAYSCALE)
shift_image(image)
spec = spectrogram(image)
cv2.imwrite("shift_spec.png", spec)

結果展示：

離散傅立葉反變換

講完傅立葉變換後，反變換基本也得到了，唯一的區別是，這一次我們需要計算一個傅立葉反變換的矩陣。這個矩陣和之前計算的矩陣 \(W_M\) 的區別只在於符號，這裏就直接給出代碼了：

def idftmtx(M):
    n = np.asmatrix(np.arange(M))
    # 下面的符號是正的
    return np.exp((2j * np.pi / M) * n.transpose() * n)

反變換的代碼如下：

def idft2d(input):
    M, N = input.shape[0], input.shape[1]
    return idftmtx(M) * input * idftmtx(N)

把之前得到的傅立葉變換的結果，輸入 idft2d 函數後，再取實部既可以得到原圖：

image = cv2.imread("your_file.png", cv2.IMREAD_GRAYSCALE)
dft = dft2d(image)
idft = idft2d(dft)
cv2.imwrite("idft.png", np.real(idft))

結果如下：

這個反變換的結果和原圖是略有差別的，因為傅立葉變換時舍棄了很多高頻成分。不過，由於圖片中高頻成分本身就比較少，所以這點差別可以忽略不計。

圖像中的傅立葉變換（三）