steve
/
winamp
mirror of https://github.com/WinampDesktop/winamp.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
winamp/Src/vlb/fft.cpp

318 lines
6.8 KiB
C++
Raw Blame History

/* $Header: /cvs/root/winamp/vlb/fft.cpp,v 1.1 2009/04/28 20:21:09 audiodsp Exp $ */
/***************************************************************************\
*
* Copyright 2000-2002 Dolby Laboratories, Inc. All Rights
* Reserved. Do not copy. Do not distribute.
* Confidential information.
*
* filename: fft.cpp
* project : MPEG-2 AAC Decoder
* contents/description: Standard FFT routine
*
\***************************************************************************/
/* Description: This function calculates a 1024,256 or 128 point fft via
decimation in frequency algorithm described in "Digital Signal
Processing" by Oppenheim and Schafer, refer to pages 304 (flow
chart) and 330-332 (Fotran program in problem 5). Modified to
get only desired block sizes. */
#include"fft.h"
#ifndef PI
#define PI 3.14159265358979
#endif
#define HOTSWAP(a,b) fTemp=a; a=b; b=fTemp;
void fftl(float *pfReal, float *pfImag,int iSize)
{
static int iInit=0;
static int iPower2=0;
static int iTransformSize=0;
static float *pfTwiddleReal=0;
static float *pfTwiddleImag=0;
if(iSize!=iTransformSize){
iInit=0;
}
if(!iInit){
int iCounter;
if((iSize-1)&iSize){
//printf("FFT Length Must be Radix 2!\n");
return;
}
free(pfTwiddleReal);
free(pfTwiddleImag);
iTransformSize=iSize;
pfTwiddleReal=(float*)malloc(iSize/2*sizeof(float));
pfTwiddleImag=(float*)malloc(iSize/2*sizeof(float));
for(iCounter=0;iCounter<iSize/2;iCounter++){
pfTwiddleReal[iCounter]=(float)(cos(2.0*PI/(double)iSize*(double)iCounter));
pfTwiddleImag[iCounter]=(float)(-1.0*sin(2.0*PI/(double)iSize*(double)iCounter));
}
iPower2=0;
for(iCounter=1;iCounter<iSize;iCounter<<=1) iPower2++;
iInit++;
}
/* Perform Fast Fourier Transform: */
/* With Help From Figure 6.11 pp 464 */
/* John G. Proakis,"Digital Signal Processing," 3rd Edt. */
/* 1996 Prentice Hall */
int iSkip;
int iFFTCounter;
int iCounter1;
int iCounter2;
int iOffset;
int iHSize;
iHSize=iSize>>1;
iFFTCounter=1;
for(iSkip=iHSize;iSkip>0;iSkip>>=1){
iOffset=0;
for(iCounter1=0;iCounter1<iFFTCounter;iCounter1++){
int iOffset1;
int iOffset2;
int iTwiddleOffset;
iOffset1=iOffset;
iOffset2=iOffset1+iSkip;
iTwiddleOffset=0;
for(iCounter2=0;iCounter2<iSkip;iCounter2++){
float fTemp;
fTemp=pfReal[iOffset1];
pfReal[iOffset1]+=pfReal[iOffset2];
pfReal[iOffset2]=fTemp-pfReal[iOffset2];
fTemp=pfImag[iOffset1];
pfImag[iOffset1]+=pfImag[iOffset2];
pfImag[iOffset2]=fTemp-pfImag[iOffset2];
if(iCounter2){
fTemp=pfReal[iOffset2];
pfReal[iOffset2]=pfReal[iOffset2]*pfTwiddleReal[iTwiddleOffset]-
pfImag[iOffset2]*pfTwiddleImag[iTwiddleOffset];
pfImag[iOffset2]=fTemp*pfTwiddleImag[iTwiddleOffset]+
pfImag[iOffset2]*pfTwiddleReal[iTwiddleOffset];
}
iOffset1++;
iOffset2++;
iTwiddleOffset+=iFFTCounter;
}
iOffset+=(iSkip<<1);
}
iFFTCounter<<=1;
}
/*Reorder Data: */
iCounter1=0;
for(iCounter2=0;iCounter2<iSize-1;iCounter2++){
int iCounter3;
if(iCounter2<iCounter1){
float fTemp;
HOTSWAP(pfReal[iCounter2],pfReal[iCounter1]);
HOTSWAP(pfImag[iCounter2],pfImag[iCounter1]);
}
iCounter3=iHSize;
while(iCounter3<=iCounter1){
iCounter1-=iCounter3;
iCounter3=iCounter3>>1;
}
iCounter1=iCounter1+iCounter3;
}
}
/*
void fftl(float* xreal,float* ximag,int N)
{
int m,mm1;
static int init=0;
int nv2,nm1,mp;
static float wreal[13][13],wimag[13][13];
int i,j,k,l;
int ip,le,le1;
float treal,timag,ureal,uimag;
if (init==0){
memset((char*) wreal,0,sizeof(wreal));
memset((char*)wimag,0,sizeof(wimag));
m=13;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[0][l]=(float)(cosl(pi/le1));
wimag[0][l]=-1.0f*(float)(sinl(pi/le1));
}
m=12;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[1][l]=(float)(cosl(pi/le1));
wimag[1][l]=-1.0f*(float)(sinl(pi/le1));
}
m=11;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[2][l]=(float)(cosl(pi/le1));
wimag[2][l]=-1.0f*(float)(sinl(pi/le1));
}
m=10;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[3][l]=(float)(cosl(pi/le1));
wimag[3][l]=-1.0f*(float)(sinl(pi/le1));
}
m=9;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[4][l]=(float)(cosl(pi/le1));
wimag[4][l]=-1.0f*(float)(sinl(pi/le1));
}
m=8;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[5][l]=(float)(cosl(pi/le1));
wimag[5][l]=-1.0f*(float)(sinl(pi/le1));
}
m=7;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[6][l]=(float)(cosl(pi/le1));
wimag[6][l]=-1.0f*(float)(sinl(pi/le1));
}
m=6;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[7][l]=(float)(cosl(pi/le1));
wimag[7][l]=-1.0f*(float)(sinl(pi/le1));
}
m=5;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[8][l]=(float)(cosl(pi/le1));
wimag[8][l]=-1.0f*(float)(sinl(pi/le1));
}
m=4;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[9][l]=(float)(cosl(pi/le1));
wimag[9][l]=-1.0f*(float)(sinl(pi/le1));
}
m=3;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[10][l]=(float)(cosl(pi/le1));
wimag[10][l]=-1.0f*(float)(sinl(pi/le1));
}
m=2;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[11][l]=(float)(cosl(pi/le1));
wimag[11][l]=-1.0f*(float)(sinl(pi/le1));
}
m=1;
for(l=0;l<m;l++){
le=1<<(m-l);
le1=le>>1;
wreal[12][l]=(float)(cosl(pi/le1));
wimag[12][l]=-1.0f*(float)(sinl(pi/le1));
}
init++;
}
switch(N){
case 8192: m=13;mp=0;break;
case 4096: m=12;mp=1;break;
case 2048: m=11;mp=2;break;
case 1024: m=10;mp=3;break;
case 512: m=9;mp=4;break;
case 256:m=8;mp=5;break;
case 128:m=7;mp=6;break;
case 64:m=6;mp=7;break;
case 32:m=5;mp=8;break;
case 16:m=4;mp=9;break;
case 8:m=3;mp=10;break;
case 4:m=2;mp=11;break;
case 2:m=1;mp=12;break;
default: printf("FFT ERROR! %d\n",N);return;
}
mm1=m-1;
nv2=N>>1;
nm1=N-1;
for(l=0;l<mm1;l++){
le=1<<(m-l);
le1=le>>1;
ureal=1;
uimag=0;
for(j=0;j<le1;j++){
for(i=j;i<N;i+=le){
ip=i+le1;
treal=xreal[i]+xreal[ip];
timag=ximag[i]+ximag[ip];
xreal[ip]=xreal[i]-xreal[ip];
ximag[ip]=ximag[i]-ximag[ip];
xreal[i]=treal;
ximag[i]=timag;
treal=xreal[ip];
xreal[ip]=xreal[ip]*ureal-ximag[ip]*uimag;
ximag[ip]=ximag[ip]*ureal+treal*uimag;
}
treal=ureal;
ureal=ureal*wreal[mp][l]-uimag*wimag[mp][l];
uimag=uimag*wreal[mp][l]+treal*wimag[mp][l];
}
}
for(i=0;i<N;i+=2){
ip=i+1;
treal=xreal[i]+xreal[ip];
timag=ximag[i]+ximag[ip];
xreal[ip]=xreal[i]-xreal[ip];
ximag[ip]=ximag[i]-ximag[ip];
xreal[i]=treal;
ximag[i]=timag;
}
j=0;
for(i=0;i<nm1;i++){
if(i<j){
treal=xreal[j];
timag=ximag[j];
xreal[j]=xreal[i];
ximag[j]=ximag[i];
xreal[i]=treal;
ximag[i]=timag;
}
k=nv2;
while(k<=j){
j=j-k;
k=k>>1;
}
j=j+k;
}
}
*/
Reference in New Issue View Git Blame Copy Permalink