Python 實現影象快速傅立葉變換和離散餘弦變換
阿新 • • 發佈：2020-07-19
影象的正交變換在數字影象的處理與分析中起著很重要的作用，被廣泛應用於影象增強、去噪、壓縮編碼等眾多領域。本文手工實現了**二維離散傅立葉變換**和**二維離散餘弦變換**演算法，並在多個影象樣本上進行測試，以探究二者的變換效果。
## 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 的整數冪的情況。
演算法步驟如下：
1. 檢查影象的尺寸，如果不是 2 的整數冪則直接退出。
2. 對影象的灰度值進行歸一化。
3. 對影象的每一行執行一維 FFT，並儲存為中間結果。
4. 對上一步結果中的每一列執行一維 FFT，返回變換結果。
5. 將零頻分量移到頻譜中心，並求絕對值進行視覺化。
6. 對中心化後的結果進行對數變換，以改善視覺效果。
### 主要程式碼
**一維 FFT**
```python
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**
```python
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
```
**零頻分量中心化**
```python
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 函式，並根據公式做了一些修改。
演算法步驟如下：
1. 檢查影象的尺寸，如果不是 2 的整數冪則直接退出。
2. 對影象的灰度值進行歸一化。
3. 對影象的每一行進行延拓，執行一維 FFT 後取實部，乘以公式中的係數，並儲存為中間結果。
4. 對上一步結果中的每一列進行延拓，執行一維 FFT 後取實部，乘以公式中的係數，返回變換結果。
5. 對結果求絕對值，並進行對數變換，以改善視覺效果。
### 主要程式碼
**二維 DCT**
```python
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
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 的整數冪的情況。
演算法步驟如下：
1. 檢查影象的尺寸，如果不是 2 的整數冪則直接退出。
2. 對影象的灰度值進行歸一化。
3. 對影象的每一行執行一維 FFT，並儲存為中間結果。
4. 對上一步結果中的每一列執行一維 FFT，返回變換結果。
5. 將零頻分量移到頻譜中心，並求絕對值進行視覺化。
6. 對中心化後的結果進行對數變換，以改善視覺效果。
### 主要程式碼
**一維 FFT**
```python
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**
```python
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
```
**零頻分量中心化**
```python
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 函式，並根據公式做了一些修改。
演算法步驟如下：
1. 檢查影象的尺寸，如果不是 2 的整數冪則直接退出。
2. 對影象的灰度值進行歸一化。
3. 對影象的每一行進行延拓，執行一維 FFT 後取實部，乘以公式中的係數，並儲存為中間結果。
4. 對上一步結果中的每一列進行延拓，執行一維 FFT 後取實部，乘以公式中的係數，返回變換結果。
5. 對結果求絕對值，並進行對數變換，以改善視覺效果。
### 主要程式碼
**二維 DCT**
```python
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