xilinx暑期学校学习笔记(二) 离散余弦变换DCT压缩算法

今天的实验里面有个涉及到DCT的，虽然不用我们具体去写DCT的代码，不过还是来看一下DCT到底是个什么玩意吧。

DCT基本算法原理

DCT： 离散余弦变换(Discrete Cosine Transform)，

一维DCT变换

基本定义
如果 { X ( m ) ∣ m = 0 , 1 , . . . , N − 1 } \{X(m)|m=0,1,...,N-1\} {X(m)∣m=0,1,...,N−1}是有限带宽信号x(t)的采样序列，则一维离散余弦变换为：
$$Y(u)=C(u)\sqrt{\frac{2}{N}}\sum _{m=0}^{N-1}X(m)\cos \frac{(2m+1)\mu \pi }{2N},u=1,2,...,N-1$$

其中

$$C(u)=\{\begin{array}{ll}1/\sqrt{2},& u=0\\ 1,& else\end{array}$$

其逆变换为

$$X(m)=\sqrt{\frac{2}{N}}\sum _{m=0}^{N-1}C(u)Y(u)\cos \frac{(2m+1)u\pi }{2N},m=1,2,...,N-1$$

二者的变换核为

$$a(u,m)=C(u)\sqrt{\frac{2}{N}}\cos \frac{(2m+1)u\pi }{2N}$$

---

二维DCT变换

对二维图像信号数据而言，有推广

$$Y(u,v)=C(u)C(v)\frac{2}{\sqrt{MN}}\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}X(m,n)\cos \frac{(2m+1)u\pi }{2M}\cos \frac{(2n+1)v\pi }{2N},u=1,2,...,M-1,v=0,1,...,N-1$$

$$C(u),C(v)=\{\begin{array}{ll}1/\sqrt{2},& u,v=0\\ 1,& else\end{array}$$

逆变换为

$$X(m,n)=\frac{2}{\sqrt{MN}}\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}C(u)C(v)Y(u,v)\cos \frac{(2m+1)u\pi }{2M}\cos \frac{(2n+1)v\pi }{2N},m=1,2,...,M-1,n=0,1,...,N-1$$

变换核是两次一维DCT变换核的乘积

$$a(u,v,m,n)=a_1(u,m)a_2(v,n)=C(u)\sqrt{\frac{2}{M}}\cos \frac{(2m+1)u\pi }{2M}C(v)\sqrt{\frac{2}{N}}\cos \frac{(2n+1)v\pi }{2N}$$

---

DCT变换在压缩中的应用

DCT的变换核构成的基向量与图像内容无关，而且变换核是可以分离的，既二维DCT可以用两次一维DCT来完成，使得数学运算难度大大简化，再配以已经发现的其它快速算法，使得DCT编码得到了广泛的应用。

特点

• DCT应用于图像数据压缩，可以减少代表图像亮度(或色度)层次数码信息，达到数据压缩的目的。
• 在编码变换过程发现图像细节的位置，以便删去或略去对视觉不敏感的部分，而更加突出视觉的敏感部分，通过选择主要数据来传输、重视图像。

信号的频谱

由图像内取出一个区块，分成8×8个像素的64格阵列，经过对逐个像素的亮度（或讨论色度）数值取样，并将像素的亮度数值列成矩阵形表格。然后利用离散余弦变换（DCT）可将各空间取样值转变为频率域的数值，即为DCT系数。

64点阵列，可得到64个DCT系数，其可以等效为一个直流平均值和63个不同频率余弦波幅值组成的64个点阵列。

这64个点代表频谱上的64个频率，相当于是64个正交基函数的系数。左上角是直流和低频部分，右下角是高频部分。（从表达式也可看出u,v越大的地方基函数频率越高）

可以按照Z字扫描顺序，将DCT系数从二维数据转变为一维数据。

该一维序列的首项是平均亮度值，后面的项代表亮度的起伏程度。

人眼睛对图像的亮度信息有较高的相对灵敏度，对图像的彩色信息不够敏感；还有，人眼睛对图像信息的低频分量具有较高的视觉灵敏度。经Z字形字扫描后所形成的数据系列，恰好与人眼睛对图像信息的敏感程度形成良好的对应关系。根据视觉生理的上述规律，可对图像数据进行压缩。

系数的量化

对每种频率分量设定不同的折算值，将转换得到的DCT系数再次进行折算，可以进一步突出视觉效果影响大的成分，而消弱或忽略视觉效果影响小的成分。这种处理方法称为量化处理，简称Q处理。

如： 将区块的64个系数分别除以量化表中对应位置量化步长，再进行四舍五入取整

对于64点阵列的64个系数来说，对应了64种不同频率，可使用64个不同的折算值。通常称这64个折算值为量化表，每个折算值称为量化步长，或称量化值。在64点阵列中，左上角的数据量化值较小，右下角的数据量化值较大。

经过量化处理后的DCT系数矩阵，可出现许多零值。一般左上角位置的数据的商数是非0，在右下角位置的数据的商数很小，经四舍五入取整值后可简写为0。

