DCT变换

最新推荐文章于 2025-02-19 16:19:49 发布

涵小呆

最新推荐文章于 2025-02-19 16:19:49 发布

阅读量3.2w

点赞数 41

分类专栏：视频编码文章标签： DCT变换

本文链接：https://blog.csdn.net/BigDream123/article/details/101426393

版权

视频编码专栏收录该内容

32 篇文章

订阅专栏

DCT变换

1.DCT变换公式及其性质

傅里叶变换表明，任何信号都能表示为多个不同振幅和频率的正弦或者余弦信号的叠加。如果采用的是余弦函数，则信号分解过程称为余弦变换；若输入信号是离散的，则称之为离散余弦变换（Discrete Cosine Transform，DCT）。数学上共存在8种类型的DCT。公式如下：

其中，在图像视频领域中，最常用的是DCT-Ⅱ，平常说的DCT一般指的是DTC-Ⅱ。DTC-Ⅲ是 DTC-Ⅱ的反变换，一般说的反DCT指的就是DCT-Ⅲ。

以DCT-Ⅱ公式为例，公式为：

$X(k) = \sqrt{\frac{2}{N}}\epsilon (k)\sum x(n)cos(\frac{k(2n+1 )\pi}{2N}),k=0,1,...N-1$

由公式可以看出，X(k)就是在频率k下，求序列上各点对应余弦项的累加和。当k=0时，各点的余弦项均为1，此时得到的结果X(0)正比于x(n)的和，被称为直流分量。当k>0时，得到个频率下的X(k)被称为交流分量。

DCT变换的C++实现可以参考https://blog.csdn.net/BigDream123/article/details/104395587

2.DCT性质

离散余弦变换具有两点性质，第一个为可分离性，第二个为能量集中性。

可分离性

以DCT-Ⅱ为例，二维DCT变换表示为：

DCT的可分离性就是可以将二维DCT拆分为两个方向的一维DCT。

能量集中

图像编码中，图像以矩阵形式存储，是一个二维数组，对图像进行二维DCT，可以将时域图像转为频域能量分布图，实现能量集中从而去除空间冗余。

例存在矩阵：

10 10 10 10
10 10 10 10
10 10 10 10
10 10 10 10

对其进行DCT变换后，变为：

40.000000 0.000001 -0.000002 0.000002
0.000001 0.000000 -0.000000 0.000000
-0.000002 -0.000000 0.000000 -0.000000
0.000002 0.000000 -0.000000 0.000000

可见，大多数能量集中在了左上角低频处。

图像中，低频分量可以看作是基本图像，高频分量为边缘轮廓细节信息（也可能是噪声）。绝大多数能量都包含了大量平坦区域，因此大部分能量集中于低频区域（左上角），从而可以达到去除空间冗余的目的。（其中，DCT变换尺寸越大，DCT去相关性能越好，但是DCT的计算复杂度会随之增大）。

3.DCT基函数和基图像

