Python 實現圖像快速傅里葉變換和離散餘弦變換

圖像的正交變換在數字圖像的處理與分析中起着很重要的做用，被普遍應用於圖像加強、去噪、壓縮編碼等衆多領域。本文手工實現了二維離散傅里葉變換和二維離散餘弦變換算法，並在多個圖像樣本上進行測試，以探究兩者的變換效果。python

1. 傅里葉變換

實驗原理

對一幅圖像進行離散傅里葉變換（DFT），能夠獲得圖像信號的傅里葉頻譜。二維 DFT 的變換及逆變換公式以下：git

DFT 儘管解決了頻域離散化的問題，但運算量太大。從公式中能夠看到，有兩個嵌套的求和符號，顯然直接計算的複雜度爲 \(O(n^2)\) 。爲了加快傅里葉變換的運算速度，後人提出快速傅里葉變換（FFT），即蝶形算法，將計算 DFT 的複雜度下降到了 \(O(n\log n)\)。github

FFT 利用傅里葉變換的數學性質，採用分治的思想，將一個 \(N\) 點的 FFT，變成兩個 \(N/2\) 點的 FFT。以一維 FFT 爲例，能夠表示以下：算法

其中，\(G(k)\) 是 \(x(k)\) 的偶數點的 \(N/2\) 點的 FFT，\(H(k)\) 是 \(x(k)\) 的奇數點的 \(N/2\) 點的 FFT。函數

這樣，經過將原問題不斷分解爲兩個一半規模的子問題，而後計算相應的蝶形運算單元，最終得以完成整個 FFT。測試

算法步驟

本次實驗中，一維 FFT 採用遞歸實現，且僅支持長度爲 2 的整數冪的狀況。編碼

算法步驟以下：spa

檢查圖像的尺寸，若是不是 2 的整數冪則直接退出。
對圖像的灰度值進行歸一化。
對圖像的每一行執行一維 FFT，並保存爲中間結果。
對上一步結果中的每一列執行一維 FFT，返回變換結果。
將零頻份量移到頻譜中心，並求絕對值進行可視化。
對中心化後的結果進行對數變換，以改善視覺效果。

主要代碼

一維 FFT3d

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

二維 FFTcode

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 倉庫