Skip to content Skip to navigation

Connexions

You are here: Home » Content » Discrete Fourier Transform

Navigation

Recently Viewed

This feature requires Javascript to be enabled.
 

Discrete Fourier Transform

Module by: Marc Balke. E-mail the author

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

;

}

Content actions

Download module as:

Add module to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

Definition of a lens

Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

Who can create a lens?

Any individual member, a community, or a respected organization.

What are tags? tag icon

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks