123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595 |
- // Copyright (c) 1996 Jens Jorgen Nielsen
- // Written by Jens Jorgen Nielsen
- // Optimizations by Janik Joire
- //
- // $History: $
- #include <math.h>
- #include <stdio.h>
- #include <stdlib.h>
- #ifndef PI
- #define PI 3.14159265358979323846F
- #endif
- /************************************************************************
- fft(long n,float xRe[],float xIm[],float yRe[],float yIm[])
- ------------------------------------------------------------------------
- NOTE : This is copyrighted material, Not public domain. See below.
- ------------------------------------------------------------------------
- Input/output:
- long n transformation length.
- float xRe[] real part of input sequence.
- float xIm[] imaginary part of input sequence.
- float yRe[] real part of output sequence.
- float yIm[] imaginary part of output sequence.
- ------------------------------------------------------------------------
- Function:
- The procedure performs a fast discrete Fourier transform (FFT) of
- a complex sequence, x, of an arbitrary length, n. The output, y,
- is also a complex sequence of length n.
-
- y[k] = sum(x[m]*exp(-i*2*pi*k*m/n), m=0..(n-1)), k=0,...,(n-1)
-
- The largest prime factor of n must be less than or equal to the
- constant maxPrimeFactor defined below.
- ------------------------------------------------------------------------
- Author:
- Jens Jorgen Nielsen For non-commercial use only.
- Bakkehusene 54 A $100 fee must be paid if used
- DK-2970 Horsholm commercially. Please contact.
- DENMARK
-
- e-mail : jnielsen@internet.dk All rights reserved. March 1996.
- ------------------------------------------------------------------------
- Implementation notes:
- The general idea is to factor the length of the DFT, n, into
- factors that are efficiently handled by the routines.
-
- A number of short DFT's are implemented with a minimum of
- arithmetical operations and using (almost) straight line code
- resulting in very fast execution when the factors of n belong
- to this set. Especially radix-10 is optimized.
- Prime factors, that are not in the set of short DFT's are handled
- with direct evaluation of the DFP expression.
- Please report any problems to the author.
- Suggestions and improvements are welcomed.
- ------------------------------------------------------------------------
- Benchmarks:
- The Microsoft Visual C++ compiler was used with the following
- compile options:
- /nologo /Gs /G2 /W4 /AH /Ox /D "NDEBUG" /D "_DOS" /FR
- and the FFTBENCH test executed on a 50MHz 486DX :
-
- Length Time [s] Accuracy [dB]
- 128 0.0054 -314.8
- 256 0.0116 -309.8
- 512 0.0251 -290.8
- 1024 0.0567 -313.6
- 2048 0.1203 -306.4
- 4096 0.2600 -291.8
- 8192 0.5800 -305.1
- 100 0.0040 -278.5
- 200 0.0099 -280.3
- 500 0.0256 -278.5
- 1000 0.0540 -278.5
- 2000 0.1294 -280.6
- 5000 0.3300 -278.4
- 10000 0.7133 -278.5
- ------------------------------------------------------------------------
- The following procedures are used :
- factorize : factor the transformation length.
- transTableSetup : setup table with sofar-, actual-, and remainRadix.
- permute : permutation allows in-place calculations.
- twiddleTransf : twiddle multiplications and DFT's for one stage.
- initTrig : initialise sine/cosine table.
- fft_4 : length 4 DFT, a la Nussbaumer.
- fft_5 : length 5 DFT, a la Nussbaumer.
- fft_10 : length 10 DFT using prime factor FFT.
- fft_odd : length n DFT, n odd.
- *************************************************************************/
- #define maxPrimeFactor 37
- #define maxPrimeFactorDiv2 (maxPrimeFactor+1)/2
- #define maxFactorCount 20
- static float c3_1 = -1.5000000000000E+00F; /* c3_1 = cos(2*pi/3)-1; */
- static float c3_2 = 8.6602540378444E-01F; /* c3_2 = sin(2*pi/3); */
-
- static float u5 = 1.2566370614359E+00F; /* u5 = 2*pi/5; */
- static float c5_1 = -1.2500000000000E+00F; /* c5_1 = (cos(u5)+cos(2*u5))/2-1;*/
- static float c5_2 = 5.5901699437495E-01F; /* c5_2 = (cos(u5)-cos(2*u5))/2; */
- static float c5_3 = -9.5105651629515E-01F; /* c5_3 = -sin(u5); */
- static float c5_4 = -1.5388417685876E+00F; /* c5_4 = -(sin(u5)+sin(2*u5)); */
- static float c5_5 = 3.6327126400268E-01F; /* c5_5 = (sin(u5)-sin(2*u5)); */
- static float c8 = 7.0710678118655E-01F; /* c8 = 1/sqrt(2); */
- static float pi;
- static long groupOffset,dataOffset,blockOffset,adr;
- static long groupNo,dataNo,blockNo,twNo;
- static float omega, tw_re,tw_im;
- static float twiddleRe[maxPrimeFactor], twiddleIm[maxPrimeFactor],
- trigRe[maxPrimeFactor], trigIm[maxPrimeFactor],
- zRe[maxPrimeFactor], zIm[maxPrimeFactor];
- static float vRe[maxPrimeFactorDiv2], vIm[maxPrimeFactorDiv2];
- static float wRe[maxPrimeFactorDiv2], wIm[maxPrimeFactorDiv2];
- void factorize(long n, long *nFact, long fact[])
- {
- long i,j,k;
- long nRadix;
- long radices[7];
- long factors[maxFactorCount];
- nRadix = 6;
- radices[1]= 2;
- radices[2]= 3;
- radices[3]= 4;
- radices[4]= 5;
- radices[5]= 8;
- radices[6]= 10;
- if (n==1)
- {
- j=1;
- factors[1]=1;
- }
- else j=0;
- i=nRadix;
- while ((n>1) && (i>0))
- {
- if ((n % radices[i]) == 0)
- {
- n=n / radices[i];
- j=j+1;
- factors[j]=radices[i];
- }
- else i=i-1;
- }
- if (factors[j] == 2) /*substitute factors 2*8 with 4*4 */
- {
- i = j-1;
- while ((i>0) && (factors[i] != 8)) i--;
- if (i>0)
- {
- factors[j] = 4;
- factors[i] = 4;
- }
- }
- if (n>1)
- {
- for (k=2; k<sqrt(n)+1; k++)
- while ((n % k) == 0)
- {
- n=n / k;
- j=j+1;
- factors[j]=k;
- }
- if (n>1)
- {
- j=j+1;
- factors[j]=n;
- }
- }
- for (i=1; i<=j; i++)
- {
- fact[i] = factors[j-i+1];
- }
- *nFact=j;
- } /* factorize */
- /****************************************************************************
- After N is factored the parameters that control the stages are generated.
- For each stage we have:
- sofar : the product of the radices so far.
- actual : the radix handled in this stage.
- remain : the product of the remaining radices.
- ****************************************************************************/
- void transTableSetup(long sofar[], long actual[], long remain[],
- long *nFact,
- long *nPoints)
- {
- long i;
- factorize(*nPoints, nFact, actual);
- if (actual[*nFact] > maxPrimeFactor)
- {
- printf("\nPrime factor of FFT length too large : %6d",actual[*nFact]);
- exit(1);
- }
- remain[0]=*nPoints;
- sofar[1]=1;
- remain[1]=*nPoints / actual[1];
- for (i=2; i<=*nFact; i++)
- {
- sofar[i]=sofar[i-1]*actual[i-1];
- remain[i]=remain[i-1] / actual[i];
- }
- } /* transTableSetup */
- /****************************************************************************
- The sequence y is the permuted input sequence x so that the following
- transformations can be performed in-place, and the final result is the
- normal order.
- ****************************************************************************/
- void permute(long nPoint, long nFact,
- long fact[], long remain[],
- float xRe[], float xIm[],
- float yRe[], float yIm[])
- {
- long i,j,k;
- long count[maxFactorCount];
- for (i=1; i<=nFact; i++) count[i]=0;
- k=0;
- for (i=0; i<=nPoint-2; i++)
- {
- yRe[i] = xRe[k];
- yIm[i] = xIm[k];
- j=1;
- k=k+remain[j];
- count[1] = count[1]+1;
- while (count[j] >= fact[j])
- {
- count[j]=0;
- k=k-remain[j-1]+remain[j+1];
- j=j+1;
- count[j]=count[j]+1;
- }
- }
- yRe[nPoint-1]=xRe[nPoint-1];
- yIm[nPoint-1]=xIm[nPoint-1];
- } /* permute */
- /****************************************************************************
- Twiddle factor multiplications and transformations are performed on a
- group of data. The number of multiplications with 1 are reduced by skipping
- the twiddle multiplication of the first stage and of the first group of the
- following stages.
- ***************************************************************************/
- void initTrig(long radix)
- {
- long i;
- float w,xre,xim;
- w=2*pi/radix;
- trigRe[0]=1; trigIm[0]=0;
- xre=(float)cos(w);
- xim=-(float)sin(w);
- trigRe[1]=xre; trigIm[1]=xim;
- for (i=2; i<radix; i++)
- {
- trigRe[i]=xre*trigRe[i-1] - xim*trigIm[i-1];
- trigIm[i]=xim*trigRe[i-1] + xre*trigIm[i-1];
- }
- } /* initTrig */
- void fft_4(float aRe[], float aIm[])
- {
- float t1_re,t1_im, t2_re,t2_im;
- float m2_re,m2_im, m3_re,m3_im;
- t1_re=aRe[0] + aRe[2]; t1_im=aIm[0] + aIm[2];
- t2_re=aRe[1] + aRe[3]; t2_im=aIm[1] + aIm[3];
- m2_re=aRe[0] - aRe[2]; m2_im=aIm[0] - aIm[2];
- m3_re=aIm[1] - aIm[3]; m3_im=aRe[3] - aRe[1];
- aRe[0]=t1_re + t2_re; aIm[0]=t1_im + t2_im;
- aRe[2]=t1_re - t2_re; aIm[2]=t1_im - t2_im;
- aRe[1]=m2_re + m3_re; aIm[1]=m2_im + m3_im;
- aRe[3]=m2_re - m3_re; aIm[3]=m2_im - m3_im;
- } /* fft_4 */
- void fft_5(float aRe[], float aIm[])
- {
- float t1_re,t1_im, t2_re,t2_im, t3_re,t3_im;
- float t4_re,t4_im, t5_re,t5_im;
- float m2_re,m2_im, m3_re,m3_im, m4_re,m4_im;
- float m1_re,m1_im, m5_re,m5_im;
- float s1_re,s1_im, s2_re,s2_im, s3_re,s3_im;
- float s4_re,s4_im, s5_re,s5_im;
- t1_re=aRe[1] + aRe[4]; t1_im=aIm[1] + aIm[4];
- t2_re=aRe[2] + aRe[3]; t2_im=aIm[2] + aIm[3];
- t3_re=aRe[1] - aRe[4]; t3_im=aIm[1] - aIm[4];
- t4_re=aRe[3] - aRe[2]; t4_im=aIm[3] - aIm[2];
- t5_re=t1_re + t2_re; t5_im=t1_im + t2_im;
- aRe[0]=aRe[0] + t5_re; aIm[0]=aIm[0] + t5_im;
- m1_re=c5_1*t5_re; m1_im=c5_1*t5_im;
- m2_re=c5_2*(t1_re - t2_re); m2_im=c5_2*(t1_im - t2_im);
- m3_re=-c5_3*(t3_im + t4_im); m3_im=c5_3*(t3_re + t4_re);
- m4_re=-c5_4*t4_im; m4_im=c5_4*t4_re;
- m5_re=-c5_5*t3_im; m5_im=c5_5*t3_re;
- s3_re=m3_re - m4_re; s3_im=m3_im - m4_im;
- s5_re=m3_re + m5_re; s5_im=m3_im + m5_im;
- s1_re=aRe[0] + m1_re; s1_im=aIm[0] + m1_im;
- s2_re=s1_re + m2_re; s2_im=s1_im + m2_im;
- s4_re=s1_re - m2_re; s4_im=s1_im - m2_im;
- aRe[1]=s2_re + s3_re; aIm[1]=s2_im + s3_im;
- aRe[2]=s4_re + s5_re; aIm[2]=s4_im + s5_im;
- aRe[3]=s4_re - s5_re; aIm[3]=s4_im - s5_im;
- aRe[4]=s2_re - s3_re; aIm[4]=s2_im - s3_im;
- } /* fft_5 */
- void fft_8()
- {
- float aRe[4], aIm[4], bRe[4], bIm[4], gem;
- aRe[0] = zRe[0]; bRe[0] = zRe[1];
- aRe[1] = zRe[2]; bRe[1] = zRe[3];
- aRe[2] = zRe[4]; bRe[2] = zRe[5];
- aRe[3] = zRe[6]; bRe[3] = zRe[7];
- aIm[0] = zIm[0]; bIm[0] = zIm[1];
- aIm[1] = zIm[2]; bIm[1] = zIm[3];
- aIm[2] = zIm[4]; bIm[2] = zIm[5];
- aIm[3] = zIm[6]; bIm[3] = zIm[7];
- fft_4(aRe, aIm); fft_4(bRe, bIm);
- gem = c8*(bRe[1] + bIm[1]);
- bIm[1] = c8*(bIm[1] - bRe[1]);
- bRe[1] = gem;
- gem = bIm[2];
- bIm[2] =-bRe[2];
- bRe[2] = gem;
- gem = c8*(bIm[3] - bRe[3]);
- bIm[3] =-c8*(bRe[3] + bIm[3]);
- bRe[3] = gem;
-
- zRe[0] = aRe[0] + bRe[0]; zRe[4] = aRe[0] - bRe[0];
- zRe[1] = aRe[1] + bRe[1]; zRe[5] = aRe[1] - bRe[1];
- zRe[2] = aRe[2] + bRe[2]; zRe[6] = aRe[2] - bRe[2];
- zRe[3] = aRe[3] + bRe[3]; zRe[7] = aRe[3] - bRe[3];
- zIm[0] = aIm[0] + bIm[0]; zIm[4] = aIm[0] - bIm[0];
- zIm[1] = aIm[1] + bIm[1]; zIm[5] = aIm[1] - bIm[1];
- zIm[2] = aIm[2] + bIm[2]; zIm[6] = aIm[2] - bIm[2];
- zIm[3] = aIm[3] + bIm[3]; zIm[7] = aIm[3] - bIm[3];
- } /* fft_8 */
- void fft_10()
- {
- float aRe[5], aIm[5], bRe[5], bIm[5];
- aRe[0] = zRe[0]; bRe[0] = zRe[5];
- aRe[1] = zRe[2]; bRe[1] = zRe[7];
- aRe[2] = zRe[4]; bRe[2] = zRe[9];
- aRe[3] = zRe[6]; bRe[3] = zRe[1];
- aRe[4] = zRe[8]; bRe[4] = zRe[3];
- aIm[0] = zIm[0]; bIm[0] = zIm[5];
- aIm[1] = zIm[2]; bIm[1] = zIm[7];
- aIm[2] = zIm[4]; bIm[2] = zIm[9];
- aIm[3] = zIm[6]; bIm[3] = zIm[1];
- aIm[4] = zIm[8]; bIm[4] = zIm[3];
- fft_5(aRe, aIm); fft_5(bRe, bIm);
- zRe[0] = aRe[0] + bRe[0]; zRe[5] = aRe[0] - bRe[0];
- zRe[6] = aRe[1] + bRe[1]; zRe[1] = aRe[1] - bRe[1];
- zRe[2] = aRe[2] + bRe[2]; zRe[7] = aRe[2] - bRe[2];
- zRe[8] = aRe[3] + bRe[3]; zRe[3] = aRe[3] - bRe[3];
- zRe[4] = aRe[4] + bRe[4]; zRe[9] = aRe[4] - bRe[4];
- zIm[0] = aIm[0] + bIm[0]; zIm[5] = aIm[0] - bIm[0];
- zIm[6] = aIm[1] + bIm[1]; zIm[1] = aIm[1] - bIm[1];
- zIm[2] = aIm[2] + bIm[2]; zIm[7] = aIm[2] - bIm[2];
- zIm[8] = aIm[3] + bIm[3]; zIm[3] = aIm[3] - bIm[3];
- zIm[4] = aIm[4] + bIm[4]; zIm[9] = aIm[4] - bIm[4];
- } /* fft_10 */
- void fft_odd(long radix)
- {
- float rere, reim, imre, imim;
- long i,j,k,n,max;
- n = radix;
- max = (n + 1)/2;
- for (j=1; j < max; j++)
- {
- vRe[j] = zRe[j] + zRe[n-j];
- vIm[j] = zIm[j] - zIm[n-j];
- wRe[j] = zRe[j] - zRe[n-j];
- wIm[j] = zIm[j] + zIm[n-j];
- }
- for (j=1; j < max; j++)
- {
- zRe[j]=zRe[0];
- zIm[j]=zIm[0];
- zRe[n-j]=zRe[0];
- zIm[n-j]=zIm[0];
- k=j;
- for (i=1; i < max; i++)
- {
- rere = trigRe[k] * vRe[i];
- imim = trigIm[k] * vIm[i];
- reim = trigRe[k] * wIm[i];
- imre = trigIm[k] * wRe[i];
-
- zRe[n-j] += rere + imim;
- zIm[n-j] += reim - imre;
- zRe[j] += rere - imim;
- zIm[j] += reim + imre;
- k = k + j;
- if (k >= n) k = k - n;
- }
- }
- for (j=1; j < max; j++)
- {
- zRe[0]=zRe[0] + vRe[j];
- zIm[0]=zIm[0] + wIm[j];
- }
- } /* fft_odd */
- void twiddleTransf(long sofarRadix, long radix, long remainRadix,
- float yRe[], float yIm[])
- { /* twiddleTransf */
- float cosw, sinw, gem;
- float t1_re,t1_im, t2_re,t2_im, t3_re,t3_im;
- float t4_re,t4_im, t5_re,t5_im;
- float m2_re,m2_im, m3_re,m3_im, m4_re,m4_im;
- float m1_re,m1_im, m5_re,m5_im;
- float s1_re,s1_im, s2_re,s2_im, s3_re,s3_im;
- float s4_re,s4_im, s5_re,s5_im;
- initTrig(radix);
- omega = 2*pi/(float)(sofarRadix*radix);
- cosw = (float)cos(omega);
- sinw = -(float)sin(omega);
- tw_re = 1.0;
- tw_im = 0;
- dataOffset=0;
- groupOffset=dataOffset;
- adr=groupOffset;
- for (dataNo=0; dataNo<sofarRadix; dataNo++)
- {
- if (sofarRadix>1)
- {
- twiddleRe[0] = 1.0;
- twiddleIm[0] = 0.0;
- twiddleRe[1] = tw_re;
- twiddleIm[1] = tw_im;
- for (twNo=2; twNo<radix; twNo++)
- {
- twiddleRe[twNo]=tw_re*twiddleRe[twNo-1]
- - tw_im*twiddleIm[twNo-1];
- twiddleIm[twNo]=tw_im*twiddleRe[twNo-1]
- + tw_re*twiddleIm[twNo-1];
- }
- gem = cosw*tw_re - sinw*tw_im;
- tw_im = sinw*tw_re + cosw*tw_im;
- tw_re = gem;
- }
- for (groupNo=0; groupNo<remainRadix; groupNo++)
- {
- if ((sofarRadix>1) && (dataNo > 0))
- {
- zRe[0]=yRe[adr];
- zIm[0]=yIm[adr];
- blockNo=1;
- do {
- adr = adr + sofarRadix;
- zRe[blockNo]= twiddleRe[blockNo] * yRe[adr]
- - twiddleIm[blockNo] * yIm[adr];
- zIm[blockNo]= twiddleRe[blockNo] * yIm[adr]
- + twiddleIm[blockNo] * yRe[adr];
-
- blockNo++;
- } while (blockNo < radix);
- }
- else
- for (blockNo=0; blockNo<radix; blockNo++)
- {
- zRe[blockNo]=yRe[adr];
- zIm[blockNo]=yIm[adr];
- adr=adr+sofarRadix;
- }
- switch(radix) {
- case 2 : gem=zRe[0] + zRe[1];
- zRe[1]=zRe[0] - zRe[1]; zRe[0]=gem;
- gem=zIm[0] + zIm[1];
- zIm[1]=zIm[0] - zIm[1]; zIm[0]=gem;
- break;
- case 3 : t1_re=zRe[1] + zRe[2]; t1_im=zIm[1] + zIm[2];
- zRe[0]=zRe[0] + t1_re; zIm[0]=zIm[0] + t1_im;
- m1_re=c3_1*t1_re; m1_im=c3_1*t1_im;
- m2_re=c3_2*(zIm[1] - zIm[2]);
- m2_im=c3_2*(zRe[2] - zRe[1]);
- s1_re=zRe[0] + m1_re; s1_im=zIm[0] + m1_im;
- zRe[1]=s1_re + m2_re; zIm[1]=s1_im + m2_im;
- zRe[2]=s1_re - m2_re; zIm[2]=s1_im - m2_im;
- break;
- case 4 : t1_re=zRe[0] + zRe[2]; t1_im=zIm[0] + zIm[2];
- t2_re=zRe[1] + zRe[3]; t2_im=zIm[1] + zIm[3];
- m2_re=zRe[0] - zRe[2]; m2_im=zIm[0] - zIm[2];
- m3_re=zIm[1] - zIm[3]; m3_im=zRe[3] - zRe[1];
- zRe[0]=t1_re + t2_re; zIm[0]=t1_im + t2_im;
- zRe[2]=t1_re - t2_re; zIm[2]=t1_im - t2_im;
- zRe[1]=m2_re + m3_re; zIm[1]=m2_im + m3_im;
- zRe[3]=m2_re - m3_re; zIm[3]=m2_im - m3_im;
- break;
- case 5 : t1_re=zRe[1] + zRe[4]; t1_im=zIm[1] + zIm[4];
- t2_re=zRe[2] + zRe[3]; t2_im=zIm[2] + zIm[3];
- t3_re=zRe[1] - zRe[4]; t3_im=zIm[1] - zIm[4];
- t4_re=zRe[3] - zRe[2]; t4_im=zIm[3] - zIm[2];
- t5_re=t1_re + t2_re; t5_im=t1_im + t2_im;
- zRe[0]=zRe[0] + t5_re; zIm[0]=zIm[0] + t5_im;
- m1_re=c5_1*t5_re; m1_im=c5_1*t5_im;
- m2_re=c5_2*(t1_re - t2_re);
- m2_im=c5_2*(t1_im - t2_im);
- m3_re=-c5_3*(t3_im + t4_im);
- m3_im=c5_3*(t3_re + t4_re);
- m4_re=-c5_4*t4_im; m4_im=c5_4*t4_re;
- m5_re=-c5_5*t3_im; m5_im=c5_5*t3_re;
- s3_re=m3_re - m4_re; s3_im=m3_im - m4_im;
- s5_re=m3_re + m5_re; s5_im=m3_im + m5_im;
- s1_re=zRe[0] + m1_re; s1_im=zIm[0] + m1_im;
- s2_re=s1_re + m2_re; s2_im=s1_im + m2_im;
- s4_re=s1_re - m2_re; s4_im=s1_im - m2_im;
- zRe[1]=s2_re + s3_re; zIm[1]=s2_im + s3_im;
- zRe[2]=s4_re + s5_re; zIm[2]=s4_im + s5_im;
- zRe[3]=s4_re - s5_re; zIm[3]=s4_im - s5_im;
- zRe[4]=s2_re - s3_re; zIm[4]=s2_im - s3_im;
- break;
- case 8 : fft_8(); break;
- case 10 : fft_10(); break;
- default : fft_odd(radix); break;
- }
- adr=groupOffset;
- for (blockNo=0; blockNo<radix; blockNo++)
- {
- yRe[adr]=zRe[blockNo]; yIm[adr]=zIm[blockNo];
- adr=adr+sofarRadix;
- }
- groupOffset=groupOffset+sofarRadix*radix;
- adr=groupOffset;
- }
- dataOffset=dataOffset+1;
- groupOffset=dataOffset;
- adr=groupOffset;
- }
- } /* twiddleTransf */
- void fft(long n,float *xRe,float *xIm,float *yRe,float *yIm)
- {
- long sofarRadix[maxFactorCount],
- actualRadix[maxFactorCount],
- remainRadix[maxFactorCount];
- long nFactor;
- long count;
- pi = PI;
- transTableSetup(sofarRadix, actualRadix, remainRadix, &nFactor, &n);
- permute(n, nFactor, actualRadix, remainRadix, xRe, xIm, yRe, yIm);
- for (count=1; count<=nFactor; count++)
- twiddleTransf(sofarRadix[count], actualRadix[count], remainRadix[count],
- yRe, yIm);
- } /* fft */
|