|
|
может функция громоздская, но хоть какие-то осмысленные действия делает
гарантия отсутствия переполнений входными данными
const short NS=2, IS=0; // >>NS+IS == /4
const int N=8;inline void ldct(short *x)
{
short x1[N], x2[N], t;const short p[]={
//13/32 11/16 3/16 13/32 13/32
13, 22, 6, 13, 13
};
const short u[]={
//11/32 15/32 3/16 11/16
11, 15, 6, 22
};
const short SH=5;// in DSP u,p[i]*2048(0x10000/32) use R.H(A*B)
//first fly
x1[0]=x[0]+x[7];
x1[1]=x[1]+x[6];
x1[2]=x[2]+x[5];
x1[3]=x[3]+x[4];x1[4]=x[3]-x[4];
x1[5]=x[2]-x[5];
x1[6]=x[1]-x[6];
x1[7]=x[0]-x[7];//first lift
t=p[3]*x1[6]>>SH;
x1[5]=x1[5]-t;
t=x1[5]*u[3]>>SH;
x1[6]=t+x1[6];
t=p[4]*x1[6]>>SH;
x1[5]=-x1[5]+t;//second bfly
x2[0]=x1[0]+x1[3];
x2[3]=x1[0]-x1[3];x2[1]=x1[1]+x1[2];
x2[2]=x1[1]-x1[2];x2[4]=x1[4]+x1[5];
x2[5]=x1[4]-x1[5];x2[6]=-x1[6]+x1[7];
x2[7]=x1[6]+x1[7];//second lift
x2[0]=x2[0]+x2[1];
x2[1]=-x2[1]+(x2[0]>>1);t=p[0]*x2[3]>>SH;
x2[2]=-x2[2]+t;
t=u[0]*x2[2]>>SH;
x2[3]=x2[3]-t;t=p[2]*x2[7]>>SH;
x2[4]=-x2[4]+t;
t=u[2]*x2[4]>>SH;
x2[7]=x2[7]-t;t=p[1]*x2[6]>>SH;
x2[5]=x2[5]+t;
t=u[1]*x2[5]>>SH;
x2[6]=x2[6]-t;//mult sin/cos in Q table
//output perm and /????
x[0]=x2[0]>>NS;
x[4]=x2[1]>>NS;
x[6]=x2[2]>>NS;
x[2]=x2[3]>>NS;
x[7]=x2[4]>>NS;
x[5]=x2[5]>>NS;
x[3]=x2[6]>>NS;
x[1]=x2[7]>>NS;
}
E-mail: 
info@telesys.ru
