IDCT的原理以及代码分析

代码版本： HM16.2

IDCT（逆DCT）变换是DCT变换的逆变换。在HEVC的解码端对图像的像素残差进行IDCT变换，再进行反量化，就可以得到最初的残差值， 原像素 = 预测值 + 残差值。这篇文章基于HM16.2版本的解码代码进行原理分析以及代码实现。

原理分析：
令： H为编码端的变换矩阵，X为像素值矩阵， Y为变换后得到的像素值矩阵。
变换公式大体为：。
n为变换矩阵的尺寸，在HEVC中存在4种大小的变换矩阵：4x4, 8x8, 16x16, 32 x32。

在解码端的IDCT公式为：。
从上式可以看出，解码端相对编码端的求解顺序是类似的，唯一额外的计算为求得变换矩阵的逆矩阵。
通过分析得知在HEVC的DCT变换矩阵为正交矩阵，正交矩阵的特性之一正交矩阵的逆矩阵等于它的转置矩阵。由此可以得到以下推倒公式：

由上式可见，在整个计算过程我们只需要知道变换矩阵H和像素值矩阵Y就可以求得矩阵X。

解码代码分析：（以4x4IDCT为例）
代码实现分两部分组成，第一步为前两个矩阵相乘，保存结果为中间变量；
第二步为中间变量与后一个矩阵相乘。
为了简化计算过程，代码采取了蝶形变换，蝶形变换之所以能够减少计算次数，它的核心思想为根据矩阵系数的特点，抽取共同或者类似部分系数，计算这些系数的结果，并由此推导出其它部分的结果。
如4x4变换矩阵：

可以看出矩阵A的第一第二列的1,3行 与 第三第四列的1,3行成对称；同理它们的2,4行成反对称。利用这一个特点，采取蝶形变换，可以大大减少计算次数。

代码：

/** 4x4 inverse transform implemented using partial butterfly structure (1D)
 *  \param src   input data (transform coefficients)
 *  \param dst   output data (residual)
 *  \param shift specifies right shift after 1D transform
 */
Void partialButterflyInverse4(TCoeff *src, TCoeff *dst, Int shift, Int line, const TCoeff outputMinimum, const TCoeff outputMaximum)
{
  Int j;
  TCoeff E[2],O[2];
  TCoeff add = (shift > 0) ? (1<<(shift-1)) : 0;

  for (j=0; j<line; j++)
  {
    /* Utilizing symmetry properties to the maximum to minimize the number of multiplications */
    //g_aiT4[TRANSFORM_INVERSE][n][m]为4x4变换矩阵n行m列的系数
    O[0] = g_aiT4[TRANSFORM_INVERSE][1][0]*src[line] + g_aiT4[TRANSFORM_INVERSE][3][0]*src[3*line];
    O[1] = g_aiT4[TRANSFORM_INVERSE][1][1]*src[line] + g_aiT4[TRANSFORM_INVERSE][3][1]*src[3*line];
    E[0] = g_aiT4[TRANSFORM_INVERSE][0][0]*src[0]    + g_aiT4[TRANSFORM_INVERSE][2][0]*src[2*line];
    E[1] = g_aiT4[TRANSFORM_INVERSE][0][1]*src[0]    + g_aiT4[TRANSFORM_INVERSE][2][1]*src[2*line];

    /* Combining even and odd terms at each hierarchy levels to calculate the final spatial domain vector */
    dst[0] = Clip3( outputMinimum, outputMaximum, (E[0] + O[0] + add)>>shift );
    dst[1] = Clip3( outputMinimum, outputMaximum, (E[1] + O[1] + add)>>shift );
    dst[2] = Clip3( outputMinimum, outputMaximum, (E[1] - O[1] + add)>>shift );
    dst[3] = Clip3( outputMinimum, outputMaximum, (E[0] - O[0] + add)>>shift );

    src   ++;
    dst += 4;
  }
}

上述函数被调用两次，计算过程与我们推导出来的公式略有不同，但是结果一样，这样做是对原有公式的一种优化，代码得到了重用。
正常的计算过程应该如下表示：


该每次调用该函数的计算过程可用公式等价如下表示（*为该函数的矩阵乘法规则）：
第一次调用：；

第二次调用：；
由上式可以看出，虽然求解过程不一样，但是最终结果一致。