圖像的正交變換在數字圖像的處理與分析中起着很重要的作用，被廣泛應用於圖像增強、去噪、壓縮編碼等眾多領域。本文手工實現了二維離散傅里葉變換和二維離散余弦變換算法，並在多個圖像樣本上進行測試，以探究二者的變換效果。
1. 傅里葉變換
實驗原理
對一幅圖像進行離散傅里葉變換（DFT），可以得到圖像信號的傅里葉頻譜。二維 DFT 的變換及逆變換公式如下：
DFT 盡管解決了頻域離散化的問題，但運算量太大。從公式中可以看到，有兩個嵌套的求和符號，顯然直接計算的復雜度為 \(O(n^2)\) 。為了加快傅里葉變換的運算速度，后人提出快速傅里葉變換（FFT），即蝶形算法，將計算 DFT 的復雜度降低到了 \(O(n\log n)\)。
FFT 利用傅里葉變換的數學性質，采用分治的思想，將一個 \(N\) 點的 FFT，變成兩個 \(N/2\) 點的 FFT。以一維 FFT 為例，可以表示如下：
其中，\(G(k)\) 是 \(x(k)\) 的偶數點的 \(N/2\) 點的 FFT，\(H(k)\) 是 \(x(k)\) 的奇數點的 \(N/2\) 點的 FFT。
這樣，通過將原問題不斷分解為兩個一半規模的子問題，然后計算相應的蝶形運算單元，最終得以完成整個 FFT。
算法步驟
本次實驗中，一維 FFT 采用遞歸實現，且僅支持長度為 2 的整數冪的情況。
算法步驟如下：
- 檢查圖像的尺寸，如果不是 2 的整數冪則直接退出。
- 對圖像的灰度值進行歸一化。
- 對圖像的每一行執行一維 FFT，並保存為中間結果。
- 對上一步結果中的每一列執行一維 FFT，返回變換結果。
- 將零頻分量移到頻譜中心，並求絕對值進行可視化。
- 對中心化后的結果進行對數變換，以改善視覺效果。
主要代碼
一維 FFT
def fft(x):
n = len(x)
if n == 2:
return [x[0] + x[1], x[0] - x[1]]
G = fft(x[::2])
H = fft(x[1::2])
W = np.exp(-2j * np.pi * np.arange(n//2) / n)
WH = W * H
X = np.concatenate([G + WH, G - WH])
return X
二維 FFT
def fft2(img):
h, w = img.shape
if ((h-1) & h) or ((w-1) & w):
print('Image size not a power of 2')
return img
img = normalize(img)
res = np.zeros([h, w], 'complex128')
for i in range(h):
res[i, :] = fft(img[i, :])
for j in range(w):
res[:, j] = fft(res[:, j])
return res
零頻分量中心化
def fftshift(img):
# swap the first and third quadrants, and the second and fourth quadrants
h, w = img.shape
h_mid, w_mid = h//2, w//2
res = np.zeros([h, w], 'complex128')
res[:h_mid, :w_mid] = img[h_mid:, w_mid:]
res[:h_mid, w_mid:] = img[h_mid:, :w_mid]
res[h_mid:, :w_mid] = img[:h_mid, w_mid:]
res[h_mid:, w_mid:] = img[:h_mid, :w_mid]
return res
運行結果
2. 余弦變換
實驗原理
當一個函數為偶函數時，其傅立葉變換的虛部為零，因而不需要計算，只計算余弦項變換，這就是余弦變換。離散余弦變換（DCT）的變換核為實數的余弦函數，因而計算速度比變換核為指數的 DFT 要快得多。
一維離散余弦變換與離散傅里葉變換具有相似性，對離散傅里葉變換進行下式的修改：
式中
由上式可見，\(\sum\limits_{x=0}^{2M-1}f_e(x)e^{\frac{-j2ux\pi}{2M}}\) 是 \(2M\) 個點的傅里葉變換，因此在做離散余弦變換時，可將其拓展為 \(2M\) 個點，然后對其做離散傅里葉變換，取傅里葉變換的實部就是所要的離散余弦變換。
算法步驟
基於上述原理，二維 DCT 的實現重用了上文中的一維 FFT 函數，並根據公式做了一些修改。
算法步驟如下：
- 檢查圖像的尺寸，如果不是 2 的整數冪則直接退出。
- 對圖像的灰度值進行歸一化。
- 對圖像的每一行進行延拓，執行一維 FFT 后取實部，乘以公式中的系數，並保存為中間結果。
- 對上一步結果中的每一列進行延拓，執行一維 FFT 后取實部，乘以公式中的系數，返回變換結果。
- 對結果求絕對值，並進行對數變換，以改善視覺效果。
主要代碼
二維 DCT
def dct2(img):
h, w = img.shape
if ((h-1) & h) or ((w-1) & w):
print('Image size not a power of 2')
return img
img = normalize(img)
res = np.zeros([h, w], 'complex128')
for i in range(h):
res[i, :] = fft(np.concatenate([img[i, :], np.zeros(w)]))[:w]
res[i, :] = np.real(res[i, :]) * np.sqrt(2 / w)
res[i, 0] /= np.sqrt(2)
for j in range(w):
res[:, j] = fft(np.concatenate([res[:, j], np.zeros(h)]))[:h]
res[:, j] = np.real(res[:, j]) * np.sqrt(2 / h)
res[0, j] /= np.sqrt(2)
return res
運行結果
完整源碼請見 GitHub 倉庫