一维 IDCT 的实现是如何工作的？

Question

我有一个逆离散余弦变换的实现，我试图弄清楚他们是如何得到这个代码的。到目前为止，我发现这可能是 DCT 而不是 DFT（离散傅里叶变换）的 Cooley-Tukey radix-2 Decimation-in-time 的优化实现。

但是，我仍然不知道每个阶段到底发生了什么。我认为 Wx 常数可能是旋转因子。

任何人都可以提供解释参考，或者对此代码提供一些解释吗？

//Twiddle factors
#define W1 2841 /* 2048*sqrt(2)*cos(1*pi/16) */
#define W2 2676 /* 2048*sqrt(2)*cos(2*pi/16) */
#define W3 2408 /* 2048*sqrt(2)*cos(3*pi/16) */
#define W5 1609 /* 2048*sqrt(2)*cos(5*pi/16) */
#define W6 1108 /* 2048*sqrt(2)*cos(6*pi/16) */
#define W7 565 /* 2048*sqrt(2)*cos(7*pi/16) */

//Discrete Cosine Transform on a row of 8 DCT coefficients.
NJ_INLINE void njRowIDCT(int* blk) {
    int x0, x1, x2, x3, x4, x5, x6, x7, x8;
    int t;
    if (!((x1 = blk[4] << 11)
        | (x2 = blk[6])
        | (x3 = blk[2])
        | (x4 = blk[1])
        | (x5 = blk[7])
        | (x6 = blk[5])
        | (x7 = blk[3])))
    {
        blk[0] = blk[1] = blk[2] = blk[3] = blk[4] = blk[5] = blk[6] = blk[7] = blk[0] << 3;
        return;
    }
    x0 = (blk[0] << 11) + 128;  //For rounding at fourth stage

    //First stage
    /*What exactly are we doing here? Do the x values have a meaning?*/
    x8 = W7 * (x4 + x5);
    x4 = x8 + (W1 - W7) * x4;
    x5 = x8 - (W1 + W7) * x5;
    x8 = W3 * (x6 + x7);
    x6 = x8 - (W3 - W5) * x6;
    x7 = x8 - (W3 + W5) * x7;

    //Second stage
    x8 = x0 + x1;
    x0 -= x1;
    x1 = W6 * (x3 + x2);
    x2 = x1 - (W2 + W6) * x2;
    x3 = x1 + (W2 - W6) * x3;
    x1 = x4 + x6;
    x4 -= x6;
    x6 = x5 + x7;
    x5 -= x7;

    //Third stage
    x7 = x8 + x3;
    x8 -= x3;
    x3 = x0 + x2;
    x0 -= x2;
    x2 = (181 * (x4 + x5) + 128) >> 8;
    x4 = (181 * (x4 - x5) + 128) >> 8;

    //Fourth stage
    blk[0] = (x7 + x1) >> 8;  //bit shift is to emulate 8 bit fixed point precision
    blk[1] = (x3 + x2) >> 8;
    blk[2] = (x0 + x4) >> 8;
    blk[3] = (x8 + x6) >> 8;
    blk[4] = (x8 - x6) >> 8;
    blk[5] = (x0 - x4) >> 8;
    blk[6] = (x3 - x2) >> 8;
    blk[7] = (x7 - x1) >> 8;

}

Answer 1

我不是 DCT 方面的专家，但我曾经写过一些 FFT 实现，所以我将尝试回答这个问题。请对以下内容持保留态度。

void njRowIDCT(int* blk)

您正确地说该算法似乎是 8 长度 Radix-2 DCT，使用 24:8 精度的定点算术。我猜测精度是因为最后一级右移 8 以获得所需的结果（以及说明性的评论；）

因为它的长度是 8，所以它的幂是 3 (2^3 = 8)，这意味着 DCT 有 3 个阶段。到目前为止，这与 FFT 非常相似。 “第四阶段”似乎只是定点运算后恢复原始精度的缩放。

据我所知，输入阶段是从输入数组 blk 到局部变量 x0-x7 的 位反转。 x8 似乎是一个临时变量。抱歉，我无法比这更具描述性。

位反转阶段

更新

查看科学家和工程师的 DSP。它对信号处理主题给出了清晰而准确的解释。本章介绍 DCT（请跳至第 497 页）。

Wn（旋转因子）对应于本章中的基函数，尽管请注意这是描述 8x8 (2D) DCT。

关于我提到的3个阶段，与8点FFT的描述进行比较：

FFT 对位反转输入数组（本质上是复杂的乘法加法）执行蝶式运算，沿途将一条路径乘以 Wn 或旋转因子。 FFT 是分阶段执行的。我还没有弄清楚你的 DCT 代码在做什么，但将其分解成这样的图表可能会有所帮助。

那个或知道自己在说什么的人站出来;-)

我将重新访问此页面并在破译更多代码时进行编辑

Answer 2

DFT 和 DCT 都只是线性变换，可以表示为单个复数矩阵乘法（有时针对严格实数输入进行修剪）。因此，您可以组合上述方程来获得每个最终项的公式，该公式最终应等于线性变换矩阵的一行（忽略舍入问题）。然后看看上面的代码序列如何手动执行公共子表达式优化或行计算之间和/或内部的重构。

Answer 3

此代码基于 Loeffler 等人的经典论文 Practical Fast 1-D DCT Algorithms With 11 Multiplications。这实际上是从右向左遍历纸张的图 7。