DCT 快速演算法實作：一維/二維正逆變換、查找表優化與多執行緒架構

離散餘弦變換 (DCT) 是 JPEG、MPEG、H.264 等影像/視訊壓縮核心。直接依定義計算二維 DCT 複雜度 O(N^4)，對 256×256 圖像需數十分鐘不可接受。本文實作實數域快速 DCT/IDCT (複雜度 O(N^2 log N))，採用查找表消除重複三角函數計算、位反序重排、多執行緒非阻塞架構，256×256 灰階圖 DCT/IDCT 僅耗 94ms (vs 普通 575 秒)，效能提升 **6000 倍以上**。

演算法理論基礎

參考《DCT 快速新演算法及濾波器結構研究與子波變換域影像去噪研究》華南理工大學博士論文。核心思想：

• 利用餘弦函數對稱性與遞迴分解，將 DCT 分解為「前向追底」 (蝶形運算疊加) + 「後向回根」 (係數還原) 兩階段。
• 全程於實數域運算，避免 FFT 引入複數運算的 4 倍乘法/2 倍加法開銷。
• 係數 cos((2k+1)π/4N) 僅依賴 N，預先計算存入查找表 C[]。

查找表初始化

double* C = nullptr;
bool bHasInit = false;

void initDCTParam(int deg) { // deg = log2(N)
    if (bHasInit) return;
    int N = 1 << deg;
    if (C) delete[] C;
    C = new double[N];
    
    int half = N >> 1;
    for (int i = 0; i < half; ++i)
        C[N - i - 1] = 2*i + 1; // 基礎奇數序列
    
    for (int g = 0; g < deg - 1; ++g) {
        int endStart = 1 << (deg - 1 - g);
        int len = endStart >> 1;
        int factor = 1 << (g + 1);
        for (int j = 0; j < len; ++j)
            C[endStart - j - 1] = factor * C[N - j - 1];
    }
    
    for (int i = 1; i < N; ++i) // C[0] 未使用
        C[i] = 2.0 * cos(C[i] * M_PI / (N << 1));
    
    bHasInit = true;
}

一維 DCT 正向變換

前向追底

void dct_forward(double* f, int deg) {
    int N = 1 << deg;
    for (int d = 0; d < deg; ++d) {
        int wings = 1 << d;
        int wingLen = N >> d;
        int halfWing = wingLen >> 1;
        for (int w = 0; w < wings; ++w) {
            for (int hw = 0; hw < halfWing; ++hw) {
                double t1 = f[w*wingLen + hw];
                double t2 = f[(w+1)*wingLen - 1 - hw];
                if (w & 1) swap(t1, t2); // 奇數翼交換
                f[w*wingLen + hw] = t1 + t2;
                f[(w+1)*wingLen - 1 - hw] = (t1 - t2) * C[wingLen - 1 - hw];
            }
        }
    }
}

後向回根

void dct_backward(double* f, int deg) {
    int N = 1 << deg;
    for (int d = deg - 1; d >= 0; --d) {
        int wings = 1 << d;
        int wingLen = 1 << (deg - d);
        int halfWing = wingLen >> 1;
        for (int w = 0; w < wings; ++w) {
            for (int hw = 0; hw < halfWing; ++hw) {
                if (hw == 0) {
                    f[halfWing + w*wingLen] *= 0.5;
                } else {
                    int r1 = bitrev(hw, deg - d - 1);
                    int r2 = bitrev(hw - 1, deg - d - 1);
                    f[halfWing + w*wingLen + r1] -= f[halfWing + w*wingLen + r2];
                }
            }
        }
    }
}

位反序重排

int bitrev(int x, int bits) {
    int rev = 0;
    for (int i = 0; i < bits; ++i) {
        rev = (rev << 1) | (x & 1);
        x >>= 1;
    }
    return rev;
}

void fbitrev(double* f, int deg) {
    int N = 1 << deg;
    double* tmp = new double[N];
    for (int i = 0; i < N; ++i) tmp[bitrev(i, deg)] = f[i];
    memcpy(f, tmp, N * sizeof(double));
    delete[] tmp;
}

完整一維 DCT (含縮放)

void fdct_1D(double* f, int deg) {
    initDCTParam(deg);
    dct_forward(f, deg);
    dct_backward(f, deg);
    fbitrev(f, deg);
    f[0] *= M_SQRT1_2; // sqrt(1/2)
    
    int N = 1 << deg;
    double scale = sqrt(2.0 / N);
    for (int i = 0; i < N; ++i) f[i] *= scale;
}

二維 DCT：列優先 + 行優先

double* temp_2D = nullptr;

void init2D_Param(int rows, int cols) {
    if (temp_2D) delete[] temp_2D;
    temp_2D = new double[max(rows, cols)];
}

void fdct_2D(double* f, int degRow, int degCol) {
    int rows = 1 << degRow, cols = 1 << degCol;
    init2D_Param(rows, cols);
    double scale = 2.0 / sqrt(rows * cols);
    
    // 1. 對每列做 1D DCT
    for (int r = 0; r < rows; ++r)
        fdct_1D_no_param(f + r*cols, degCol);
    
    // 2. 對每行做 1D DCT (需轉置暫存)
    for (int c = 0; c < cols; ++c) {
        for (int r = 0; r < rows; ++r)
            temp_2D[r] = f[r*cols + c];
        fdct_1D_no_param(temp_2D, degRow);
        for (int r = 0; r < rows; ++r)
            f[r*cols + c] = temp_2D[r] * scale;
    }
    bHasInit = false; // 允許下次不同尺寸重新初始化
}

IDCT 逆變換：對稱流程

查找表係數改為倒數：C[i] = 1.0 / (2.0 * cos(...))。前向追底改為疊加差分，後向回根改為蝶形平均。完整 fidct_1D、fidct_2D 程式碼見附件。

多執行緒非阻塞架構

避免 UI 凍結，將 DCT/IDCT 放入背景執行緒：

struct RUNINFO {
    int row, col;       // log2 尺寸
    double* data;       // 輸入輸出緩衝 (呼叫端配置)
    bool finished;
    DWORD startTick, endTick;
};

UINT ThreadProc_DCT(LPVOID p) {
    RUNINFO* info = (RUNINFO*)p;
    info->startTick = GetTickCount();
    fdct_2D(info->data, GetTwoIndex(info->row), GetTwoIndex(info->col));
    info->endTick = GetTickCount();
    info->finished = true;
    return 0;
}

// UI 執行緒輪詢 info->finished 或用 PostMessage 通知完成

效能實測 (WinXP, VC6, 256×256 灰階圖)

圖 1-3 為程式介面截圖：普通 DCT 耗時、快速 DCT 耗時、快速 IDCT 耗時。

應用場景與擴充

• JPEG 編碼器核心模組 (量化表乘除、Zigzag 掃描、霍夫曼編碼另行實作)
• 影像特徵提取 (DCT 係數統計、低頻能量集中特性)
• 視訊編碼器區塊內預測殘差變換
• 彩色圖像：YCbCr 三通道並行處理 (可擴充 RUNINFO 含 dataY/dataCb/dataCr)

總結

實數域快速 DCT/IDCT 結合查找表與位反序，將 O(N^4) 降為 O(N^2 log N)，配合多執行緒非阻塞架構，滿足即時影像處理效能需求。完整原始碼含一維/二維正逆變換、查找表管理、多執行緒包裝、MFC 測試專案 (載入 BMP、顯示 DCT 係數圖、IDCT 還原驗證)，請下載附件。