Summary: Part of a larger work on time series analysis of the financial markets, especially the NYSE and the FRED files produced by the Federal Reserve Bank of St Louis, a time series analysis is a important tool in finding cyclical and redundant information in a data stream. The DJIA data stream has a history of 80+ years, so there is enough data ( 20,000+ data items per stream ) to get a 'good' DFT analysis. I do not know of any other publically available data streams that have the complexity of the DJIA. I had taken a stability approach, using a Z transform, and discovered that the Z transform signaled a stability problem the day after Ronald Reagan's 500 point 'correction'. Otherwise, the Z transform provided advanced warning of most all the major market meltdowns. Since the Z transform only works for a DFT when the data lies on the unit circle, this seems to have been a odd bedfellow in the quest for a market analysis tool. A better approach , in my opinion, would be a analysis of the phase angle output of the DFT. There seems to be a apparent correlation between phase angle changes and market changes (up to 21 days in advance) in the DJIA data stream. A listing containing the DFT's ( 32 and 64 bit versions for MS C++/x86 for VS 2008 ) with header file listing is included. These DFT's are not optimized for speed and/or efficiency. however, they do work !!!
#include "math.h"
void DFT64( double * , double * , double * , double * , int )
void DFT32(float* , float* , float* , float* , int );
void DFT64( double * dxx, double * cxx, double * th, double * fr, int cnt)
{
double dxe, cxe, theta, freq, rad, pi2, dne, dnx, cv, cvc;
int i, j, n;
n = cnt; // number of data items
rad = 360.0; // degrees in one cycle
__asm
{
;
finit; // init FPU
;
fldpi;
fldpi;
fadd;
fstp pi2; // pi2 = pi + pi : 2 * pi
;
fld rad;
fdiv pi2;
fstp rad; // rad = 360.0 / pi2 : rad = 57.3 deg/rad
;
fld pi2;
fchs; // change sign pi2 = -1.0 * pi2
fidiv n;
fstp pi2; // pi2 = (-1.0 * 2.0 * pi) / n
;
mov i,0; // i = 0 : outer loop counter
;
Loop1: ; // Top of outer Loop1
;
fild i;
fmul pi2;
fstp dnx; // dnx = pi2 * i
;
}
dxe = *(dxx + i); // dxe = dxx[i]
__asm
{
;
fldz;
fstp cv; // cv = 0 : power series accumulator
;
mov j,0; // j = 0 : inner loop counter
;
Loop2: ; // Top of inner Loop2
;
fild j;
fmul dnx;
fstp dne; // dne = dnx * j : dne = ( -1.0 * 2.0 * pi * i * j ) / n
;
}
cxe = exp(dne); // cxe = e^dne : c++ call ( math.h )
__asm
{
;
fld cxe;
fmul dxe;
fadd cv;
fstp cv; // cv = cv + ( cxe * dxe ) : calc power series
;
inc j; // j = j + 1 : increment inner loop counter
mov eax,j; // eax = j : cpu register eax = j
cmp eax,n; // eax ? n : n compared to cpu register eax
jl Loop2; // j < n ? if so jump to Loop2 label
;
fld cv;
fcos;
fmul cv;
fstp cvc; // cvc = (cos( cv ) * cv) : calc adjacent side
;
fld cvc;
fdiv cv;
fstp cvc; // cvc = cvc / cv : cvc = adjacent / hypoteneus
}
theta = acos(cvc); // theta = Acos( cvc )
__asm
{
;
fld theta;
fmul rad;
fstp theta; // theta = theta * rad : convert radians to degrees
;
fild i;
fstp freq; // freq = i : frequency of DFT at this point in time
;
}
*(cxx + i) = cv; // cxx[i] = cv : amplitude
*(th + i) = theta; // th[i] = theta : phase shift in degrees
*(fr + i) = freq; // fr[i] = freq : frequency
__asm
{
;
inc i; // i = i + 1
mov eax,i; // eax = i
cmp eax,n; // eax ? n
jl Loop1; // jump to Loop1 label if j < n
;
}
}
void DFT32(float* dxx, float* cxx, float* th, float* fr, int cnt)
{
float dxe, cxe, theta, freq, rad, pi2, dne, dnx, cv, cvc;
int i, j, n;
n = cnt; // number of data items
rad = 360.0; // degrees in one cycle
__asm
{
;
finit; // init FPU
;
fldpi;
fldpi;
fadd;
fstp pi2; // pi2 = pi + pi : 2 * pi
;
fld rad;
fdiv pi2;
fstp rad; // rad = 360.0 / pi2 : rad = 57.3 deg/rad
;
fld pi2;
fchs; // change sign pi2 = -1.0 * pi2
fidiv n;
fstp pi2; // pi2 = (-1.0 * 2.0 * pi) / n
;
mov i,0; // i = 0 : outer loop counter
;
Loop1: ; // Top of outer Loop1
;
fild i;
fmul pi2;
fstp dnx; // dnx = pi2 * i
;
}
dxe = *(dxx + i); // dxe = dxx[i]
__asm
{
;
fldz;
fstp cv; // cv = 0 : power series accumulator
;
mov j,0; // j = 0 : inner loop counter
;
Loop2: ; // Top of inner Loop2
;
fild j;
fmul dnx;
fstp dne; // dne = dnx * j : dne = ( -1.0 * 2.0 * pi * i * j ) / n
;
}
cxe = exp(dne); // cxe = e^dne : c++ call ( math.h )
__asm
{
;
fld cxe;
fmul dxe;
fadd cv;
fstp cv; // cv = cv + ( cxe * dxe ) : calc power series
;
inc j; // j = j + 1 : increment inner loop counter
mov eax,j; // eax = j : cpu register eax = j
cmp eax,n; // eax ? n : n compared to cpu register eax
jl Loop2; // j < n ? if so jump to Loop2 label
;
fld cv;
fcos;
fmul cv;
fstp cvc; // cvc = (cos( cv ) * cv) : calc adjacent side
;
fld cvc;
fdiv cv;
fstp cvc; // cvc = cvc / cv : cvc = adjacent / hypoteneus
}
theta = acos(cvc); // theta = Acos( cvc )
__asm
{
;
fld theta;
fmul rad;
fstp theta; // theta = theta * rad : convert radians to degrees
;
fild i;
fstp freq; // freq = i : frequency of DFT at this point in time
;
}
*(cxx + i) = cv; // cxx[i] = cv : amplitude
*(th + i) = theta; // th[i] = theta : phase shift in degrees
*(fr + i) = freq; // fr[i] = freq : frequency
__asm
{
;
inc i; // i = i + 1
mov eax,i; // eax = i
cmp eax,n; // eax ? n
jl Loop1; // jump to Loop1 label if j < n
;
}