tc.y

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

	      $$.dim[1] = max($1.dim[1], $3.dim[1]);

	      for (i=0; i<min($1.dim[0], $3.dim[0]); i++) {
		for (j=0; j<min($1.dim[1], $3.dim[1]); j++)
		  $$.v[i * $$.dim[0] + j] = $1.v[i * $1.dim[0] + j] * 
		    	$3.v[ i* $3.dim[0] + j];
		for (j; j<max($1.dim[1], $3.dim[1]); j++)
		  $$.v[i * $$.dim[0] + j] = 0;
	      }
	      for (; i<max($1.dim[0], $3.dim[0]); i++) {
		for (j=0; j<max($1.dim[1], $3.dim[1]); j++)
		  $$.v[i* $$.dim[0] + j] = 0;
	      }
	    }
	| VAR
	    {
	      $$ = varnames[$1];
	    }
	| ICON RANDOM ICON
	    {
	      int i;
	      $$.dim[0] = $$.dim[1] = 8;
	      for(i=0; i<64; i++)
		$$.v[i] = rand_dct(-$2, $3) ;
	    }
	| '{' scalar_expr_list '}'		{ $$ = $2; }
	;

scalar_expr_list:
	  fcon				{  $$.dim[0] = 1; $$.dim[1] = 1; $$.v[0] = $1; }
	| scalar_expr_list fcon		{ $$ = $1; $$.v[$$.dim[1]++] = $2; }
	;					  

o_mode:
	  _4POINT			{ $$ = O_4POINT; }
	| _8POINT			{ $$ = O_8POINT; }
	| SEP				{ $$ = O_SEP; }
	|				{ $$ = 0; }
	;
					  
o_fwd:
	  INV				{ $$ = 0; }
	|				{ $$ = O_FWD; }
	;
					  

	


fcon:
	  INF
	     {
	       $$ = 100000000;
	     }
	| FCON	{ $$ = $1; }
	| ICON	{ $$ = $1; }
         ;

%%

#include <math.h>

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;
   }
}

#define	THETA	  (M_PI/8)		/* 22.5 degrees */ 	
#define SIN_THETA (sqrt(2)*sin(M_PI/8))	/* 0.2706*2   */
#define COS_THETA (sqrt(2)*cos(M_PI/8))	/* 0.6533*2   */
void
dct4(double *tabval, double *ptres, int step, int dim, int fwd)
{
     double x0, x1, x2, x3, y0, y1, y2, y3;

     if (dim != 4)
       fprintf(stderr, "dct requires four points\n");

     /* stage 1 */
     x0 = tabval[0] + tabval[3];
     x1 = tabval[1] + tabval[2];
     x2 = tabval[1] - tabval[2];
     x3 = tabval[0] - tabval[3];

     /* stage 2 */
     y0 = x0 + x1;
     y1 = x0 - x1;
/*
     y2 =  SIN_THETA * x2 + COS_THETA * x3;
     y3 = -COS_THETA * x2 + SIN_THETA * x3; 
*/
     y2 =  (COS_THETA-SIN_THETA)*x3 + SIN_THETA*(x2+x3);
     y3 = -(SIN_THETA+COS_THETA)*x2 + SIN_THETA*(x2+x3);

     ptres[0*step] = y0 ;
     ptres[2*step] = y1 ;
     ptres[1*step] = y2;
     ptres[3*step] = y3;
}


#undef THETA
#undef SIN_THETA
#undef COS_THETA

#define	THETA	  (M_PI/16)		/* 11.25 degrees */ 	
#define SIN_THETA (sqrt(2)*sin(2*THETA))
#define COS_THETA (sqrt(2)*cos(2*THETA))

#define	COEFF_a	sqrt(2.)*(-cos(1*THETA)+cos(3*THETA)+cos(5*THETA)-cos(7*THETA))
#define	COEFF_b	sqrt(2.)*( cos(1*THETA)+cos(3*THETA)-cos(5*THETA)+cos(7*THETA))
#define	COEFF_c	sqrt(2.)*( cos(1*THETA)+cos(3*THETA)+cos(5*THETA)-cos(7*THETA))
#define	COEFF_d	sqrt(2.)*( cos(1*THETA)+cos(3*THETA)-cos(5*THETA)-cos(7*THETA))
#define	COEFF_e	sqrt(2.)*(-cos(3*THETA)+cos(7*THETA))  
#define	COEFF_f	sqrt(2.)*(-cos(1*THETA)-cos(3*THETA))
#define	COEFF_g	sqrt(2.)*(-cos(3*THETA)-cos(5*THETA))
#define	COEFF_h	sqrt(2.)*(-cos(3*THETA)+cos(5*THETA))
#define	COEFF_i	sqrt(2.)*( cos(3*THETA))

