|
PureBytes Links
Trading Reference Links
|
>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];
}
}
}
|