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
```
### 執行結果
![lena_fft](https://ipichub.oss-cn-hangzhou.aliyuncs.com/2020-07-19-023959.jpg)
![baboon_fft](https://ipichub.oss-cn-hangzhou.aliyuncs.com/2020-07-19-024003.jpg)
![circle_fft](https://ipichub.oss-cn-hangzhou.aliyuncs.com/2020-07-19-024005.jpg)
![avicii_fft](https://ipichub.oss-cn-hangzhou.aliyuncs.com/2020-07-19-024009.jpg)
## 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
```
### 執行結果
![lena_dct](https://ipichub.oss-cn-hangzhou.aliyuncs.com/2020-07-19-024010.jpg)
![baboon_dct](https://ipichub.oss-cn-hangzhou.aliyuncs.com/2020-07-19-024017.jpg)
![circle_dct](https://ipichub.oss-cn-hangzhou.aliyuncs.com/2020-07-19-024019.jpg)
![avicii_dct](https://ipichub.oss-cn-hangzhou.aliyuncs.com/2020-07-19-024020.jpg)
**完整原始碼請見 GitHub