fdct.c

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

/*
 * Copyright (c) 1998, 2000 TriMedia Technologies Inc.          
 *
 * +------------------------------------------------------------------+
 * | This software is furnished under a license and may only be used  |
 * | and copied in accordance with the terms and conditions of  such  |
 * | a license and with the inclusion of this copyright notice. This  |
 * | software or any other copies of this software may not be provided|
 * | or otherwise made available to any other person.  The ownership  |
 * | and title of this software is not transferred.                   |
 * |                                                                  |
 * | The information in this software is subject  to change without   |
 * | any  prior notice and should not be construed as a commitment by |
 * | TriMedia Technologies .                                          |
 * |                                                                  |
 * | this code and information is provided "as is" without any        |
 * | warranty of any kind, either expressed or implied, including but |
 * | not limited to the implied warranties of merchantability and/or  |
 * | fitness for any particular purpose.                              |
 * +------------------------------------------------------------------+
 *
 *  Module name              : fdct.c    1.2
 *
 *  Last update              : 18:31:00 - 2000/05/12
 *
 *  Description              : forward 2-D DCT reference implementation
 *
 * The reference implementation is provided along with the optimized TriMedia
 * implementation of the two dimensional forward Discrete Cosine Transform
 * (DCT). It can be used to validating the exactness of the
 * results or to compare with a customer implementations. Both separable (n^3)
 * and standard (n^4) algorithms are provided.
 *
 * Compile :       tmcc fdct.c -lm -o fdct.out
 *     
 * Run :           tmsim fdct.out [-sep] [-onepass] [-scale factor] <input>
 *
 * <input> : 	   ASCII file, 64 integer values, free format
 *
 * Options:	   -sep		Use separable algorithm (one pass)
 *		-onepass	Only do first pass. Implies -sep
 *		-scale		Multiply outputs by scale factor, for
 *				1 -D test only at this point.	 
 *				Use -scale 2.82856 ( 2 \/2  ) to
 *				obtain the same results as custom ops.
 *

 */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define     M_PI            3.14159265358979323846
#define     M_SQRT1_2       0.70710678118654752440

#define	ONEPASS	1
#define	SEP	2


void fdct2d  (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);
void debugprint(double *, int, int);

double	data[64] ;
double	result[64];
double scale = 1.0;
int flags;

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

  while (argc > 1) {
    if (!strcmp(argv[1], "-scale")) {
      sscanf(argv[2], "%lf", &scale);
      argc--;
      argv++;
    } else if (!strcmp(argv[1], "-onepass"))
      flags = ONEPASS|SEP;
    else if (!strcmp(argv[1], "-sep"))
      flags = SEP;
    else
      break;

    argc--;
    argv++;
  }

  if (!freopen(argv[1], "r", stdin)) {
    fprintf(stderr, "Cannot open input file %s\n", argv[1]);
    exit(1);
  }

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

  if (flags&SEP)
    dct2dsep(data, result, 8, 8, 1);
  else
    fdct2d(data, result, 8, 8);

  for (i=0; i<8; i++) {				/* print output */
    for (j=0; j<8; j++)
      printf("%6d ", (int)result[i*8+j]);
    printf("\n");
  }
  exit(EXIT_SUCCESS);
}

/*
 * 1-D Forward/Reverse Discrete Cosine Transform (standard algorithm)
 */

void
dct1d(double *tabval, double *ptres, int step, int dim, int fwd)
 {
 int 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;
   }
}

/* 
 * 2-D Discrete Cosine Transform (standard algorithm)
 */
void
fdct2d(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) * M_SQRT1_2;
  facu[1] = sqrt(2./dim1);
  facv[0] = sqrt(2./dim2) * M_SQRT1_2;
  facv[1] = sqrt(2./dim2);

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

/*
 * 2-D Forward/Reverse Discrete Cosine Transform (separable algorithm)
 * Uses 1-D algorithm
 */
#define MAXVAL 1000

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

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

  if (flags&ONEPASS) {			/* print intermediate results */
    for (v=0; v<dim2; v++)
      debugprint (&tabtemp[v], dim2, dim1);
    exit(0);
  }

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


/*
 * print out intermediate results here
 */ 
void debugprint(double *ptr, int step, int dim)
{
  int i;

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