dct.c

DSP实现FDCT算法

#include <math.h>
#define N 8
#define PI 3.1415926
float DCTinput[N]={1.0,2,3,4,5,6,7,8},FDCTinput[N]={1.0,2,3,4,5,6,7,8};
float DCTout[N],IDCTout[N];
float x[N];
float X[N];
float sp[N];
float xp, Xp;
#define  C1 0.9808
#define  C2 0.9239
#define  C3 0.8315
#define  C4 0.7071
#define  C5 0.5556
#define  C6 0.3827
#define  C7 0.1951
float blk[N];
float blk1[N];
float blk2[N];
	

dct(float s[])
{
	int n, k;
	float alpha;
	float pin, sqrtn, sqrt2n;
	
	pin    = 0.5*PI/N;
	sqrtn  = sqrt(1.0/N);
	sqrt2n = sqrt(2.0/N);	

	for(k=0; k<N; k++)
		{
			    if(k==0)
				alpha=sqrtn;
				else
					alpha=sqrt2n;
			for(n=0;n<N;n++){
				x[n]=s[n];
				     }
			xp=0.0;
			for(n=0; n<N; n++){
				xp+=x[n]*cos((2*n+1)*k*pin);
			}
			xp*=alpha;
		DCTout[k]=xp;
		}
	
  
 } 
 idct(float c[])
    {
	int n, k;
	float alpha;
	float pin, sqrtn, sqrt2n;
	
	pin    = 0.5*PI/N;
	sqrtn  = sqrt(1.0/N);
	sqrt2n = sqrt(2.0/N);	
    for(n=0; n<N; n++)
		{
          
		  Xp=0.0;
		  for(k=0;k<N;k++)
		  {
               if(k==0)
				alpha=sqrtn;
			  	else
				alpha=sqrt2n;
				Xp+=alpha*DCTout[k]*cos((2*n+1)*k*pin);
		  }
		  IDCTout[n]=Xp;
		}
		
	}
	



fdct(float *blk)//先做行DCT
  { 
    float S07,S16,S25,S34,S0734,S1625;
    float D07,D16,D25,D34,D0734,D1625;
    S07=blk[0]+blk[7];
    S16=blk[1]+blk[6];
    S25=blk[2]+blk[5];
    S34=blk[3]+blk[4];
    S0734=S07+S34;
    S1625=S16+S25;
    D07=blk[0]-blk[7]; 
    D16=blk[1]-blk[6];
    D25=blk[2]-blk[5];
    D34=blk[3]-blk[4];
    D0734=S07-S34;
    D1625=S16-S25;
    blk1[0]=0.5*(C4*(S0734+S1625));
    blk1[1]=0.5*(C1*D07+C3*D16+C5*D25+C7*D34);
    blk1[2]=0.5*(C2*D0734+C6*D1625);
    blk1[3]=0.5*(C3*D07-C7*D16-C1*D25-C5*D34);
    blk1[4]=0.5*(C4*(S0734-S1625));
    blk1[5]=0.5*(C5*D07-C1*D16+C7*D25+C3*D34);
    blk1[6]=0.5*(C6*D0734-C2*D1625);
    blk1[7]=0.5*(C7*D07-C5*D16+C3*D25-C1*D34);
  }


  
fidct (float *blk)
{ double tmp[16]; 

 tmp[0]=blk1[0]*C4+blk1[2]*C2; 
 tmp[1]=blk1[4]*C4+blk1[6]*C6; 
 tmp[2]=blk1[0]*C4+blk1[2]*C6; 
 tmp[3]=-blk1[4]*C4-blk1[6]*C2; 
 tmp[4]=blk1[0]*C4-blk1[2]*C6; 
 tmp[5]=-blk1[4]*C4+blk1[6]*C2; 
 tmp[6]=blk1[0]*C4-blk1[2]*C2; 
 tmp[7]=blk1[4]*C4-blk1[6]*C6; 
 tmp[8]=blk1[1]*C7-blk1[3]*C5; 
 tmp[9]=blk1[5]*C3-blk1[7]*C1; 
 tmp[10]=blk1[1]*C5-blk1[3]*C1; 
 tmp[11]=blk1[5]*C7+blk1[7]*C3; 
 tmp[12]=blk1[1]*C3-blk1[3]*C7; 
 tmp[13]=-blk1[5]*C1-blk1[7]*C5; 
 tmp[14]=blk1[1]*C1+blk1[3]*C3; 
 tmp[15]=blk1[5]*C5+blk1[7]*C7; 
 //second step 
 tmp[0]=0.5*(tmp[0]+tmp[1]); 
 tmp[1]=0.5*(tmp[2]+tmp[3]); 
 tmp[2]=0.5*(tmp[4]+tmp[5]); 
 tmp[3]=0.5*(tmp[6]+tmp[7]); 
 tmp[4]=0.5*(tmp[8]+tmp[9]); 
 tmp[5]=0.5*(tmp[10]+tmp[11]); 
 tmp[6]=0.5*(tmp[12]+tmp[13]); 
 tmp[7]=0.5*(tmp[14]+tmp[15]); 
 //third step 
 blk2[0]=tmp[0]+tmp[7]; 
 blk2[1]=tmp[1]+tmp[6]; 
 blk2[2]=tmp[2]+tmp[5]; 
 blk2[3]=tmp[3]+tmp[4]; 
 blk2[4]=tmp[3]-tmp[4]; 
 blk2[5]=tmp[2]-tmp[5]; 
 blk2[6]=tmp[1]-tmp[6]; 
 blk2[7]=tmp[0]-tmp[7]; 
 
 }



main()
{
    dct(DCTinput);
    idct(DCTout);
    fdct(FDCTinput);
    fidct(blk1);
    while(1);
 }