www.pudn.com > AVS_M_ver10.rar > fft.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/fft.h"
#define NUM 17 //number of imput data
#define FFT_ORDER 32 //fft transform order
//#define INFFT_ORDER 8 //invfft transform order
/************************************************************************/
/* premyfft
* prepare coefficient for the fft transform
* argument:
* n(i): point number of fft
* mode(i): 1 represents reverse fft transform
* 0 represents fft transform
* nexp(o): the power of 2 as for the number of point
* w(o): the coefficient of fft transform
*/
/************************************************************************/
void premyfft(int n,int mode,int *nexp,Complex *w)
{
int i,k,nt,nexp1;
float s;
Complex c1,c2;
nexp1=1;nt=1;
do{
nt=1;
for(i=1;i<=nexp1;i++)
nt=2*nt;
if(nt>=n)
break;
nexp1=nexp1+1;
}while(1);
if(nt==n)
{
s=(float)(8*atan(1)/(float)(nt));
c1=set(cos(s),-sin(s));
if(mode!=0) c1=conjg(c1);
c2=set(1.,0);
for(k=1;k<=nt;k++)
{
w[k]=c2;
c2=cmulc(c2,c1);
}
}
else
{
nexp1=-1;
}
*nexp=nexp1;
}
/************************************************************************/
/* myfft
* fft transform
* argument:
* n(i): the number of fft
* t(i): 1.0?
* mode(i): 1 represents reverse fft transform
* 0 represents fft transform
* nexp(i): the power of 2 as for the number of point
* w(i): the coefficient of fft transform
* x(i/o): original signal/fft result(the index of signal is
from 1 to n)
*/
/************************************************************************/
void myfft(int n,int mode,float t,int nexp,Complex *w,Complex *x)
{
Complex c1,c2;
int k,mm,ll,j,jj,kk,i,nn,nv2,nm1;
float s;
mm=1;
ll=n;
for(k=1;k<=nexp;k++)
{
nn=ll/2;
jj=mm+1;
for(i=1;i<=n;i=i+ll)
{
kk=i+nn;
c1=caddc(x[i],x[kk]);
x[kk]=csubc(x[i],x[kk]);
x[i]=c1;
}
if(nn==1)
continue;
else
{
for(j=2;j<=nn;j++)
{
c2=w[jj];
for(i=j;i<=n;i=i+ll)
{
kk=i+nn;
c1=caddc(x[i],x[kk]);
x[kk]=cmulc(csubc(x[i],x[kk]),c2);
x[i]=c1;
}
jj=jj+mm;
}
ll=nn;
mm=mm*2;
}
}
nv2=n/2;
nm1=n-1;
j=1;
for(i=1;i=j);
else
{
c1=x[j];
x[j]=x[i];
x[i]=c1;
}
k=nv2;
do
{
if(k>=j)
break;
else
{
j=j-k;
k=k/2;
}
}while(1);
j=j+k;
}
if(mode==0)
s=t;
else
s=(float)(1./(t*(float)n));
for(i=1;i<=n;i++)
x[i]=cmulf(x[i],s);
}
Complex set(double r,double i)
{
Complex temp;
temp.real=r;
temp.imag=i;
return temp;
}
Complex conjg(Complex x)
{
Complex temp;
temp.real=x.real;
temp.imag=-x.imag;
return temp;
}
Complex cmulc(Complex c1,Complex c2)
{
Complex temp;
temp.real=c1.real*c2.real -c1.imag *c2.imag;
temp.imag =c1.imag*c2.real +c1.real *c2.imag ;
return temp;
}
Complex cmulf(Complex c1,float c2)
{
Complex temp;
temp.real =c1.real *c2;
temp.imag =c1.imag *c2;
return temp;
}
Complex caddc(Complex c1,Complex c2)
{
Complex temp;
temp.real =c1.real +c2.real ;
temp.imag =c1.imag +c2.imag ;
return temp;
}
Complex csubc(Complex c1,Complex c2)
{
Complex temp;
temp.real =c1.real -c2.real ;
temp.imag =c1.imag -c2.imag ;
return temp;
}
Complex ctoc(Complex c1)
{
Complex temp;
temp.real =c1.real;
temp.imag =c1.imag;
return temp;
}
double cabsv(Complex c)
{
return sqrt( pow(c.real,2)+pow(c.imag,2.0) );
}
void output(Complex *x,int n)
{
int i;
for(i=0;i0)
{
*gamma1=0.96f;
*gamma2=0.54f;
*gamma3=0.4f;
/* *gamma1=0.95f;
*gamma2=0.58f;
*gamma3=0.5f;*/
}
else if( temp>-10.0 && temp<=0.0 )
{
*gamma1=0.95f;//0.95f;
*gamma2=0.74f;
*gamma3=0.18f;//0.1f;
}
else if(temp>-20.0 && temp<=-10.0)
{
/* *gamma1=0.95f;
*gamma2=0.5f;
*gamma3=0.1f;*/
*gamma1=0.92f;
*gamma2=0.5f;
*gamma3=0.13f;
}
else if(temp>-25.0 && temp<=-20.0)
{
*gamma1=0.91f;
*gamma2=0.45f;
*gamma3=0.2f;
}
else if(temp<=-25.0)
{
*gamma1=0.90f;
*gamma2=0.28f;
*gamma3=0.16f;
}
}