fidct.c

用于TM1300/PNX1300系列DSP（主要用于视频处理）的各种滤波器源码

#define     M_PI            3.14159265358979323846
#define     M_SQRT1_2       0.70710678118654752440
#include <math.h>
#include <stdio.h>

double	data[64] ;
double	result[64];

void idct2d  (double *tabval, double *ptres, int dim1, int dim2);
void dct1d   (double *tabval, double *ptres, int step, int dim, int fwd);
void dct2dsep(double *tabval, double *ptres, int dim1, int dim2, int fwd);


int idct_debug = 0;	
int idct_sep   = 0;
double scale = 1.0;
int clip_val = 10000;

double clip(double a)
{
   if (a < -clip_val)
     return -clip_val;
   if (a >  clip_val)
     return  clip_val;
   return a;
}

main(int argc, char **argv)
{
  int i, j, dataword;
  double atof();

  while (argc > 1) {
    if (!strcmp(argv[1], "-scale")) {
      scale = atof(argv[2]);
      argc--;
      argv++;
    } else if (!strcmp(argv[1], "-onepass"))
      idct_debug = 1;
    else if (!strcmp(argv[1], "-sep"))
      idct_sep = 1;
    else if (!strcmp(argv[1], "-clip")) {
      clip_val = atoi(argv[2]);
      argc--;
      argv++;
    }

    argc--;
    argv++;
  }

  for (i=0; i<64; i++) {
    scanf("%d", &dataword);
    data[i] = dataword;
  }

  if (idct_sep)
    dct2dsep(data, result, 8, 8, 0);
  else
    idct2d(data, result, 8, 8);

  for (i=0; i<8; i++) {
    for (j=0; j<8; j++)
      printf("%5.0f ", clip(result[8*i+j] * scale));
    printf("\n");
  }
  exit(0);
}

void
dct1d(double *tabval, double *ptres, int step, int dim, int fwd)
 {
 register v, i;
 double sum, *ptval, coef;
 double facmul[2] ;

 facmul[0] = sqrt(2./dim) * M_SQRT1_2;
 facmul[1] = sqrt(2./dim);

 for (v=0;v<dim;v++) {
     ptval = tabval ;
     sum = 0;
     for (i=0;i<dim;i++) {
       if (fwd)
	 coef  = cos(((2*i+1)*v*M_PI)/(2*dim)) * facmul[v!=0];
       else 
         coef  = cos(((2*v+1)*i*M_PI)/(2*dim)) * facmul[i!=0];
       sum +=  *ptval * coef ; 
       ptval += step;
     }
     *ptres =  sum;
     ptres +=  step;
   }
}

void idctprint(double *ptr, int step)
{
  int i;

  for (i=0; i<8; i++)
    printf("%-5.0f ", clip(ptr[step*i] * scale));
  printf("\n");
}

#define MAXVAL 1000

void
dct2dsep(double *tabval, double *ptres, int dim1, int dim2, int fwd)
{
  int u,v;
  double tabtemp[MAXVAL]; 	

  memset(ptres, 0, MAXVAL*sizeof(double));

  for (u=0;u<dim1;u++)
    dct1d(&tabval[u], &tabtemp[u], dim2, dim1, fwd);

  if (idct_debug) {
    for (v=0; v<dim1; v++)
      idctprint (&tabtemp[dim2*v], 1);
    exit(0);
  }

  for (v=0; v<dim1; v++)
    dct1d(&tabtemp[dim2*v], &ptres[dim2*v], 1, dim2, fwd);
}




void
idct2d(double *tabval, double *ptres, int dim1, int dim2)
{		
 int u,v,i,j ;
 double facu[2] , facv[2], sum, value, coef;

 facu[0] = sqrt(2./dim1);
 facu[1] = sqrt(2./dim1) * M_SQRT1_2;
 facv[0] = sqrt(2./dim2);
 facv[1] = sqrt(2./dim2) * M_SQRT1_2;

 for (i=0;i<dim1;i++) {
   for (j=0;j<dim2;j++) {
     sum = 0;  
     for (u=0;u<dim1;u++) {
       for (v=0;v<dim2;v++) {
	 value = tabval[(u*dim1) + v];
         coef  = cos(((2*i+1)*u*M_PI)/(2*dim1))* cos(((2*j+1)*v*M_PI)/(2*dim2)) *
	   facu[!u]*facu[!v];
         sum += value * coef;
       }
     }
     ptres[(i*dim1)+j] =  sum ;
   }
 }
}