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

TimDyh發表於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 的整數冪的情況。

演算法步驟如下：

檢查影像的尺寸，如果不是 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

執行結果

lena_fft

baboon_fft

circle_fft

avicii_fft

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