void
floatdct8(double *tab, double *res, int step, int dim, int fwd)
{
     double x0, x1, x2, x3, x4, x5, x6, x7;
     double y0, y1, y2, y3;
     double z1, z2, z3, z4, z5;
     double a, b, c, d, e, f, g, h, i;
     double o0, o1, o2, o3, o4, o5, o6, o7; 

     if (dim != 8 || fwd!=1)
       fprintf(stderr, "dct must be forward, requires eight points\n");

     /* stage 1, common */
     x0 = tab[0] + tab[7]; x1 = tab[1] + tab[6];
     x2 = tab[2] + tab[5]; x3 = tab[3] + tab[4];
     x4 = tab[3] - tab[4]; x5 = tab[2] - tab[5];
     x6 = tab[1] - tab[6]; x7 = tab[0] - tab[7];

     /* stages 2 and 3, even part */
     y0 = x0 + x3;  y1 = x1 + x2; y2 = x1 - x2;  y3 = x0 - x3;
     o0 = y0 + y1; o2 =  (COS_THETA-SIN_THETA)*y3 + SIN_THETA*(y2+y3);
     o4 = y0 - y1; o6 = -(SIN_THETA+COS_THETA)*y2 + SIN_THETA*(y2+y3);

     /* stages 2 and 3, odd part */
     z1 = x4 + x7;   z2 = x5 + x6;  z3 = x4 + x6;
     z4 = x5 + x7;   z5 = (x4 + x6) + (x5 + x7);
     a = x4*COEFF_a;  b = x5*COEFF_b;  c = x6*COEFF_c;
     d = x7*COEFF_d;  e = z1*COEFF_e;  f = z2*COEFF_f;
     g = z3*COEFF_g;  h = z4*COEFF_h;  i = z5*COEFF_i;

     /* stage 4, odd part */
     o1 = d + e + (h + i);  o3 = c + f + (g + i);
     o5 = b + f + (h + i);  o7 = a + e + (g + i);

     res[0*step] = o0; res[2*step] = o2;  
     res[1*step] = o1; res[3*step] = o3;
     res[4*step] = o4; res[6*step] = o6;
     res[5*step] = o5; res[7*step] = o7;
}


#define	HROUND(x)	(((x) + 0x4000) >> 15)

#define	C0	0xa73b4b42
#define	C1	0x11a8cdb7
#define	C2	0xcdb7a73b
#define	C3	0x4b42ee58
#define	C4	0x4b4211a8
#define	C5	0x3249a73b
#define	C6	0x11a83249
#define	C7	0x58c54b42

#define C8	0x40004000
#define C9	0x4000c000
#define C10	0x539e22a3
#define C11     0x22a3ac62



#define PACKWORD(a,b) PACK16LSB(((int)(a)),((int)b))


void
fixdct8(double *tab, double *res, int step, int dim, int fwd)
{
	int s0, s1, s2, s3, s4, s5, s6, s7;
	int tmp0, tmp1, tmp2, tmp3, tmp101, tmp132, tmp176, tmp145, tmp201, tmp232; 

	tmp0 =  PACKWORD(tab[0], tab[1]);
	tmp1 =  PACKWORD(tab[3], tab[2]);
	tmp2 =  PACKWORD(tab[4], tab[5]);
	tmp3 =  PACKWORD(tab[7], tab[6]);
 				
	tmp101 = DSPIDUALADD(tmp0,tmp3);
	tmp132 = DSPIDUALADD(tmp1,tmp2);
	tmp176 = DSPIDUALSUB(tmp0,tmp3);
	tmp145 = DSPIDUALSUB(tmp1,tmp2);


	tmp101 = DSPIDUALADD(tmp0,tmp3);
	tmp132 = DSPIDUALADD(tmp1,tmp2);
	tmp176 = DSPIDUALSUB(tmp0,tmp3);
	tmp145 = DSPIDUALSUB(tmp1,tmp2);
				
	tmp201 = DSPIDUALADD(tmp101,tmp132);
	tmp232 = DSPIDUALSUB(tmp101,tmp132);
				
	s0 = IFIR16(tmp201,C8);
	s4 = IFIR16(tmp201,C9);
	s2 = IFIR16(tmp232,C10);
	s6 = IFIR16(tmp232,C11);

	s7 = IFIR16(tmp145,C0) + IFIR16(tmp176,C1);
	s3 = IFIR16(tmp145,C2) + IFIR16(tmp176,C3);
	s5 = IFIR16(tmp145,C4) + IFIR16(tmp176,C5);
	s1 = IFIR16(tmp145,C6) + IFIR16(tmp176,C7);

	res[0] = HROUND(s0);
	res[1] = HROUND(s1);	
	res[2] = HROUND(s2);
	res[3] = HROUND(s3);	
	res[4] = HROUND(s4);
	res[5] = HROUND(s5);	
	res[6] = HROUND(s6);
	res[7] = HROUND(s7);
}