图像信息的基函数表示方法中，一般将图像I(x,y）表示成一些特征或者基函数线性叠加的形式：

$I(x,y) = \sum a_{i}(x,y)s_{i}$

其中，基函数 $a_{i}$ (x,y)可以看作是图像的基本组成成分， $s_{i}$ 表示系数，可看成是每个基函数对应的幅度，是原始图像在基函数张成的高维空间的投影系数。

将上式可以表示为矩阵形式：

x=As

这可以看作是图像的重构过程，通过一组被适当加权的基函数求和重构原图像，即从变换域到空域变换的映射。

考虑到基函数ai(x,y)构成一个可逆线性系统，矩阵A是一个方阵，这样可以对系统求逆得到：

$s_{i} = \sum w_{i}(x,y)I(x,y)$

上式是图像的变换过程：将图像从空域映射到变换域用基函数表示，求原图像在各基函数所占的量si。

式中Wi即为变换基。上式可以用矩阵表示：

s = Wx

当变换系数s的个数和原图像像素x的个数相同时，这样在变换前后的维度相同，保证变换过程中没有引入新的信息，也未破坏原有的信息，因此通过变换系数s可以无损的恢复原图像x。

每一种变换都有其对应的基函数，各变换的不同之处，也仅是基函数不同而已。因此在分析不同变换的性能时，经常会对其基函数进行分析。

以图像变换为例，定义原图像f(x,y)，通过g(x,y,u,v)变换可以得到变换域函数F(u,v)。

$F(u,v) = \sum_{x}^{M-1} \sum_{y}^{M-1} f(x,y)g(x,y,u,v)$

$f(x,y) = \sum_{u}^{M-1} \sum_{v}^{M-1} F(u,v)h(x,y,u,v)$

这里g(x,y,u,v)为正变换核，反变换核h(x,y,u,v)就是基函数。f(x,y)就等于h(x,y,u,v)以F(u,v)为权重加权求和.

下面以DCT-Ⅱ为例，二维变换如下：

反变换为：

基函数如下表示：

$B(m,n,k,l) = C(k)C(l)cos(\frac{(2m+1)k\pi}{2N})cos(\frac{(2n+1)l\pi}{2N})$

从反变换公式看出，x(m,n)是由NxN个频率分量组成的，每个频率分量与一个变换系数相对应。对于每个变换系数，遍历(m,n)所有值时，构成了一幅基图像。将变换系数看作权重，对基函数加权求和，就可以得到时域信号x(m,n)。因此基函数是构成时域信号的基本单元。

N=4时基函数图像如下：（DCT变换基图像的MATLAB绘制代码https://download.csdn.net/download/BigDream123/12022161）

当k=0，l=0时，表示水平垂直分量的频率都等于0，即左上角的基图像，此时基函数B(m,n)=1/4，图像平坦，没有变化；右上角表示水平频率k最大而垂直频率为0的基图像，图像在水平方向上发生了连续变化。每个(k,l)对应的基函数图像，就是该(k,l)频率下的构成原始图像的基本单元。

简而言之，对于离散余弦变换可以解释为：将任一4x4像素块表示为图所示的16个基图像的加权和，其权值即为对应位置的DCT系数。

4.整数DCT变换

自H.264中首次采用了整数DCT变换，下面以4x4DCT为例推导H.265中使用的4x4整数DCT变换。二维4x4DCT变换公式如下：

其中，

将上式进行变形，可以得到

将括号中的部分记为Z(m,l)，则上式可以分解为一下二式：

可以看出，二维DCT变换可以分解为两个一维DCT，即可以先对像素块的行（或列）做一维DCT，再对像素块的列（或行）做一维DCT。现将一维DCT写成如下相乘的形式：

其中，X表示原始像素块，Y表示变换后的DCT系数矩阵，变换矩阵A由如下定义：

将A写成矩阵形式：

令

则有

对上式的a,b,c同时乘以128并整数化（四舍五入），得

为保持正交性，即满足需对b，c进行微调

最终的变换矩阵为

其对应的一维修正矩阵E较为简单，其所有元素为1/128，即

H.265/HEVC中还使用了更大尺寸的整数DCT：8x8，16x16，32x32这三种整数DCT推导方法与4x4基本相同，唯一的区别在于矩阵元素整数化时放大倍数不同。对于NxN的变换块，放大倍数为 $64\sqrt{N}$

5.蝶形变换

蝶形变换的思想就是提取矩阵的相关部分，定义中间变量，减少运算次数。

以4x4的一维DCT变换为例：

x

分析结果矩阵的第一列：

Y00 = 64*X00 + 64*X10 + 64*X20 + 64*X30

Y10 = 83*X00 + 36*X10 - 36*X20 - 83*X30

Y20 = 64*X00 - 64*X10 - 64*X20 + 64*X30

Y30 = 36*X00 - 83*X10 + 83*X20 - 36*X30

可以的定义四个中间变量

T1 = X00 + X30

T2 = X00 - X30

T3 = X10 + X20

T4 = X10 - X20

从而可以将第一列的结果转换为

Y00 = 64*X00 + 64*X10 + 64*X20 + 64*X30 = 64 * T1 + 64 * T3

Y10 = 83*X00 + 36*X10 - 36*X20 - 83*X30 = 83 * T2 + 36 * T4

Y20 = 64*X00 - 64*X10 - 64*X20 + 64*X30 = 64 * T1 - 64 * T3

Y30 = 36*X00 - 83*X10 + 83*X20 - 36*X30 = 36 * T2 - 83 * T4

通过定义临时变量，减少了重复的计算，以此加快计算速度，这就是蝶形变换。

VTM中的DCT变换就是使用蝶形快速算法，以DCT4x4为例：

/** 4x4 forward transform implemented using partial butterfly structure (1D)
*  \param src   input data (residual)
*  \param dst   output data (transform coefficients)
*  \param shift specifies right shift after 1D transform
*  \param line
*/
void fastForwardDCT2_B4(const TCoeff *src, TCoeff *dst, int shift, int line, int iSkipLine, int iSkipLine2)
{
  int j;
  TCoeff E[2], O[2];
  TCoeff add = (shift > 0) ? (1 << (shift - 1)) : 0;

  const TMatrixCoeff *iT = g_trCoreDCT2P4[TRANSFORM_FORWARD][0];

  TCoeff *pCoef = dst;
  const int  reducedLine = line - iSkipLine;
  for (j = 0; j<reducedLine; j++)
  {
    /* E and O */
    E[0] = src[0] + src[3];
    O[0] = src[0] - src[3];
    E[1] = src[1] + src[2];
    O[1] = src[1] - src[2];

    dst[0] = (iT[0] * E[0] + iT[1] * E[1] + add) >> shift;
    dst[2 * line] = (iT[8] * E[0] + iT[9] * E[1] + add) >> shift;
    dst[line] = (iT[4] * O[0] + iT[5] * O[1] + add) >> shift;
    dst[3 * line] = (iT[12] * O[0] + iT[13] * O[1] + add) >> shift;

    src += 4;
    dst++;
  }
  if (iSkipLine)
  {
    dst = pCoef + reducedLine;
    for (j = 0; j<4; j++)
    {
      memset(dst, 0, sizeof(TCoeff)*iSkipLine);
      dst += line;
    }
  }
}

其中，g_trCoreDCT2P4数组定义如下

const TMatrixCoeff g_trCoreDCT2P4[TRANSFORM_NUMBER_OF_DIRECTIONS][4][4] =
{
  DEFINE_DCT2_P4_MATRIX(64,    83,    36),
  DEFINE_DCT2_P4_MATRIX(64,    83,    36)
};
#define DEFINE_DCT2_P4_MATRIX(a,b,c) \
{ \
  { a,  a,  a,  a}, \
  { b,  c, -c, -b}, \
  { a, -a, -a,  a}, \
  { c, -b,  b, -c}  \
}