果冻君的提问：所以可以压缩成稀疏矩阵吗？？

实际操作归纳

• 要将输入的图像分解为 8*8 的块
• 后对每个块进行二维离散变换把每个块转变成 64 个 DCT 系数值
• 后将变换得到的 DCT 系数进行编码和传送
• 解码时对每个块进行二维 DCT 反变换
• 将反变换后的块组合成图像

将8*8的大小代入原公式，有

$$Y(u,v)=\frac{1}{4}C(u)C(v)\sum _{m=0}^7\sum _{n=0}^{7}X(m,n)\cos \frac{(2m+1)u\pi }{16}\cos \frac{(2n+1)v\pi }{16},u=1,2,...,M-1,v=0,1,...,N-1$$

可以等效为两个N=8的一维DCT变换。

对课程案例的分析

1维DCT变换子函数

其源代码为

void dct_1d(dct_data_t src[DCT_SIZE], dct_data_t dst[DCT_SIZE])
{
   unsigned int k, n;
   int tmp;
   const dct_data_t dct_coeff_table[DCT_SIZE][DCT_SIZE] = {
#include "dct_coeff_table.txt"
   };
DCT_Outer_Loop:
   for (k = 0; k < DCT_SIZE; k++) {
DCT_Inner_Loop:
      for(n = 0, tmp = 0; n < DCT_SIZE; n++) {
         int coeff = (int)dct_coeff_table[k][n];
         tmp += src[n] * coeff;
      }
      dst[k] = DESCALE(tmp, CONST_BITS);
   }
}

引用的系数表为

  8192,  8192,  8192,  8192,  8192,  8192,  8192,  8192,
 11363,  9633,  6436,  2260, -2260, -6436, -9632,-11362,
 10703,  4433, -4433,-10703,-10703, -4433,  4433, 10703,
  9633, -2260,-11362, -6436,  6436, 11363,  2260, -9632,
  8192, -8192, -8192,  8192,  8192, -8191, -8191,  8192,
  6436,-11362,  2260,  9633, -9632, -2260, 11363, -6436,
  4433,-10703, 10703, -4433, -4433, 10703,-10703,  4433,
  2260, -6436,  9633,-11362, 11363, -9632,  6436, -2260

由代码可知，输入数据和输出数据都是一维数组，其首先将系数表文件导入数组中，然后对数据进行遍历。
Outer循环对k，也就是系数表的行进行遍历
Inner循环对数据和系数表的列n进行遍历
因此其每次Inner循环都是将数据src中的每个值与系数表第k行的每个数值相乘再相加，就得到了输出数据的第k个值。

从一维DCT变换的公式可知，内层循环做累加时所乘系数应该为

$$a(u,m)=C(u)\sqrt{\frac{2}{N}}\cos \frac{(2m+1)\mu \pi }{2N}$$

在次代码中，u为j，m为n。

最后的Descale的宏定义为

#define DESCALE(x,n)  (((x) + (1 << ((n)-1))) >> n)

其效果相当于给x加一个2^(n-1)的常量便宜再除以2ⁿ
很多本来就不够大的数字会被移位变成0。

二维DCT子函数

其源码为

void dct_2d(dct_data_t in_block[DCT_SIZE][DCT_SIZE],
      dct_data_t out_block[DCT_SIZE][DCT_SIZE])
{
   dct_data_t row_outbuf[DCT_SIZE][DCT_SIZE];
   dct_data_t col_outbuf[DCT_SIZE][DCT_SIZE], col_inbuf[DCT_SIZE][DCT_SIZE];
   unsigned i, j;

   // DCT rows
Row_DCT_Loop:
   for(i = 0; i < DCT_SIZE; i++) {
      dct_1d(in_block[i], row_outbuf[i]);
   }
   // Transpose data in order to re-use 1D DCT code
Xpose_Row_Outer_Loop:
   for (j = 0; j < DCT_SIZE; j++)
Xpose_Row_Inner_Loop:
      for(i = 0; i < DCT_SIZE; i++)
         col_inbuf[j][i] = row_outbuf[i][j];
   // DCT columns
Col_DCT_Loop:
   for (i = 0; i < DCT_SIZE; i++) {
      dct_1d(col_inbuf[i], col_outbuf[i]);
   }
   // Transpose data back into natural order
Xpose_Col_Outer_Loop:
   for (j = 0; j < DCT_SIZE; j++)
Xpose_Col_Inner_Loop:
      for(i = 0; i < DCT_SIZE; i++)
         out_block[j][i] = col_outbuf[i][j];
}

输入数据为二维数组，输出数据也是二维数组。
其首先进行Row_DCT_Loop，对ROWs进行DCT变换
接下来的Xpose循环实际上是对数组进行了转置，这样新的数组相当于是进行了Cols的DCT变换
然后再进行一次DCT循环。
最后再转置变回原来的方向。

参考文献

[1]孙娟. DCT原理及其在视频压缩编码中的实现[EB/OL]. 北京：中国科技论文在线 [2008-07-23].
[2] https://www.cnblogs.com/stnlcd/p/7261842.html