程序

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


static const double PI=3.1415;
#define N 4


float Y[4][4];
float H[4][4],Ht[4][4];
float Q1[4][4],Q2[4][4];
float R[N][N];
float D[4][4],E[4][4];
float X[4][4]=
{  
        {17, 2, 31, 49},  
        {9, 8, 7, 6},  
        {15, 30, 22, 7},  
        {18, 50, 100, 5}  
    };  


float A[4][4]=
{  
        {1, 1,  1, 1},  
        {2, 1, -1, -2},  
        {1,-1, -1, 1},  
        {1,-2,  2, -1}
    };  
    
float At[4][4]=
{  
    {1, 2,  1, 1 },  
    {1, 1, -1, -2},  
    {1,-1, -1, 2},  
    {1,-2,  1, -1}  
};   
/*整数dct变换中需要点乘的系数*/ 
float G1[4][4]={
    {0.25,0.25*0.6325,0.25,0.25*0.6325},
    {0.25*0.6325,0.1,0.25*0.6325,0.1},
    {0.25,0.25*0.6325,0.25,0.25*0.6325},
    {0.25*0.6325,0.1,0.25*0.6325,0.1},
};


void init_dct(float A[N][N]){
    int i,j;
    for(j=0;j<N;j++){
        A[0][j]=sqrt(1.0/N);
    }
    for(i=1;i<N;i++){
        for(j=0;j<N;j++){
            A[i][j] = sqrt(2.0/N)*cos((2*j+1)*i*PI/(2*N));
        }
    }
}


void tran(float a[N][N], float b[N][N])
{
    int i,j;
    for(i=0;i<N;i++)
    for(j=0;j<N;j++)
    a[i][j]=b[j][i];
}


void dct(float X[4][4]){
    int i,j,k;
    float B[4][4];
    float M[4][4];
    for(j=0;j<4;j++){
        B[0][j]=X[0][j]+X[3][j];
        B[1][j]=X[1][j]+X[2][j];
        B[2][j]=X[1][j]-X[2][j];
        B[3][j]=X[0][j]-X[3][j];
    }
    for(j=0;j<4;j++){
        M[0][j]=B[0][j]+B[1][j];
        M[1][j]=2*B[3][j]+B[2][j];
        M[2][j]=B[0][j]-B[1][j];
        M[3][j]=B[3][j]-2*B[2][j];
    }
    for(i=0;i<4;i++){
        B[i][0]=M[i][0]+M[i][3];
        B[i][1]=M[i][1]+M[i][2];
        B[i][2]=M[i][1]-M[i][2];
        B[i][3]=M[i][0]-M[i][3];
    }
    for(i=0;i<4;i++){
        Y[i][0]=B[i][0]+B[i][1];
        Y[i][1]=2*B[i][3]+B[i][2];
        Y[i][2]=B[i][0]-B[i][1];
        Y[i][3]=B[i][3]-2*B[i][2];
    }
    for(i=0;i<4;i++){
        for(j=0;j<4;j++){
            Y[i][j]*=G1[i][j];
        }
    }
}
 
void display_martix(float T[N][N]){
    int i,j;
    for(i=0;i<N;i++){
        for(j=0;j<N;j++){
            printf("%8.3f\t",T[i][j]);
            if(j!=0 && j%3==0)printf("\n");
        }
    }
}


void dct1(float X[4][4]){
    int i,j,k;
    init_dct(H);
    tran(Ht,H);
    for(i=0;i<4;i++){
        for(j=0;j<4;j++){
            for(k=0;k<4;k++){
                Q1[i][j]+=H[i][k]*X[k][j];
            }
        }
    }  
    for(i=0;i<4;i++){
        for(j=0;j<4;j++){
            for(k=0;k<4;k++){
                Q2[i][j]+=Q1[i][k]*Ht[k][j];
            }
        }
    }  
}


void idct(float Y[4][4]){
    int i,j,k;
    for(i=0;i<4;i++){
        for(j=0;j<4;j++){
            for(k=0;k<4;k++){
                D[i][j]+=Ht[i][k]*Y[k][j];
            }
        }
    }  
    for(i=0;i<4;i++){
        for(j=0;j<4;j++){
            for(k=0;k<4;k++){
                E[i][j]+=D[i][k]*H[k][j];
            }
        }
    }  
}


void minus_matrix(float A[N][N],float B[N][N],float C[N][N]){
    int i,j;
    for(i=0;i<N;i++){
        for(j=0;j<N;j++){
            C[i][j]=A[i][j]-B[i][j];
        }
    }
} 


 void main(){


  /*采用普通的dct变换，矩阵相乘的方式,Q2为变换后的结果*/
  dct1(X);
  printf("dct系数H[i][j]\n");
  display_martix(H);
  printf("dct转置系数Ht[i][j]\n");
  display_martix(Ht);
  printf("普通dct结果Q2[i][j]\n");
  display_martix(Q2);
  
  /*采用整数dct变换，矩阵相乘的方式,Y为变换后的结果*/
  dct(X);
  printf("整数dct结果Y[i][j]\n");
  display_martix(Y);
 
  /*采用dct逆变换,E为变换后的结果*/
  idct(Y);  
  printf("逆变换结果E[i][j]\n");
  display_martix(E);
 
  /*计算dct和整数dct的差值矩阵*/
  minus_matrix(X,E,R);
  printf("整数dct和dct的差值矩阵R[i][j]\n");
  display_martix(R); 

}