void
fdct2d(double *tabval, double *ptres, int dim1, int dim2, int fwd)
 {
 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<dim2;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]*facu[v!=0];
         sum += value * coef;
       }
     ptres[(u*dim1)+v] =  sum ;
     }
   }
}

void
idct2d(double *tabval, double *ptres, int dim1, int dim2, int fwd)
 {		
 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 ;
   }
 }
}


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

  if (state->onepass) {
    for (u=0;u<dim1;u++)
      dct1d(&tabval[dim2*u], &ptres[dim2*u], 1, dim2, fwd);
  } else {
    for (u=0;u<dim1;u++)
      dct1d(&tabval[dim2*u], &tabtemp[dim2*u], 1, dim2, fwd);
    memset(ptres, 0, MAXVAL*sizeof(double));
    
    for (v=0; v<dim2; v++)
      dct1d(&tabtemp[v], &ptres[v], dim2, dim1, fwd);
  }
}

void printfloat(double *block, char *lab, int dim1, int dim2, FILE *fp)
{
  int i, j;
  double val;

  if (lab)
    fprintf(fp, "\n%s\n\n", lab);

  for (i=0; i<dim1; i++) {
    fprintf(fp, "\t");	
    for (j=0; j<dim2; j++) {
      val = block[(i*dim2) + j];
      if (fabs(val) < 0.0000001)
	val = 0;
      fprintf(fp, "%c%-*.*f%s", val<0 ? '-' : ' ', j==dim2-1 ? 0 : state->width, state->scale, fabs(val), j==dim2-1 ? "" : " ");
    }
    fprintf(fp, "\n");
  }
}

const	int	I_j[] = {0, 2, 4, 6, 1, 3, 5, 7};
void shuffle(double *swap, double *block)
{
  int i, j;

	/* transpose input block with shuffling and scaling up by factor 8 */
	/* for accomodating funny Trimedia format for input of IDCT : note */
	/* that this is done costless in normal code by iquant/zz.	   */

	/* endianess has been taken care of by definition of H_ROUNDING */
	for	(j = 0; j < 8; j++)
	{
		for	(i = 0; i < 8; i++)
		{
			swap[(i * 8) | I_j[j]] = *block++ * 8;
		}
	}
}

myround(double *tabval, double *ptres, double fact)
{ 
  int i;

  for (i=0; i<64; i++)
	ptres[i] = tcsirint(tabval[i]/fact) * fact;
}

pause(n) { }


blockmin(double *tabval, double *ptres, double val)
{
  int i;
  for (i=0; i<64; i++) {
    ptres[i] = tabval[i];
    if (ptres[i] < val)
      ptres[i] = val;
  }
}

blockmax(double *tabval, double *ptres, double val)
{
  int i;
  for (i=0; i<64; i++) {
    ptres[i] = tabval[i];
    if (ptres[i] > val)
      ptres[i] = val;
  }
}

void not_a_dct() { }

reflect(double *ptres, int dim)
{
  double swap;
  int i;

  for (i=0; i<dim/2; i++) {
    swap = ptres[i];
    ptres[i] = ptres[dim-i-1] ;
    ptres[dim-i-1] = swap;
  }
}


void fastshuffle(double *block, short *swap)
{
  int i, j;

	/* transpose input block with shuffling and scaling up by factor 8 */
	/* for accomodating funny Trimedia format for input of IDCT : note */
	/* that this is done costless in normal code by iquant/zz.	   */

	/* endianess has been taken care of by definition of H_ROUNDING */
	for	(j = 0; j < 8; j++)
	{
		for	(i = 0; i < 8; i++)
		{
			swap[(i << 3) | I_j[j]] = *block++ * 8;
		}
	}
}

void
fast_8x8dct(double *tabval, double *ptres, int dim1, int dim2, int fwd)
{
  int i;
  short tabtemp[64], tabres[64];

  for (i=0; i<64; i++)
    tabtemp[i] = tabval[i];

  dct8x8fix (tabtemp, tabres);

  for (i=0; i<64; i++)
    ptres[i] = tabres[i];
}

void	fast_8x8idct(double *tabval,  double *ptres, int dim1, int dim2, int fwd)
{
  short	datain [64];
  short dataout[64];
  int i;

  for (i=0; i<64; i++)
    datain[i] = tabval[i];

  idct8x8fix(datain, dataout, state->clip);

  for (i = 0; i<64; i++)
    ptres[i] = dataout[i];
}


char *
toname(int n)
{
  if (n==1) return "positive";
  if (n==-1) return "negative";
  return "**gok**";
}