DCT algorithm can be said to be the first step of lossy compression, more used in video compression.

Transfer from http://m.blog.csdn.net/blog/zhu530548851/21444459

DCT algorithm can be said to be the first step of lossy compression, more used in video compression. It transforms the two dimensions into one-dimensional data and aggregates the energy into the upper-left corner. About DCT, there are many fast algorithms, I personally think that my algorithm is still relatively fast. Here's My Code:

#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;iC[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{
Temp1=f[wing*winglen+i_halfwing];
temp2=f[(wing+1) *winglen-1-i_halfwing];
if (wing%2)
Swap (&AMP;TEMP1,&AMP;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{//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&LT;&LT;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;iC[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{
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&AMP;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;
}