H264___DCT蝶形算法____理解

转自：http://fellowher.blog.163.com/blog/static/5076351520099279115651/  中间有自己的理解。

 为什么要进行变换

空间图像数据通常是很难压缩的：相邻的采样点具有很强的相关性（相互关联的），而且能量一般平均分布在一幅图像中，从而要想丢掉某些数据和降低数据精度而不明显影响图像质量，就要选择合适的变换，方法，使图像易于被压缩。适合压缩的变换方法要有这样几个性质：

1、可以聚集图像的能量（将能量集中到少数有意义的数值上）；如下图：

举个例子说明，下图是变换前的数据：

以下是变换后的数据：

可以看出，经变换后，数据的能量基本上集中到左上方（低频信号）了，而变换后的数据完全可以通过反变换还原成原来的数据。为了达到压缩文件的目的，我们就可以丢弃掉一些能量低的数据（高频信号），而对图像质量影响很小。

2、可以除去数据之间的相关性（以使丢掉不重要的数据对图像的质量影响很少）。

3、变换方法应该适合用软硬件实现。

 

以下将描述4x4整数DCT推导过程

 

 

将上面的DCT公式转变为矩阵形式，为了说明标准中整数变换和反变换，先设d=c/b,

[sno: 这里两个［］中间是矩阵乘的关系]

我们可以把DCT变换写成下式：

［sno:这是三个矩阵相乘，就是DCT转换的公式：Cf * X *Cft  ］

其中，为了保证正交，a ,b, c, d的值可分别如下设置:

根据假设，上式可化为：

把两边的对称矩阵移到左边可得：（对角阵移项）

这一步也有点问题：（等价的乘数移项）

 

此时，d的值为0.4142。这样的话，还是实数运算。如我们令d=1/2,则b=sqrt(2/5).

c=sqrt(1/10),同样，可以保证矩阵的正交，同时，可以把运算变为整数运算。

1/2可以提到矩阵外面，并与右边的点乘合并，得

其中，

在JM编码器中，变换过程只包括了

后面的点乘实际上是在量化过程中进行，因为后面的点乘还有实数运算，实数运算将不可避免地产生精度误差,而且运算量巨大。而量化本身就会丢失一些信号，因些，这些实数运算放在量化过程中将大大的降低变换的运算率同时又不明显影响精度？[sno: 之所有把后面的Ef和这里分割开，是因为在转换和量化这一步，不可避免的会遇到乘法操作，这里说的只用加减和移位操作，是在这三个矩阵乘法大数据量运算之间只使用加减和移位，可以把64次的乘法运算转为只有两次的乘法运算，极大的减少消耗，在量化的时候，再对每个W[i][j]分别进行乘法操作。]

 

然而，4X4的矩阵运算如果按常规算法的话仍要进行64次乘法运算和48次加法，这将耗费较多的时间，于是在H.264中，有一种改进的算法（蝶形算法）可以减少运算的次数。这种矩阵运算算法构造非常巧妙，利用构造的矩阵的整数性质和对称性，可完全将乘法运算转化为加法运算。

 

变换过程在JM中代码实现如下：

  //  Horizontal transform水平变换，其实就是左乘Cf,

  for (j=0; j < BLOCK_SIZE && !lossless_qpprime; j++)

  {

            for (i=0; i < 2; i++)

    {

      i1=3-i;

      m5[i]=img->m7[i][j]+img->m7[i1][j];

      m5[i1]=img->m7[i][j]-img->m7[i1][j];

    }

    img->m7[0][j]=(m5[0]+m5[1]);

    img->m7[2][j]=(m5[0]-m5[1]);

    img->m7[1][j]=m5[3]*2+m5[2];

    img->m7[3][j]=m5[3]-m5[2]*2;

  }

  //  Vertical transform垂直变换，其实就是右乘CfT

  for (i=0; i < BLOCK_SIZE && !lossless_qpprime; i++)

  {

            for (j=0; j < 2; j++)

    {

      j1=3-j;

      m5[j]=img->m7[i][j]+img->m7[i][j1];

      m5[j1]=img->m7[i][j]-img->m7[i][j1];

    }

    img->m7[i][0]=(m5[0]+m5[1]);

    img->m7[i][2]=(m5[0]-m5[1]);

    img->m7[i][1]=m5[3]*2+m5[2];

    img->m7[i][3]=m5[3]-m5[2]*2;

  }

上面的内容载抄自网络，我做了少量更改。多数内容来自毕厚杰的书。

这里来分析一下蝶形算法，这个蝶形算法和一般FFT的蝶形算法不同，由于我没有找到相关论文，能找到的书和网络资料又语焉不详，只好自己推导。上面的JM代码就是计算下面三个4x4矩阵的过程。

 

分析一下前两个矩阵的乘法，只分析他们结果矩阵的第一行。有什么办法可以减少运算量呢？首先采用传统方法计算，得到结果：

X[0] = x[00]+x[10]+x[20]+x[30]

X[1] = 2*x[00]+x[10]-x[20]-2*x[30]

X[2]= x[00]-x[10]-x[20]+x[30]

X[3] = x[00]-2x[10]+2x[20]-x[30]

 

计算代价是16次乘法12次加法，考虑到矩阵的1的乘法可以省略，去除8个乘1，还需要8次乘法和12次加法。那么我们再仔细思考他们的相关性，从一般算法意义上来说，可以用空间代价换时间代价，比如设置中间变量来减少计算次数。用不同的颜色把需要重复运算的部分标上，作为中间变量。

X[0] = x[00]+x[10]+x[20]+x[30]

X[1] = 2*x[00]+x[10]-x[20]-2*x[30]

X[2]= x[00]-x[10]-x[20]+x[30]

X[3] = x[00]-2x[10]+2x[20]-x[30]

那么提取出来的中间变量将是：

x[00]+x[30]

x[00]-x[30]

x[10]+x[20]

x[10]-x[20]

存储了这四个中间变量，我们对比看看蝶形图，和图中第一层的算式相符合。用这些中间变量来组合，就可以把最终的X[0]..X[3], 计算出来。这样，就把运算量降低到2个乘法和8个加法，剩余的运算就是叠代这个算法。

所以，可以得出以下结论：

1. 这个蝶形图和一般意义的FFT或FDCT蝶形图不同，是对H.264在整数DCT基础上的具体算法优化，只对于以上Cf矩阵。
2. 计算过程是把上面的三个4x4矩阵乘法分成两两矩阵相乘。再把残差矩阵和后来的中间结果Cf x X一行行分别输入蝶形图进行一维整数DCT计算。
3. 蝶形图优化思想就是提取矩阵的相关部分，定义中间变量，减少运算次数。［sno：这里很重要，之所以使用蝶形算法，就是统计一下多次运行的规律，不重复计算相关的量。］
4. sno: 矩阵相乘 A*B   就是把A矩阵当前元素所在的这一行乘以B矩阵当前元素所在的这一列。代码实现如下
   int  matrix_multiply(int a[4][4],int b[4][4],int c[4][4])
   {
           int i,j,k;
           for(i=0;i<4;i++)
                   for(j=0;j<4;j++)
                   {
                           c[i][j]=0;
                           for(k=0;k<4;k++)        
                                   c[i][j]+=a[i][k]*b[k][j];  ///c[i][j]=a的i这一行*b的j这一列的总和 ,这里是累计和                          
                   }
           return 0;        
   }