【转】H.264，HEVC下整数4*4DCT，量化，反量化，IDCT过程简单讲解

通过此贴可以加深了解神术DCT变换以及量化的过程，下面以4X4为变换单元来做。
//DCT变换、量化、反DCT变换、反量化。    


#include <stdio.h>
#include <math.h>


int QP=27;        //quantize parameter

int D[4][4];　　　　
int W[4][4];
int Z[4][4];
int Wi[4][4];
int Di[4][4];
int Xi[4][4];
//左乘数组
int Cf[4][4]=
        {
                {1, 1,  1, 1},
                {2, 1, -1, -2},
                {1,-1, -1, 1},
                {1,-2,  2, -1}
        };
        //右乘转置数组int Cft[4][4]=
        {
                {1, 2,  1, 1 },
                {1, 1, -1, -2},
                {1,-1, -1, 2},
                {1,-2,  1, -1}
        };
　　　　　//逆变换数组int Ci[4][4]=
        {
                {2,  2,  2, 1 },
                {2,  1, -2, -2},
                {2, -1, -2, 2 },
                {2, -2,  2, -1}
        };
        
        int Cit[4][4]=
        {
                {2 ,  2 ,  2,  2 },
                {2 ,  1 , -1, -2 },
                {2 , -2 , -2,  2 },
                {1 , -2 ,  2, -1 }
        };

        int MF[6][6]=
        {
                {13107, 5243, 8066},
                {11916, 4660, 7490},
                {10082, 4194, 6554},
                {9362,  3647, 5825},
                {8192,  3355, 5243},
                {7282,  2893, 4559}
        };

        int V[6][6]=
        {
                {10, 16, 13 },
                {11, 18, 14 },
                {13, 20, 16 },
                {14, 23, 18 },
                {16, 25, 20 },
                {18, 29, 23 }
        };


        //X就是输入的残差矩阵
        
         int X[4][4]=
         {
                 {-4 , 4, -14 , -25},
                 {-4 , 4 , -13 , -23},
                 {-4 , 3, -14, -21 },
                 {-4, 3,  -7, -14 }
 
         };

void main()
{
        int i,j;

        //函数声明
        int matrix_multiply(int a[4][4],int b[4][4],int c[4][4]);
        int matrix_substract(int a[4][4],int b[4][4]);
        void matrix_display(int a[4][4]);
        void prescale_quantize(int W[4][4]);
        void rescale(int Z[4][4]);

        //X就是输入的残差矩阵
        printf("\nDisplay the Resident matrix X[4][4]:\n");
        // 单纯的显示所输入的数据，没有任何意义
        matrix_display(X);


        matrix_multiply(Cf,X,D);        


        
        matrix_multiply(D,Cft,W);        

        printf("\nDisplay the Transformed matrix W[4][4]:\n");
        matrix_display(W);

        //量化
        prescale_quantize(W);  //quantize

        printf("\nDisplay the quantized matrix Z[4][4]:\n");
        matrix_display(Z);
        //z是DCT变换+量化之后的结果 
        

        // decode
        //逆量化
        //inverse quantize
　　　　　rescale(Z); 
        printf("\nDisplay the I-quantized matrix Wi[4][4]:\n");
        matrix_display(Wi);

        //inverse "core" DCT transform
        matrix_multiply(Ci,Wi,Di);                //        
        matrix_multiply(Di,Cit,Xi);                //
 
        for(i=0;i<4;i++)
                for(j=0;j<4;j++)
                        Xi[i][j]=(int)(((float)Xi[i][j]/256)+0.5);        //divide 64*4=256

        printf("\nDisplay the I-transformed matrix Xi[4][4]:\n");
        //XI = IDCT后数据
        matrix_display(Xi);

        // IDCT后的数据 减去原数据之差
        matrix_substract(Xi,X);
        printf("\nDisplay the difference matrix Xi[4][4]-X[4][4]:\n");
        matrix_display(Xi);
}

// multiply the matrix of a and 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];                                
                }
        return 0;
        
}

// multiply the matrix of a and b
int  matrix_substract(int a[4][4],int b[4][4])
{
        int i,j;
        
        for(i=0;i<4;i++)
                for(j=0;j<4;j++)
                {
                        a[i][j]-=b[i][j];
                }
        return 0;
}

// 单纯的显示所输入的数据，没有任何意义。
void matrix_display(int a[4][4])
{
        int i,j;
        printf("=========================================\n");
    for(i=0;i<4;i++)
                for(j=0;j<4;j++)
                {
                        printf("%d\t",(int)a[i][j]);
                        if(j==3) 
                                printf("\n");
                }
        printf("=========================================\n");
}

//量化
void prescale_quantize(int W[4][4])
{

        int i,j;
        for(i=0;i<4;i++)
                for(j=0;j<4;j++)
                {
                        int sign=1;
                        int qbits=15+floor(QP/6.0);
                        //int qbits=16;  //2005.6.21:: why the qbits is 16(H.264 white-paper) not 15?
                        int mf,k;        //
                        
                        int f=(int)(pow(2.0,qbits)/3);
                        
                        if(W[i][j]<0)
                        {
                                sign=-1;
                                W[i][j]*=-1;
                        }
                        
                        if((i==0||i==2)&&(j==0||j==2))
                                k=0;
                        else
                                if((i==1||i==3)&&(j==1||j==3))
                                        k=1;
                                else
                                        k=2;
                                
                        mf=MF[QP%6][k];
                                
                        Z[i][j]=(W[i][j]*mf+f)>>qbits;
                        Z[i][j]*=sign;
                }        
}

//逆量化
void rescale(int Z[4][4])
{
        int i,j;
        for(i=0;i<4;i++)
                for(j=0;j<4;j++)
                {
                        int v,k,f,t;        

                        t = floor(QP/6.0);                
                        f=(int)pow(2.0,t);                        

                        if((i==0||i==2)&&(j==0||j==2))
                                k=0;
                        else
                                if((i==1||i==3)&&(j==1||j==3))
                                        k=1;
                                else
                                        k=2;

                        v=V[QP%6][k];
                        Wi[i][j]=Z[i][j]*v*f;
                }
}