|
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];
> }
> }
> }
>
>
>
>
>
>
>
>
>
|