DCT算法可以说是有损压缩的第一步，多用在视频压缩方面。

转自http://m.blog.csdn.net/blog/zhu530548851/21444459

 DCT算法可以说是有损压缩的第一步，多用在视频压缩方面。它将二维转化为一维数据，并且将能量聚集到左上角。关于DCT，有很多快速算法，我个人认为，我这个算法还是比较快速的。下面是我的代码：

#include<stdio.h>
#include <math.h>
#include<string.h>
#include<stdlib.h>
#define PI 3.1415926
double *C=NULL;
double *temp_2D=NULL;
#define N 256




void fdct_2D(double *f,int deg_row,int deg_col);
void fidct_2D(double *F,int deg_row,int deg_col);
void init2D_Param(int rows,int cols);
void fdct_1D_no_param(double *f,int deg);
void initDCTParam(int deg);
void dct_forward(double *f,int deg);
void dct_backward(double *f,int deg);
void fidct_1D_no_param(double *F,int deg);
void initIDCTParam(int deg);
void idct_backward(double *F,int deg);
void idct_forward(double *F,int deg);
void fbitrev(double *f,int deg);
int bitrev(int bi,int deg);
void swap(double *a,double *b);
int GetTwoIndex(int num);




void fdct_2D(double *f,int deg_row,int deg_col)
{    

int rows,cols,i_row,i_col;
double two_div_sqrtcolrow;
rows=1<<deg_row;
cols=1<<deg_col;
init2D_Param(rows,cols);
two_div_sqrtcolrow=2.0/(sqrt((double)(rows*cols)));

for(i_row=0;i_row<rows;i_row++)
{
fdct_1D_no_param(f+i_row*cols,deg_col);
}
for(i_col=0;i_col<cols;i_col++)
{
for(i_row=0;i_row<rows;i_row++)
{
temp_2D[i_row]=f[i_row*cols+i_col];
}
fdct_1D_no_param(temp_2D,deg_row);
for(i_row=0;i_row<rows;i_row++)
{
f[i_row*cols+i_col]=temp_2D[i_row]*two_div_sqrtcolrow;
}

}

}
void fidct_2D(double *F,int deg_row,int deg_col)
{    

int rows,cols,i_row,i_col;
double sqrtcolrow_div_two;
rows=1<<deg_row;
cols=1<<deg_col;
init2D_Param(rows,cols);
sqrtcolrow_div_two=(sqrt((double)(rows*cols)))/2.0;

for(i_row=0;i_row<rows;i_row++)
{
fidct_1D_no_param(F+i_row*cols,deg_col);
}
for(i_col=0;i_col<cols;i_col++)
{
for(i_row=0;i_row<rows;i_row++)
{
temp_2D[i_row]=F[i_row*cols+i_col];
}
fidct_1D_no_param(temp_2D,deg_row);
for(i_row=0;i_row<rows;i_row++)
{
F[i_row*cols+i_col]=temp_2D[i_row]*sqrtcolrow_div_two;
}

}

}
void init2D_Param(int rows,int cols)
{
if(temp_2D!=NULL) free(temp_2D);
temp_2D=(double *)malloc(sizeof(double)*rows);
}
void fdct_1D_no_param(double *f,int deg)
{
initDCTParam(deg);
dct_forward(f,deg);
dct_backward(f,deg);
fbitrev(f,deg);
f[0]=1/(sqrt(2.0))*f[0];
}
void initDCTParam(int deg)
{
int total,halftotal,i,group,endstart,factor,j,len;
total=1<<deg;
if(C!=NULL) free(C);
C=(double *)malloc(sizeof(double)*total);
halftotal=total>>1;
for(i=0;i<halftotal;i++)
C[total-i-1]=(double)(2*i+1);

for(group=0;group<deg-1;group++)
{ 
endstart=1<<(deg-1-group);
len=endstart>>1;
factor=1<<(group+1);
for(j=0;j<len;j++)
C[endstart-j-1]=factor*C[total-j-1];
}
for(i=1;i<total;i++)
C[i]=2.0*cos(C[i]*PI/(total<<1));

}
void dct_forward(double *f,int deg)
{
int i_deg,i_halfwing,total,wing,wings,winglen,halfwing;
double temp1,temp2;
total=1<<deg;
for(i_deg=0;i_deg<deg;i_deg++)
{
wings=1<<i_deg;
winglen=total>>i_deg;
halfwing=winglen>>1;
for(wing=0;wing<wings;wing++)
{
for(i_halfwing=0;i_halfwing<halfwing;i_halfwing++)
{
temp1=f[wing*winglen+i_halfwing];
temp2=f[(wing+1)*winglen-1-i_halfwing];
if(wing%2)
swap(&temp1,&temp2);
f[wing*winglen+i_halfwing]=temp1+temp2;
f[(wing+1)*winglen-1-i_halfwing]=(temp1-temp2)*C[winglen-1-i_halfwing];
}
}
}
}
void swap(double *a,double *b)
{
double temp;
temp=*a;
*a=*b;
*b=temp;
}
void dct_backward(double *f,int deg)
{
int total,i_deg,wing,wings,halfwing,winglen,i_halfwing,temp1,temp2;
    total=1<<deg;
for(i_deg=deg-1;i_deg>=0;i_deg--)
{
wings=1<<i_deg;
winglen=1<<(deg-i_deg);
halfwing=winglen>>1;
for(wing=0;wing<wings;wing++)
{
for(i_halfwing=0;i_halfwing<halfwing;i_halfwing++)
{  //f[i_halfwing+wing*winglen]=f[i_halfwing+wing*winglen];
if(i_halfwing==0)
f[halfwing+wing*winglen+i_halfwing]=0.5*f[halfwing+wing*winglen+i_halfwing];
else
{
temp1=bitrev(i_halfwing,deg-i_deg-1);
temp2=bitrev(i_halfwing-1,deg-i_deg-1);
f[halfwing+wing*winglen+temp1]=f[halfwing+wing*winglen+temp1]-f[halfwing+wing*winglen+temp2];
}
}
}
}
}
void fidct_1D_no_param(double *F,int deg)
{
initIDCTParam(deg);
F[0]=F[0]*sqrt(2.0);
fbitrev(F,deg);
idct_forward(F,deg);
idct_backward(F,deg);

}
void fbitrev(double *f,int deg)
{   
int i,ii,len;
double temp;
if(deg==1) return;
len=(1<<deg)-1;
i=1;
while(i<len)
{   ii=bitrev(i,deg); 
if(ii>i)
{  temp=f[ii];
 f[ii]=f[i];
 f[i]=temp;
}
i++;
}
}
void initIDCTParam(int deg)
{
int total,halftotal,i,group,endstart,factor,j,len;
total=1<<deg;
if(C!=NULL) free(C);
C=(double *)malloc(sizeof(double)*total);
halftotal=total>>1;
for(i=0;i<halftotal;i++)
C[total-i-1]=(double)(2*i+1);

for(group=0;group<deg-1;group++)
{ 
endstart=1<<(deg-1-group);
len=endstart>>1;
factor=1<<(group+1);
for(j=0;j<len;j++)
C[endstart-j-1]=factor*C[total-j-1];
}
for(i=1;i<total;i++)
C[i]=1.0/(2.0*cos(C[i]*PI/(total<<1)));

}
void idct_backward(double *F,int deg)
{
int i_deg,i_halfwing,total,wing,wings,winglen,halfwing;
double temp1,temp2;
total=1<<deg;
for(i_deg=deg-1;i_deg>=0;i_deg--)
{
wings=1<<i_deg;
winglen=total>>i_deg;
halfwing=winglen>>1;
for(wing=0;wing<wings;wing++)
{
for(i_halfwing=0;i_halfwing<halfwing;i_halfwing++)
{
temp1=F[wing*winglen+i_halfwing];
temp2=F[(wing+1)*winglen-1-i_halfwing]*C[winglen-1-i_halfwing];
                if(wing%2)
{
F[wing*winglen+i_halfwing]=(temp1-temp2)*0.5;
F[(wing+1)*winglen-1-i_halfwing]=(temp1+temp2)*0.5;
}
else
{
F[wing*winglen+i_halfwing]=(temp1+temp2)*0.5;
F[(wing+1)*winglen-1-i_halfwing]=(temp1-temp2)*0.5;
}

}
}
}

}
void idct_forward(double *F,int deg)
{
int total,i_deg,wing,wings,halfwing,winglen,i_halfwing,temp1,temp2;
    total=1<<deg;
for(i_deg=0;i_deg<deg;i_deg++)
{
wings=1<<i_deg;
winglen=1<<(deg-i_deg);
halfwing=winglen>>1;
for(wing=0;wing<wings;wing++)
{
for(i_halfwing=halfwing-1;i_halfwing>=0;i_halfwing--)
{
if(i_halfwing==0)
F[halfwing+wing*winglen+i_halfwing]=2.0*F[halfwing+wing*winglen+i_halfwing];
else
{ 
temp1=bitrev(i_halfwing,deg-i_deg-1);
temp2=bitrev(i_halfwing-1,deg-i_deg-1);
F[halfwing+wing*winglen+temp1]=F[halfwing+wing*winglen+temp1]+F[halfwing+wing*winglen+temp2];
}

}
}

}

} 
int bitrev(int bi,int deg)
{
int j=1,temp=0,degnum,halfnum;
degnum=deg;
//if(deg<0) return 0;
if(deg==0) return bi;
halfnum=1<<(deg-1);
while(halfnum)
{
if(halfnum&bi)
temp+=j;
j<<=1;
halfnum>>=1;
}
return temp; 
}
int GetTwoIndex(int num)
{
int temp=0;
if(num<=0) return -1;
while(1<<temp<num)
temp++;
return temp;


}
int main()
{
int i;
double *data,*data2;
int rol;
data=(double *)malloc(N*sizeof(double));

rol=sqrt((double)N);
data2=(double *)malloc(rol*sizeof(double));
for(i=0;i<N;i++)
{
data[i]=-i;
if(i%rol==0&&i!=0)
{
printf("\n");
}
printf("%4d",(int)data[i]);
}
printf("\nFDCT begin!\n");
fdct_2D(data,GetTwoIndex(rol),GetTwoIndex(rol));
for(i=0;i<N;i++)
{
if(i%rol==0&&i!=0)
{
printf("\n");
}
printf("%4d",(int)data[i]);
}


printf("\n2nd FDCT begins!\n");
for(i=0;i<N;i++)
{
if(i%(int)sqrt(N)==0)
{
data2[(int)(i/sqrt(N))]=data[i];
printf("%4d",(int)data2[(int)(i/sqrt(N))]);
if(((int)(i/sqrt(N)))%4==3)
{
printf("\n");
}
}


}
printf("\n");
fdct_2D(data2,GetTwoIndex(sqrt(rol)),GetTwoIndex(sqrt(rol)));
for(i=0;i<rol;i++)
{
printf("%4d",(int)data2[i]);
if(i%4==3)
{
printf("\n");
}
}
printf("\n");
fdct_2D(data2,GetTwoIndex(sqrt(rol)),GetTwoIndex(sqrt(rol)));
for(i=0;i<rol;i++)
{
printf("%4d",(int)data2[i]);
if(i%4==3)
{
printf("\n");
}
}




printf("\nFIDCT begins!\n");
fidct_2D(data,GetTwoIndex(rol),GetTwoIndex(rol));
for(i=0;i<N;i++)
{
if(i%rol==0&&i!=0)
{
printf("\n");
}
printf("%4d",(int)data[i]);
}
printf("\n");
return 0;
}