[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Code for Fast Fourier Transform



PureBytes Links

Trading Reference Links


You won't be able to do it in EL with TS 4 due to the 64k program size
limit.  The only way to do it with TS 4 is with a .dll.  Suggestion:
Compile the C code provided into a TS4-style (16-bit) .dll ...


The Omega Man



----- Original Message -----
From: Jasper Riedel <jrbox@xxxxxxx>
To: omega-list <omega-list@xxxxxxxxxx>
Sent: Thursday, June 17, 1999 5:26 PM
Subject: Re: Code for Fast Fourier Transform


> >Can you send me code for the FFT?  Thanks.
>
> Here you go ... but its in C ...
> What i'm looking for is a TESTED code for EL to use in TS4.
>
> Regards - jr
>
>
> ---- cut here ---
>
> // FLT_TYPE is either typedef'd as double or float
>
> BOOL DoFFT(int nPow2Data, FLT_TYPE* pInData, FLT_TYPE* fftAmpdata,
FLT_TYPE*
> fftPhidata, FLT_TYPE fEffValue, FLT_TYPE fLowBorder)
> {
> FLT_TYPE totalcos;
> FLT_TYPE totalsin;
> FLT_TYPE totalamplit;
>
> FLT_TYPE* cosdata = (FLT_TYPE*) malloc(nPow2Data * sizeof(FLT_TYPE));
> FLT_TYPE* sindata = (FLT_TYPE*) malloc(nPow2Data * sizeof(FLT_TYPE));
> if (!cosdata || !sindata) {
>   free(cosdata);
>   free(sindata);
>   return FALSE;
> }
> memset(sindata, 0, nPow2Data * sizeof(FLT_TYPE));
>
> int n, nPow2Exp;
> for (n=1, nPow2Exp=0; n < nPow2Data; ++nPow2Exp, n *= 2);
>
> for (n=0; n<nPow2Data; ++n)
>   cosdata[n] = pInData[n];
>
> RawFft(cosdata, sindata, nPow2Exp, 1);
>
> totalcos = totalsin = 0;
>
> for (n=0; n<nPow2Data / 2; ++n) {
>
>   fftAmpdata[n] = (FLT_TYPE) sqrt(cosdata[n] * cosdata[n] + sindata[n] *
> sindata[n]);
>
>   totalcos += cosdata[n];
>   totalsin += sindata[n];
> }
>
> totalamplit = (FLT_TYPE) sqrt(totalcos * totalcos + totalsin * totalsin);
>
> for (n=0; n<nPow2Data / 2; ++n) {
>
>   fftAmpdata[n] /= totalamplit; // scales from 0 to 1
>
>   if (fEffValue > 0) {
>    fftAmpdata[n] *= fEffValue; // scales to effVal
>   }
>
>   if (fLowBorder > 0) {
>    fftAmpdata[n] = max(fftAmpdata[n], fLowBorder);
>   }
> }
>
>
> static FLT_TYPE lastValid = 0;
>
> if (fftPhidata != NULL) {
>
>   lastValid = 0;
>
>   for (n=0; n<nPow2Data / 2; ++n) {
>
>    if (sindata[n] != 0 && cosdata[n] != 0) {
>
> //    fftPhidata[n] = (FLT_TYPE) (atan(cosdata[n] / sindata[n]) * 180. /
> PI);
>     fftPhidata[n] = (FLT_TYPE) ((atan2(cosdata[n], sindata[n])) * 90. /
PI);
>
>     lastValid = fftPhidata[n];
>    }
>    else fftPhidata[n] = lastValid;
>
>   }
> }
>
> free(cosdata);
> free(sindata);
>
> return TRUE;
> }
>
>
> /* PCFFT.C  -- by J.G.G. Dobbe -- Performs an FFT on two arrays (Re, Im)
of
>                type float (can be changed). This unit is written in C and
>                doesn't call assembler routines.
> */
>
> /* --------------------- Include directive ------------------------ */
> // #include "pcfft.h" - gone
>
> /* --------------------- Local variables -------------------------- */
>
> static FLT_TYPE CosArray[28] =
> { /* cos{-2pi/N} for N = 2, 4, 8, ... 16384 */
> -1.00000000000000,  0.00000000000000,  0.70710678118655,
>   0.92387953251129,  0.98078528040323,  0.99518472667220,
>   0.99879545620517,  0.99969881869620,  0.99992470183914,
>   0.99998117528260,  0.99999529380958,  0.99999882345170,
>   0.99999970586288,  0.99999992646572,
>   /* cos{2pi/N} for N = 2, 4, 8, ... 16384 */
> -1.00000000000000,  0.00000000000000,  0.70710678118655,
>   0.92387953251129,  0.98078528040323,  0.99518472667220,
>   0.99879545620517,  0.99969881869620,  0.99992470183914,
>   0.99998117528260,  0.99999529380958,  0.99999882345170,
>   0.99999970586288,  0.99999992646572
> };
> static FLT_TYPE SinArray[28] =
> { /* sin{-2pi/N} for N = 2, 4, 8, ... 16384 */
>   0.00000000000000, -1.00000000000000, -0.70710678118655,
> -0.38268343236509, -0.19509032201613, -0.09801714032956,
> -0.04906767432742, -0.02454122852291, -0.01227153828572,
> -0.00613588464915, -0.00306795676297, -0.00153398018628,
> -0.00076699031874, -0.00038349518757,
>   /* sin{2pi/N} for N = 2, 4, 8, ... 16384 */
>   0.00000000000000,  1.00000000000000,  0.70710678118655,
>   0.38268343236509,  0.19509032201613,  0.09801714032956,
>   0.04906767432742,  0.02454122852291,  0.01227153828572,
>   0.00613588464915,  0.00306795676297,  0.00153398018628,
>   0.00076699031874,  0.00038349518757
> };
>
> /* --------------------- Function implementations ----------------- */
> /* --------------------- ShuffleIndex ----------------------------- */
> static inline unsigned int ShuffleIndex(unsigned int i, int WordLength)
>
> /* Function     : Finds the shuffle index of array elements. The array
> length
>                   must be a power of two; The power is stored in
> "WordLength".
>    Return value : With "i" the source array index, "ShuffleIndex"
>                   returns the destination index for shuffling.
>    Comment      : -
> */
> {
>   unsigned int  NewIndex;
>   unsigned char BitNr;
>   NewIndex = 0;
>   for (BitNr = 0; BitNr <= WordLength - 1; BitNr++)
>   {
>     NewIndex = NewIndex << 1;
>     if ((i & 1) != 0) NewIndex = NewIndex + 1;
>     i = i >> 1;
>   }
>   return NewIndex;
> }
> /* --------------------- Shuffle2Arr ------------------------------ */
> static void Shuffle2Arr(FLT_TYPE *a, FLT_TYPE *b, int bitlength, unsigned
> int N)
> /* Function     : Shuffles both arrays "a" and "b". This function is
called
>                  before performing the actual FFT so the array elements
>                  are in the right order after FFT.
>   Return value : -
>   Comment      : -
> */
> {
>   unsigned int  IndexOld, IndexNew;
>   FLT_TYPE      temp;
>
> // //  int bitlengthtemp = bitlength;                  /* Save for later
use
> */
> //  N = 1;                                       /* Find array-length */
> //  do
> //  {
> //    N = N * 2;
> //    bitlengthtemp = bitlengthtemp - 1;
> //  } while (bitlengthtemp > 0) ;
>
>                                             /* Shuffle all elements */
>   for (IndexOld = 0; IndexOld <= N - 1; IndexOld++)
>   {                              /* Find index to exchange elements */
>     IndexNew = ShuffleIndex(IndexOld, bitlength);
>     if (IndexNew > IndexOld)
>     {                                         /* Exchange elements: */
>       temp = a[IndexOld];                             /* Of array a */
>       a[IndexOld] = a[IndexNew];
>       a[IndexNew] = temp;
>       temp = b[IndexOld];                             /* Of array a */
>       b[IndexOld] = b[IndexNew];
>       b[IndexNew] = temp;
>     }
>   }
> }
> /* --------------------- Fft -------------------------------------- */
> void RawFft(FLT_TYPE *Re, FLT_TYPE *Im, int Pwr, int Dir)
>
> /* Function     : Actual FFT algorithm. "Re" and "Im" point to start of
real
>                  and imaginary arrays of numbers, "Pwr" holds the array
> sizes
>                  as a power of 2 while "Dir" indicates whether an FFT
> (Dir>=1)
>                  or an inverse FFT must be performed (Dir<=0).
>   Return value : The transformed information is returned by "Re"
>                  and "Im" (real and imaginary part respectively).
>   Comment      : -
> */
> {
>   int pwrhelp;
>   int N;
>   int Section;
>   int AngleCounter;
>   int FlyDistance;
>   int FlyCount;
>   int index1;
>   int index2;
>   FLT_TYPE tempr, tempi;
> // FLT_TYPE Re1, Re2, Im1, Im2;
> // FLT_TYPE sqrtn;
>   FLT_TYPE c, s;
>   FLT_TYPE scale;
>   FLT_TYPE temp;
>   FLT_TYPE Qr, Qi;
>
>
>   pwrhelp = Pwr;    /* Determine size of arrs */
>   N = 1;
>   do
>   {
>     N = N * 2;
>     pwrhelp--;
>   } while (pwrhelp > 0) ;
>
>   Shuffle2Arr(Re, Im, Pwr, N); /* Shuffle before (i)FFT */
>
>   if (Dir >= 1) AngleCounter = 0;  /* FFT */
>   else          AngleCounter = 14;  /* Inverse FFT */
>   Section = 1;
>   while (Section < N)
>   {
>     FlyDistance = 2 * Section;
>     c = CosArray[AngleCounter];
>     s = SinArray[AngleCounter];
>     Qr = 1; Qi = 0;
>     for (FlyCount = 0; FlyCount <= Section - 1; FlyCount++)
>     {
>       index1 = FlyCount;
>       do
>       {
>         index2 = index1 + Section;
>           /* Perform 2-Point DFT  */
>
>         tempr = (FLT_TYPE) 1.0 * Qr * Re[index2] - (FLT_TYPE) 1.0 * Qi *
> Im[index2];
>         tempi = (FLT_TYPE) 1.0 * Qr * Im[index2] + (FLT_TYPE) 1.0 * Qi *
> Re[index2];
>
>         Re[index2] = Re[index1] - tempr;  /* For Re-part */
>         Re[index1] = Re[index1] + tempr;
>         Im[index2] = Im[index1] - tempi;  /* For Im-part */
>         Im[index1] = Im[index1] + tempi;
>
>         index1 = index1 + FlyDistance;
>       } while (index1 <= (N - 1));
>
>       /*                 k                                  */
>       /*  Calculate new Q = cos(ak) + j*sin(ak) = Qr + j*Qi */
>       /*          -2*pi                                     */
>       /*  with: a = -----                                   */
>       /*            N                                       */
>       temp = Qr;
>       Qr   = Qr*c - Qi*s;
>       Qi   = Qi*c + temp*s;
>     }
>     Section = Section * 2;
>     AngleCounter = AngleCounter + 1;
>   }
>   if (Dir <= 0)        /* Normalize for */
>   {           /* inverse FFT only */
>     scale = (FLT_TYPE) 1.0/N;
>     for (index1 = 0; index1 <= N - 1; index1++)
>     {
>       Re[index1] = scale * Re[index1];
>       Im[index1] = scale * Im[index1];
>     }
>   }
> }
>
>
>
>
>
>
>
>
>