1. 程式人生 > >Python 實現影象快速傅立葉變換和離散餘弦變換

Python 實現影象快速傅立葉變換和離散餘弦變換

影象的正交變換在數字影象的處理與分析中起著很重要的作用,被廣泛應用於影象增強、去噪、壓縮編碼等眾多領域。本文手工實現了**二維離散傅立葉變換**和**二維離散餘弦變換**演算法,並在多個影象樣本上進行測試,以探究二者的變換效果。 ## 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