00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <math.h>
00004 
00005 static const double PI=3.14159265358979323;
00006 
00007 /*
00008  * 1-d implemented directly from the formulas.
00009  * Very accurate, very slow.
00010  *
00011  * Modified to compute DCT scaled by sqrt(2)*4
00012  */
00013 
00014 static void dct8_ref(int *data)
00015 {
00016   double output[8] = {0};
00017   short x,n;
00018   for(x=0;x<8;x++) {
00019     for(n=0;n<8;n++)
00020       output[x] += data[n] * cos(PI * x * (2*n+1)/16.0);
00021   }
00022   for(x=0;x<8;x++) {
00023     output[x] /= 4.0;  /* Apply typical weighting to output */
00024     if(x==0) output[x] /= sqrt(2.0);
00025 
00026     output[x] *= 4.0 * sqrt(2);    /* Scale by sqrt(2)*4 */
00027 
00028     data[x] = floor(output[x] + 0.5); /* Round accurately */
00029   }
00030 }
00031 
00032 /*
00033  * 1-d implemented directly from the formulas.
00034  * Very accurate, very slow.
00035  *
00036  * Modified to compute IDCT scaled by sqrt(2)*4
00037  */
00038 
00039 static void idct8_ref(int *data)
00040 {
00041   double output[8] = {0};
00042   short x,n;
00043   for(x=0;x<8;x++) {
00044     output[x]= data[0] / sqrt(2.0);
00045     for(n=1;n<8;n++)
00046       output[x] += data[n] * cos(PI * n * (2*x+1)/16.0);
00047   }
00048   for(x=0;x<8;x++) {
00049     output[x] /= 2.0;
00050     output[x] *= sqrt(2.0);
00051     data[x] = floor(output[x] + 0.5); /* Round accurately */
00052   }
00053 }
00054 
00055 /*--------------------------------------------------------------------------
00056  * From Figure 1 of Loeffler, Ligtenberg, and Moschytz.
00057  * ("Practical Fast 1-D DCT Algorithms with 11 Multiplications,"
00058  * Acoustics, Speech, and Signal Processing, 1989. ICASSP-89, 1989.
00059  * pp 988-991.)
00060  *
00061  * Note: Output is regular DCT scaled by sqrt(2)*4.
00062  *
00063  * The choice of 10-bit fixed-point constants here is arbitrary.
00064  * Using more bits gives better accuracy, but with increased risk
00065  * of overflow.  Note that with 10 bit accuracy, sqrt(2) is 1448/1024,
00066  * which is an exact multiple of 8.  Hence, 181/128 is just as accurate,
00067  * and reduces overflow.
00068  *
00069  * output[x] = sqrt(2)*((x==0)?sqrt(2):1) SUM(n=0..n-1) input[n] * cos(pi*x*(2n+1)/16)
00070  *----------------------------------------------------------------------------*/
00071 
00072 static void dct8(int *dctBlock)
00073 {
00074   static const int c1=1004 ; /* cos(pi/16)<<10 */
00075   static const int s1=200  ; /* sin(pi/16)<<10 */
00076   static const int c3=851  ; /* cos(3pi/16)<<10 */
00077   static const int s3=569  ; /* sin(3pi/16)<<10 */
00078   static const int r2c6=554; /* sqrt(2)*cos(6pi/16)<<10 */
00079   static const int r2s6=1337;
00080   static const int r2=181;   /* sqrt(2)<<7 */
00081   int x0=dctBlock[0], x1=dctBlock[1], x2=dctBlock[2], x3=dctBlock[3],
00082     x4=dctBlock[4], x5=dctBlock[5], x6=dctBlock[6], x7=dctBlock[7];
00083   int x8;
00084 
00085   /* Stage 1 */
00086   x8=x7+x0; x0-=x7;  x7=x1+x6; x1-=x6;
00087   x6=x2+x5; x2-=x5;  x5=x3+x4; x3-=x4;
00088 
00089   /* Stage 2 */
00090   x4=x8+x5; x8-=x5;  x5=x7+x6; x7-=x6;
00091   x6=c1*(x1+x2); x2=(-s1-c1)*x2+x6; x1=(s1-c1)*x1+x6;
00092   x6=c3*(x0+x3); x3=(-s3-c3)*x3+x6; x0=(s3-c3)*x0+x6;
00093 
00094   /* Stage 3 */
00095   x6=x4+x5; x4-=x5;
00096   x5=r2c6*(x7+x8); x7=(-r2s6-r2c6)*x7+x5; x8=(r2s6-r2c6)*x8+x5;
00097   x5=x0+x2;x0-=x2; x2=x3+x1; x3-=x1;
00098 
00099   /* Stage 4, round, and output */
00100   dctBlock[0]=x6;  dctBlock[4]=x4;
00101   dctBlock[2]=(x8+512)>>10; dctBlock[6] = (x7+512)>>10;
00102   dctBlock[7]=(x2-x5+512)>>10; dctBlock[1]=(x2+x5+512)>>10;
00103   dctBlock[3]=(x3*r2+65536)>>17; dctBlock[5]=(x0*r2+65536)>>17;
00104 }
00105 
00106 /***************************************************************************/
00107 
00108 /*
00109  * From Figure 1 of Loeffler, Ligtenberg, and Moschytz.
00110  * ("Practical Fast 1-D DCT Algorithms with 11 Multiplications,"
00111  * Acoustics, Speech, and Signal Processing, 1989. ICASSP-89, 1989.
00112  * pp 988-991.)
00113  *
00114  * Note: Output is regular IDCT scaled by sqrt(2)*4.
00115  *
00116  * output[x] = sqrt(2)*((x==0)?sqrt(2):1) SUM(n=0..n-1) input[n] * cos(pi*x*(2n+1)/16)
00117  */
00118 static void idct8(int *dctBlock)
00119 {
00120   static const int c1=251 ; /* cos(pi/16)<<8 */
00121   static const int s1=50  ; /* sin(pi/16)<<8 */
00122   static const int c3=213 ; /* cos(3pi/16)<<8 */
00123   static const int s3=142 ; /* sin(3pi/16)<<8 */
00124   static const int r2c6=277; /* cos(6pi/16)*sqrt(2)<<9 */
00125   static const int r2s6=669;
00126   static const int r2=181; /* sqrt(2)<<7 */
00127 
00128   /* Stage 4 */
00129   int x0=dctBlock[0]<<9, x1=dctBlock[1]<<7, x2=dctBlock[2],
00130     x3=dctBlock[3]*r2, x4=dctBlock[4]<<9, x5=dctBlock[5]*r2,
00131     x6=dctBlock[6], x7=dctBlock[7]<<7;
00132   int x8=x7+x1; x1 -= x7;
00133 
00134   /* Stage 3 */
00135   x7=x0+x4; x0-=x4; x4=x1+x5; x1-=x5; x5=x3+x8; x8-=x3;
00136   x3=r2c6*(x2+x6);x6=x3+(-r2c6-r2s6)*x6;x2=x3+(-r2c6+r2s6)*x2;
00137 
00138   /* Stage 2 */
00139   x3=x7+x2; x7-=x2; x2=x0+x6; x0-= x6;
00140   x6=c3*(x4+x5);x5=(x6+(-c3-s3)*x5)>>6;x4=(x6+(-c3+s3)*x4)>>6;
00141   x6=c1*(x1+x8);x1=(x6+(-c1-s1)*x1)>>6;x8=(x6+(-c1+s1)*x8)>>6;
00142 
00143   /* Stage 1, rounding and output */
00144   x7+=512; x2+=512;x0+=512;x3+=512;
00145   dctBlock[0]=(x3+x4)>>10;  dctBlock[1]=(x2+x8)>>10;
00146   dctBlock[2]=(x0+x1)>>10;  dctBlock[3]=(x7+x5)>>10;
00147   dctBlock[4]=(x7-x5)>>10;  dctBlock[5]=(x0-x1)>>10;
00148   dctBlock[6]=(x2-x8)>>10;  dctBlock[7]=(x3-x4)>>10;
00149 }
00150 
00151 /*---------------------------------------------------------*/
00152 
00153 #define C1    0.980785  /* cos(PI/16)           */
00154 #define S1    0.19509   /* sin(PI/16)           */
00155 #define C3    0.83147   /* cos(3*PI/16)         */
00156 #define S3    0.55557   /* sin(3*PI/16)         */
00157 #define R2C6  0.541196  /* sqrt(2)*cos(6*PI/16) */
00158 #define R2S6  1.30656   /* sqrt(2)*sin(6*PI/16) */
00159 #define R2    1.41421   /* sqrt(2) (duh)        */
00160 
00161 static void fdct8(float *vv)
00162 {
00163   float x0=vv[0], x1=vv[1], x2=vv[2], x3=vv[3],
00164         x4=vv[4], x5=vv[5], x6=vv[6], x7=vv[7], x8 ;
00165 
00166   x8=x7+x0; x0-=x7;  x7=x1+x6; x1-=x6;
00167   x6=x2+x5; x2-=x5;  x5=x3+x4; x3-=x4;
00168 
00169   x4=x8+x5; x8-=x5;  x5=x7+x6; x7-=x6;
00170   x6=C1*(x1+x2); x2=(-S1-C1)*x2+x6; x1=(S1-C1)*x1+x6;
00171   x6=C3*(x0+x3); x3=(-S3-C3)*x3+x6; x0=(S3-C3)*x0+x6;
00172 
00173   x6=x4+x5; x4-=x5;
00174   x5=R2C6*(x7+x8); x7=(-R2S6-R2C6)*x7+x5; x8=(R2S6-R2C6)*x8+x5;
00175   x5=x0+x2;x0-=x2; x2=x3+x1; x3-=x1;
00176 
00177   vv[0]=x6     ; vv[4]=x4     ;
00178   vv[2]=x8     ; vv[6]=x7     ;
00179   vv[7]=(x2-x5); vv[1]=(x2+x5);
00180   vv[3]=x3*R2  ; vv[5]=x0*R2  ;
00181 }
00182 static void ifdct8(float *vv)
00183 {
00184   float x0=vv[0]   , x1=vv[1], x2=vv[2]   ,
00185         x3=vv[3]*R2, x4=vv[4], x5=vv[5]*R2,
00186         x6=vv[6]   , x7=vv[7], x8=x7+x1    ;
00187 
00188   x1 -= x7;
00189 
00190   x7=x0+x4; x0-=x4; x4=x1+x5; x1-=x5; x5=x3+x8; x8-=x3;
00191   x3=R2C6*(x2+x6);x6=x3+(-R2C6-R2S6)*x6;x2=x3+(-R2C6+R2S6)*x2;
00192 
00193   x3=x7+x2; x7-=x2; x2=x0+x6; x0-= x6;
00194   x6=C3*(x4+x5);x5=(x6+(-C3-S3)*x5);x4=(x6+(-C3+S3)*x4);
00195   x6=C1*(x1+x8);x1=(x6+(-C1-S1)*x1);x8=(x6+(-C1+S1)*x8);
00196 
00197   vv[0]=(x3+x4)/8.0;  vv[1]=(x2+x8)/8.0;
00198   vv[2]=(x0+x1)/8.0;  vv[3]=(x7+x5)/8.0;
00199   vv[4]=(x7-x5)/8.0;  vv[5]=(x0-x1)/8.0;
00200   vv[6]=(x2-x8)/8.0;  vv[7]=(x3-x4)/8.0;
00201 }
00202 
00203 int main( int argc , char * argv[] )
00204 {
00205    int xx[8] = { 1,2,3,4,4,4,3,2 } ;
00206    int yy[8] ; float zz[8] ; int ii ;
00207 
00208    for( ii=0 ;ii < 8 ; ii++ ) zz[ii] = xx[ii] ;
00209 
00210    memcpy(yy,xx,sizeof(int)*8) ;
00211    printf("y = %d %d %d %d %d %d %d %d\n",
00212           yy[0],yy[1],yy[2],yy[3],yy[4],yy[5],yy[6],yy[7]) ;
00213    dct8_ref(yy) ;
00214    printf("dct8_ref = %d %d %d %d %d %d %d %d\n",
00215           yy[0],yy[1],yy[2],yy[3],yy[4],yy[5],yy[6],yy[7]) ;
00216    idct8_ref(yy) ;
00217    printf("idct8_ref = %d %d %d %d %d %d %d %d\n",
00218           yy[0],yy[1],yy[2],yy[3],yy[4],yy[5],yy[6],yy[7]) ;
00219 
00220    memcpy(yy,xx,sizeof(int)*8) ;
00221    dct8(yy) ;
00222    printf("dct8 = %d %d %d %d %d %d %d %d\n",
00223           yy[0],yy[1],yy[2],yy[3],yy[4],yy[5],yy[6],yy[7]) ;
00224    idct8(yy) ;
00225    printf("idct8 = %d %d %d %d %d %d %d %d\n",
00226           yy[0],yy[1],yy[2],yy[3],yy[4],yy[5],yy[6],yy[7]) ;
00227 
00228    fdct8(zz) ;
00229    printf("fdct8 = %g %g %g %g %g %g %g %g\n",
00230           zz[0],zz[1],zz[2],zz[3],zz[4],zz[5],zz[6],zz[7]) ;
00231    ifdct8(zz) ;
00232    printf("ifdct8 = %g %g %g %g %g %g %g %g\n",
00233           zz[0],zz[1],zz[2],zz[3],zz[4],zz[5],zz[6],zz[7]) ;
00234 
00235    exit(0) ;
00236 }