H264-整数DCT变换和蝶形变换代码实现

背景

为压缩数据，H264标准需要使用整数DCT变换把图像信号转换到频域。变换过程实际上就是矩阵乘法，使用蝶形变换可以减少运算量，提升算法性能。

原理

整数DCT变换公式推导过程不再赘述，可以参考《深入理解视频编解码技术 --基于H264标准及参考模型》。
以4x4为例，变换公式如下：

  Y = Cf * X * Cf.T
  其中 Cf = 
  	[1, 1,  1,  1,
     2, 1, -1, -2,
     1,-1, -1,  1,
     1,-2,  2, -1]

展开如下：

• 先以Cf*X 为例，矩阵乘需要进行64次乘法和48次加法，这将耗费较多的时间。蝶形算法则是利用构造的矩阵的整数性质和对称性，将乘法运算转化为加法运算。
  先分析两个相乘矩阵的第一列，结果如下：
  p0 = $x00+x30$ + $x10+x20$
  p1 = $2\ast x00+2\ast x30$ + $x10-x20$
  p2 = $x00-x30$ - $x10-x20$
  p3 = $x00-x30$ - $2\ast x10+2\ast x20$
• 从上面公式我们发现有重复运算（使用不同颜色标识），从一般算法意义上来说，可以用空间代价换时间代价，比如设置中间变量来减少计算次数。
  提取中间结果：
  M0 = x00 + x30
  M3 = x00 - x30
  M1 = x10 + x20
  M2 = x10 - x20
• 存储了这四个中间变量，我们对比看看蝶形图，和图中第一层的算式相符合。用这些中间变量来组合，就可以把最终的p0…p3, 计算出来。这样，就把运算量降低到2个乘法和8个加法，剩余的运算就是叠代这个算法。
  

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

• 这个蝶形图和一般意义的FFT或FDCT蝶形图不同，是对H.264在整数DCT基础上的具体算法优化，只对于以上Cf矩阵。
• 计算过程是把上面的三个4x4矩阵乘法分成两两矩阵相乘。再把残差矩阵和后来的中间结果Cf x X一行行分别输入蝶形图进行一维整数DCT计算。
• 蝶形图优化思想就是提取矩阵的相关部分，定义中间变量，减少运算次数。

代码示例

以下代码分别用传统矩阵运算和蝶形运算计算结果并对比结果

#include <iostream>
using namespace std;

//Cf矩阵
int Cf[4][4]=
{
    1, 1,  1,  1,
    2, 1, -1, -2,
    1,-1, -1,  1,
    1,-2,  2, -1
};

void matrixTransform(int a[][4])
{
	int i, j;
	int tmp;
	for(i=0; i<4; i++)
	{
		for(j=i; j<4; j++)
		{
			tmp = a[i][j];
			a[i][j] = a[j][i];
			a[j][i] = tmp;
		}
	}
}

//矩阵求积
// c = a * b
void  matrixMultiply_4x4(int a[][4],int b[][4],int c[][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];  
			}
		}
	}      
}

//矩阵显示
void matrixShow(int a[][4])
{
    int i, j;
    cout << "*****************************" << endl;
    for(i = 0; i < 4; i++)
	{
        for(j = 0; j < 4; j++)
        {
            cout << a[i][j] << "\t";
        }
		cout << endl;
	}
 
    cout << "*****************************" << endl << endl;
}

// 蝶形变换求解 Cf* orgYUV * Cf.T
// b为输出matrix， a为输入 
void forward4x4(int a[][4], int b[][4])
{
	int tmp[16];
	int p0,p1,p2,p3;
	int t0,t1,t2,t3;
	int i, j;
	int *ptmp = tmp;
	
	//horizontal， Cf*orgYUV 
	for(i=0; i<4; i++)
	{
		t0 = a[0][i];
		t1 = a[1][i];
		t2 = a[2][i];
		t3 = a[3][i];
		
		p0 = t0+t3;
		p1 = t1+t2;
		p2 = t1-t2;
		p3 = t0-t3;
		
		*(ptmp++) = p0+p1;
		*(ptmp++) = p2+ (p3<<1); // 加运算符优先级要于左移 
		*(ptmp++) = p0-p1;
		*(ptmp++) = p3- (p2<<1);
	}	
	
	//vertical
	for(int i=0; i<4; i++)
	{
		t0 = tmp[i];
		t1 = tmp[i+4];
		t2 = tmp[i+2*4];
		t3 = tmp[i+3*4];
		
		p0 = t0+t3;
		p1 = t1+t2;
		p2 = t1-t2;
		p3 = t0-t3;
		
		b[i][0] = p0+p1;
		b[i][1] = p2+ (p3<<1);
		b[i][2] = p0-p1;
		b[i][3] = p3- (p2<<1);
	}
}

int main(int argc, char** argv) {
	//原始的YUV矩阵
	int orgYUV[4][4] =
	{
		-85, 88, 126, 121,
	    -79, 70,  65,  83,
        -80, 66,  49,  43,
        -82, 86,  97,  41
	};
	
	int D[4][4];
	int W[4][4];
	// W = Cf* orgYUV * Cf.T
	matrixMultiply_4x4(Cf, orgYUV, D);
	matrixTransform(Cf);	
	matrixMultiply_4x4(D, Cf, W);

	cout << "====== W =====" << endl;
	matrixShow(W);
	
	forward4x4(orgYUV, W);
	cout << "====== butter W =====" << endl;
	matrixShow(W);
	
	return 0;
}