www.pudn.com > AVS_M_ver10.rar > fft3.c
/*
***********************************************************************
* COPYRIGHT AND WARRANTY INFORMATION
*
* Copyright 2007 Audio Video Coding Standard, Part ¢ú
*
* This software module was developed by AVS Audio sub-group
*
* DISCLAIMER OF WARRANTY
*
* These software programs are available to the users without any
* license fee or royalty on an "as is" basis. The AVS disclaims
* any and all warranties, whether express, implied, or statutory,
* including any implied warranties of merchantability or of fitness
* for a particular purpose. In no event shall the contributors or
* the AVS be liable for any incidental, punitive, or consequential
* damages of any kind whatsoever arising from the use of this program.
*
* This disclaimer of warranty extends to the user of this program
* and user's customers, employees, agents, transferees, successors,
* and assigns.
*
* The AVS does not represent or warrant that the program furnished
* hereunder are free of infringement of any third-party patents.
* Commercial implementations of AVS, including shareware, may be
* subject to royalty fees to patent holders. Information regarding
* the AVS patent policy is available from the AVS Web site at
* http://www.avs.org.cn
*
* THIS IS NOT A GRANT OF PATENT RIGHTS - SEE THE AVS PATENT POLICY.
************************************************************************
*/
#include
#include
#include
#include "../include/amr_plus.h"
/*_____________________________________________________________________
| |
| Constant Definitions |
|_____________________________________________________________________|
*/
#define ORDER_MAX 7
void fft3(float X[], float Y[], short n)
{
float Z[N_MAX/6];
float *Z0, *Z1, *Z2;
float *z0, *z1, *z2;
float *x;
float *yre, *yim, *zre, *zim, *wre, *wim;
short m, step, sign, order;
short i, j, k;
short Ind[3];
short *ind;
/* Determine the order of the transform, the length of decimated */
/* transforms m, and the step for the sine and cosine tables. */
switch(n) {
case 48:
order = 4;
m = 16;
step = 24;
break;
case 96:
order = 5;
m = 32;
step = 12;
break;
case 192:
order = 6;
m = 64;
step = 6;
break;
default:
printf(" invalid fft3 size!\n");
exit(0);
}
/* Compose decimated sequences X[3i], X[3i+1],X[3i+2] */
/* compute their FFT of length m. */
Z0 = &Z[0]; z0 = &Z0[0];
Z1 = &Z0[m]; z1 = &Z1[0]; /* Z1 = &Z[ m]; */
Z2 = &Z1[m]; z2 = &Z2[0]; /* Z2 = &Z[2m]; */
x = &X[0];
for (i = 0; i < n/3; i++)
{
*z0++ = *x++; /* Z0[i] = X[3i]; */
*z1++ = *x++; /* Z1[i] = X[3i+1]; */
*z2++ = *x++; /* Z2[i] = X[3i+2]; */
}
fft_rel(&Z0[0], m, order);
fft_rel(&Z1[0], m, order);
fft_rel(&Z2[0], m, order);
/* Compute the DC coefficient and store it into Y[0]. Note that */
/* the variables Z0, ..., Z8, z0, ..., z8 are not needed after */
/* this. */
*Y = *Z0 + *Z1 + *Z2;
/* Initialize the index table, which points to the sine and */
/* cosine tables. */
ind = &Ind[0];
for (k = 1; k < 3; k++) {
*ind++ = (short) (k*step);
}
/* Butterflies of order 3. */
sign = 1;
zre = &Z[1]; zim = &Z[m-1];
yre = &Y[1]; yim = &Y[n/2+1];
for (i = 0; i < 3; i++) {
for (j = 1; j < m/2; j++) {
wre = &zre[0];
wim = &zim[0];
ind = &Ind[0];
*yre = *wre;
*yim = sign*(*wim);
for (k = 1; k < 3; k++) {
wre += m;
wim += m;
*yre += (*wre)*t_cos[*ind] + sign*(*wim)*t_sin[*ind];
*yim += -(*wre)*t_sin[*ind] + sign*(*wim)*t_cos[*ind];
*ind = (short) (*ind+k*step);
if (*ind >= N_MAX)
{
*ind -= N_MAX;
}
ind++;
}
yre++; zre += sign;
yim++; zim -= sign;
}
wre = &zre[0];
ind = &Ind[0];
*yre = *wre;
/* *yim = 0.0; */
if (i<2)
{
*yim = 0.0;
}
for (k = 1; k < 3; k++) {
wre += m;
*yre += (*wre)*t_cos[*ind];
/* *yim += -(*wre)*t_sin[*ind]; */
if (i<2)
{
*yim += -(*wre)*t_sin[*ind];
}
*ind = (short)(*ind+k*step);
if (*ind >= N_MAX)
{
*ind -= N_MAX;
}
ind++;
}
sign = (short) -sign;
yre++; zre += sign;
yim++; zim -= sign;
}
for (i=0; i= N_MAX)
{
*ind -= N_MAX;
}
ind++;
}
z = &Z[k*m];
*z = *yr0; ind = &Ind[1];
*z += *yr1*t_cos[*ind] - *yi1*t_sin[*ind]; ind++;
*z += *yr1*t_cos[*ind] + *yi1*t_sin[*ind];
zre = &z[1]; zim = &z[m-1];
yr0f = &yr0[ 1]; yi0f = &yi0[0];
yr1f = &yr1[ 1]; yi1f = &yi0f[m];
yr0b = &yr1[-1]; yi0b = &yi1[-1];
for (i = 1; i < m/2; i++) {
ind = &Ind[0];
for (j = 0; j < 3; j++) {
*ind = (short) (*ind+step);
if (*ind >= N_MAX)
{
*ind -= N_MAX;
}
ind++;
}
ind = &Ind[0];
*zre = *yr0f*t_cos[*ind] - *yi0f*t_sin[*ind];
*zim = *yr0f*t_sin[*ind] + *yi0f*t_cos[*ind]; ind++;
*zre += *yr1f*t_cos[*ind] - *yi1f*t_sin[*ind];
*zim += *yr1f*t_sin[*ind] + *yi1f*t_cos[*ind]; ind++;
*zre += *yr0b*t_cos[*ind] + *yi0b*t_sin[*ind];
*zim += *yr0b*t_sin[*ind] - *yi0b*t_cos[*ind];
yr0b--; yi0b--;
yr0f++; yi0f++;
yr1f++; yi1f++;
zre++; zim--;
}
ind = &Ind[0];
for (j = 0; j < 3; j++) {
*ind = (short) (*ind+step);
if (*ind >= N_MAX)
{
*ind -= N_MAX;
}
ind++;
}
ind = &Ind[0];
*zre = *yr0f*t_cos[*ind] - *yi0f*t_sin[*ind]; ind++;
*zre += *yr1f*(t_cos[*ind] - t_sin[*ind]); ind++;
*zre += *yr0b*t_cos[*ind] + *yi0b*t_sin[*ind]; /* I have some doubts here, must check !!!!*/
/* Update the table step. */
step = (short)(step+3*(1<<(ORDER_MAX-order))); /* triple the step if we want to use va's cosine table*/
}
/* Compute the inverse FFT for all nine blocks. */
z0 = &Z[0];
z1 = &z0[m]; /* z1 = &Z[ m]; */
z2 = &z1[m]; /* z2 = &Z[2m]; */
ifft_rel(&z0[0], m, order);
ifft_rel(&z1[0], m, order);
ifft_rel(&z2[0], m, order);
/* Decimation and scaling, scale = 1/9m. */
scale =((float)1/3)/((float)m);
for (i = 0; i < n/3; i++) {
*X++ = (*z0++)*scale;
*X++ = (*z1++)*scale;
*X++ = (*z2++)*scale;
}
}