Fast fourier transform FFT

hi, anyone use the FFT function successfully yet?

what I was trying is the calculate the FFT of the last 2^5 data points and then reconstruct the time series from the frequency spectrum (power spectrum).

I was not successful yet. Hope somebody figured that out. Here is what I was playing with.

My aim is to calculate the frequency spectrum, do some filtering and see if I can extrapolate the reconstructed filtered time series. This is how frequency filtering is often done, calculate the FFT, filter and then calculate the IFFT. But I am mainly interested in extrapolation eventually. See how that looks like

len = 2^5;
SetBarsRequired( len, len );

SetChartBkColor( colorBlack );
SetChartOptions( 1, chartShowDates, chartGridMiddle, 0, 0, 0 );
SetBarFillColor( IIf( C > O, ColorRGB( 0, 75, 0 ), IIf( C <= O, ColorRGB( 75, 0, 0 ), colorLightGrey ) ) );
Plot( C, "", IIf( C > O, ColorRGB( 0, 255, 0 ), IIf( C <= O, ColorRGB( 255, 0, 0 ), colorLightGrey ) ), 64, Null, Null, 0, 0, 1 );

ffc = FFT( C, len );

for( i = 0; i < Len - 1; i = i + 2 )
{
    amp[ i ] = amp[ i + 1 ] = sqrt( ffc[ i ] ^ 2  +  ffc[ i + 1 ] ^ 2 );
    phase[ i ] = phase[ i + 1 ] = atan2( ffc[ i + 1], ffc[ i ] );
}

ff = 0;

for( k = 0; k < len; k = k + 2 )
{
    for( i = 0; i < len; i++ )
    {
        ff[BarCount - len + i] = ff[BarCount - len + i] + amp[k] * cos( ( 2 * 3.1415926 * i / len ) + phase[k] );
    }
}

ff = IIf( ff, ff, Null );
Plot(ff,"",colorWhite,1);

2 Likes

@empottasch your interests and knowledge seem encyclopedic! I have no first hand experience with FFT but a fellow named Dennis Meyers used to publish a lot of trading related FFT papers in the 1990’s and early 2000’s that may be of interest to you.

http://meyersanalytics.com/papers.html

Good luck.

1 Like

@empottasch and I don’t know if you have seen this afl,

http://www.amibroker.com/members/traders/01-2007.html

It should be something like this.
There is still a problem of scaling but the signal shape is good. Will try to correct it later.

Best regards.

len = 128;

// DFT
ffc = FFT( C, len );

// Real and Imaginary Parts in Frequency Domain
Re= Im = 0;
for( k = 0; k < (len/2) + 1; k ++ )
{
fRe[k] = ffc[ 2*k ];
fIm[k] = ffc[2*k+1];
}

// Real and Imaginary Parts in Time Domain
tRe  = 2*fRe/len;
tRe[0]  = fRe[k]/len;
tRe[len/2]  = fRe[len/2]/len;
tIm  = -2*fIm/len;

// Inverse DFT
ff = 0;
for( i = 0; i < len; i++)
{
	ff[BarCount - len + i] = 0;
    for( k = 0; k < (len/2) + 1; k ++ )
    {
        ff[BarCount - len + i] = ff[BarCount - len + i] + tRe[k] * cos(  2 * 3.1415926 * k * i / len ) + tIm[k]*sin(  2 * 3.1415926 * k * i / len );
    }
}

ff = IIf(ff != 0 , ff, Null );

Plot( C, "", colorblue);
Plot(ff,"Reconstructed Signal",colorblack, styleOwnScale);

With scaling correction.

len = 128;


// DFT
ffc = FFT( C, len );

// Real and Imaginary Parts in Frequency Domain
fRe= fIm = 0;
for( k = 0; k < (len/2) + 1; k ++ )
{
fRe[k] = ffc[ 2*k ];
fIm[k] = ffc[2*k+1];
}

// Real and Imaginary Parts in Time Domain
tRe  = 2*fRe/len;
tRe[0]  = fRe[0]/len;
tRe[len/2]  = fRe[len/2]/len;
tIm  = -2*fIm/len;

// Inverse DFT
ff = 0;
for( i = 0; i < len; i++)
{
	ff[BarCount - len + i] = 0;
    for( k = 0; k < (len/2) + 1; k ++ )
    {
        ff[BarCount - len + i] = ff[BarCount - len + i] + tRe[k] * cos(  2 * 3.1415926 * k * i / len ) + tIm[k]*sin(  2 * 3.1415926 * k * i / len );
    }
}

ff = IIf(ff != 0 , ff, Null );

Plot(ff,"Reconstructed Signal",colorblack);
Plot( C, "", colorblue);

3 Likes

thanks for the replies, guys. Great work Waz, seems to work. I will check it out a bit and will back on this later

yes I’ve seen that site before but forgot about it. Will check it out again, thanks

the idea was that I use the first few frequency bins to make some crude predictions. But I am not sure if it is possible to implement this correctly.

I have some code below that show the basic idea. With the mouse you can select a date on the chart and the 128 data points before that date will be used to calculate the frequency spectrum, the data after that date should be “predicted”.

I do not use all the frequency bins, in this example only the first 6.

But what I calculate now is just a repeat of the original curve. But post it anyway, maybe somebody has an idea

bi = BarIndex();
fvb = FirstVisibleValue( bi );
lvb = LastVisibleValue( bi );

dn = DateTime();
sd = SelectedValue( dn );
start = dn == sd;
stidx = LastValue( ValueWhen( start, BarIndex() ) );

shift =  BarCount - 1 - stidx;

len = 2 ^ 7; //128;
data = ref( C, -( shift ) );

// DFT
ffc = FFT( data, len );

// Real and Imaginary Parts in Frequency Domain
fRe = fIm = 0;

for( k = 0; k < ( len / 2 ) + 1; k ++ )
{
    fRe[k] = ffc[ 2 * k ];
    fIm[k] = ffc[2 * k + 1];
}

// Real and Imaginary Parts in Time Domain
tRe  = 2 * fRe / len;
tRe[0]  = fRe[0] / len;
tRe[len / 2]  = fRe[len / 2] / len;
tIm  = -2 * fIm / len;

// Inverse DFT
ff = 0;

for( i = 0; i < len + shift; i++ )
{
    //for( k = 0; k < ( len / 2 ) + 1; k ++ )
    for( k = 0; k <= 5; k ++ )
    {
        ff[BarCount - len + i - shift] = ff[BarCount - len + i - shift] + tRe[k] * cos( 2 * 3.1415926 * k * i / len ) + tIm[k] * sin( 2 * 3.1415926 * k * i / len );
    }
}

ff = IIf( ff != 0 , ff, Null );

SetChartBkColor( colorBlack );
SetChartOptions( 1, chartShowDates, chartGridMiddle, 0, 0, 0 );
SetBarFillColor( IIf( C > O, ColorRGB( 0, 75, 0 ), IIf( C <= O, ColorRGB( 75, 0, 0 ), colorLightGrey ) ) );
Plot( C, "", IIf( C > O, ColorRGB( 0, 255, 0 ), IIf( C <= O, ColorRGB( 255, 0, 0 ), colorLightGrey ) ), 64, Null, Null, 0, 0, 1 );

Plot( ff, "Reconstructed Signal", colorWhite, styleLine | styleNoRescale, Null, Null, 0, 0, 2 );
Plot( bi == stidx, "", colorDarkBlue, styleHistogram | styleOwnScale | styleNoLabel | styleNoRescale, 0, 1, 0, -2, 10 );

Actually, the FFT bins do not have any physical meaning, especially in the time domain. Therefore, removing some of them will do some filtering but you have no idea about what have been filtered. Filtering should be done on an amplitude spectrum computed in the frequency domain.
The procedure is a bit more complicated:
1 - preprocess the data to remove mean and trend otherwise the spectrum will be useless
2 - compute the FFT on the detrended data
3 - compute amplitudes and phases
4 - perform the filtering on the amplitude spectrum (for example: identify maximum amplitude and put to 0 all amplitudes smaller than 10% * max)
5 - compute real and imaginary parts from the filtered spectrum
6 - perform an inverse FFT

2 Likes

hi Waz,

well the lower bins contain the low frequencies. The first bin is very large because it contains most of the power since that is the large trend. So I agree one should detrend if you want to only observe higher frequency waves in the data. But this was not my intention.

Actually I started all this because I was working on this TB-F code and I needed a very good top or bottom estimate, but I will use something else for that.

But I will work a bit more on this so you can toggle to the power spectrum view. How is that calculated in your opinion? Like this?

Power = tRe^2 + tIm^2;

In my opinion detrending is not necessary. It just removes the lower frequencies but you can just take them out the the power spectrum or FFT also (at least that is my understanding)

Still it would maybe be interesting to look if a certain periodicity at higher frequency also has an unusual amount of power

some form of detrending is probably a better way to go. Will let it sink in a bit. But I doubt frequency analysis will ever be of any use.

I added some stuff. In the param window you can toggle to the power spectrum. The bin start position allows you to remove the first bins so you can have a better look at the higher frequencies.

Also, when you are in the frequency domain and click on the chart it will still move the chosen data array in the time domain so when in the frequency domain you should not do that else the spectrum will change.

Frequency or the wavelength is shown as a mouse tip. Made this for daily data

So now will go do some other stuff again first :smile:

SetBarsRequired( sbrAll, sbrAll );

bi = BarIndex();
fvb = FirstVisibleValue( bi );
lvb = LastVisibleValue( bi );

dn = DateTime();
sd = SelectedValue( dn );
start = dn == sd;
stidx = LastValue( ValueWhen( start, BarIndex() ) );
shift =  Max( 0, BarCount - 1 - stidx );

len = Param( "Number of Data", 512, 8, 2048, 2 );
bins = Param( "Number of Frequency Bins to show (Time Domain)", 20, 1, len / 2 + 1, 1 );
view = ParamToggle( "Display", "Frequency Domain| Time domain", 1 );
st = Param( "Start bin position (Power Domain)", 3, 0, len / 2 + 1, 1 );

data = ref( C, -( shift ) );

// DFT
ffc = FFT( data, len );

// Real and Imaginary Parts in Frequency Domain
fRe = fIm = 0;

for( k = 0; k < ( len / 2 ) + 1; k++ )
{
    fRe[k] = ffc[ 2 * k ];
    fIm[k] = ffc[2 * k + 1];
}

// Real and Imaginary Parts in Time Domain
tRe  = 2 * fRe / len;
tRe[0]  = fRe[0] / len;
tRe[len / 2]  = fRe[len / 2] / len;
tIm  = -2 * fIm / len;

// Inverse DFT
if( view )
{
    ff = 0;

    for( i = 0; i < len + shift; i++ )
    {
        //for( k = 0; k < ( len / 2 ) + 1; k ++ )
        for( k = 0; k <= bins; k ++ )
        {
            ff[BarCount - len + i - shift] = ff[BarCount - len + i - shift] + tRe[k] * cos( 2 * 3.1415926 * k * i / len ) + tIm[k] * sin( 2 * 3.1415926 * k * i / len );
        }
    }

    ff = IIf( ff != 0 , ff, Null );

    SetChartBkColor( colorBlack );
    SetChartOptions( 1, chartShowDates, chartGridMiddle, 0, 0, 0 );

    SetBarFillColor( IIf( C > O, ColorRGB( 0, 75, 0 ), IIf( C <= O, ColorRGB( 75, 0, 0 ), colorLightGrey ) ) );
    Plot( C, "", IIf( C > O, ColorRGB( 0, 255, 0 ), IIf( C <= O, ColorRGB( 255, 0, 0 ), colorLightGrey ) ), 64, Null, Null, 0, 0, 1 );

    Plot( ff, "Reconstructed Signal", colorWhite, styleLine | styleNoRescale, Null, Null, 0, 0, 2 );
    Plot( bi == stidx, "", colorDarkBlue, styleHistogram | styleOwnScale | styleNoLabel | styleNoRescale, 0, 1, 0, -2, 10 );
}
else
{
    // power spectrum calculation
    power = tRe ^ 2 + tIm ^ 2;

    GfxSetOverlayMode( 2 );
    GfxSetCoordsMode( 0 );
    GfxSelectPen( ColorRGB( 0, 255, 255 ), 2 );

    pxwidth = Status( "pxwidth" );
    pxheight = Status( "pxheight" );
    margin = int( pxwidth / 20 );

    // chart border
    GfxMoveTo( margin, margin );
    GfxLineTo( pxwidth - margin, margin );
    GfxLineTo( pxwidth - margin, pxheight - margin );
    GfxLineTo( margin, pxheight - margin );
    GfxLineTo( margin, margin );

    GfxSelectFont( "Tahoma", 15, 700, False, False, 0 );
    GfxTextOut( "Power Spectrum" , pxwidth / 2, margin / 2 );
    GfxTextOut( "Frequency" , pxwidth / 2, pxheight - margin / 2 );
    GfxSelectFont( "Tahoma", 15, 700, False, False, 900 );
    GfxTextOut( "Power" , margin / 2, pxheight / 2 );

    maxpower = 0;

    for( i = st; i < ( len / 2 ) + 1; i++ )
    {
        if( power[i] > maxpower )
        {
            maxpower = power[i];
        }
    }

    xbin = ( pxwidth - 2 * margin ) / ( ( len / 2 ) + 1 );
    ybin = ( pxheight - 2 * margin ) / maxpower;

    if( st > 0 )
    {
        GfxMoveTo( margin, pxheight - margin );

        for( i = 0; i < st; i++ )
        {
            GfxLineTo( margin + i * xbin, pxheight - margin );
        }
    }

    if( st == 0 ) GfxMoveTo( margin, pxheight - margin - power[st]*ybin );

    for( i = st + 1; i < ( len / 2 ) + 1; i++ )
    {
        GfxLineTo( margin + i * xbin, pxheight - margin - power[i]*ybin );
    }

    // axis
    nyquist = 1 / ( 2 * Interval() );
    RequestMouseMoveRefresh();
    // retrieve co-ordinates in pixel units
    px = GetCursorXPosition( 1 );
    py = GetCursorYPosition( 1 );
    GfxSelectFont( "Tahoma", 10, 700, False, False, 0 );

    if( px > margin AND px < pxwidth - margin )
    {
        pxr = ( px - margin ) / ( pxwidth - 2 * margin ) ;

        freq = pxr * nyquist;
        wavelength = ( 1 / freq ) / ( 24.*60 * 60 );

        GfxSetTextColor( ColorRGB( 0, 0, 0 ) );
        GfxSetBkColor( ColorRGB( 255, 255, 0 ) );

        if( px < pxwidth / 2 )
            GfxTextOut( "Frequency: " + pxr * nyquist + " Hz, Wavelength: " + wavelength + " days" , px + 15, py );

        if( px > pxwidth / 2 )
            GfxTextOut( "Frequency: " + pxr * nyquist + " Hz, Wavelength: " + wavelength + " days" , px - 400, py );
    }

    GfxSelectFont( "Tahoma", 15, 700, False, False, 0 );
    GfxSetTextColor( ColorRGB( 255, 0, 0 ) );
    GfxSetBkColor( ColorRGB( 255, 255, 255 ) );
    GfxTextOut( "Do not mouseclick on chart => will change spectrum" , margin, margin / 2 );
}

Here is my interpretation with detrending and spectrum.
Result looks correct.
The spectrum is plotted in reverse order.
The red line on the spectrum is the amplitude cut-off.

I still have problem if I want to increase the length to more than 200 bars. It gives an error message when computing fRe and fIm from the fft bins. It gives an error 10: trying to access a non-existing 200-th elements of the array. I cannot figure out why.

image

th = Param("Threshold for Filtering", 10, 5, 50, 1);
len = 128;
pi = 3.1415926;

// Detrend Data
x0 = BarCount - len - 1; 
a = LastValue( LinRegIntercept( C, len) ); 
b = LastValue( LinRegSlope( C, len) ); 
lr = a + b * (BarIndex() - x0); 
data = C - lr; 

// DFT
ffc = FFT( data, len );

// Real and Imaginary Parts in Frequency Domain
fRe= fIm = 0;
for( k = 0; k < len/2; k ++ )
{
fRe[k] = ffc[2 * k];
fIm[k] = ffc[2 * k + 1];
}

// Compute Amplitude and Phase
amp = sqrt(fRe^2 + fIm^2);
ph = atan2(fIm, fRe);

// Plot Amplitude Spectrum
spectrum = Reverse(amp);

// Filter the amplitude spectrum
maxamp = Highest(amp);
amp = IIf(amp < maxamp*th/100, 0, amp);

// Compute real and imaginary parts of the filtered spectrum
fRe = amp * cos(ph);
fIm = amp * sin(ph);


// Real and Imaginary Parts in Time Domain
tRe  = 2*fRe/len;
tRe[0]  = fRe[0]/len;
tRe[(len/2) - 1]  = fRe[(len/2) - 1]/len;
tIm  = -2*fIm/len;

// Inverse DFT
ff = 0;
for( i = 0; i < len; i++)
{
	ff[BarCount - len + i] = 0;
    for( k = 0; k < len/2; k ++ )
    {
        ff[BarCount - len + i] = ff[BarCount - len + i] + tRe[k] * cos(  2 * pi * k * i / len ) + tIm[k]*sin(  2 * pi * k * i / len );
    }
}

ff = IIf(ff != 0 , ff, Null );

// Retrend data
ff = ff + lr;

plotchoice = ParamList("Plot : ", "Price & Filter|Spectrum");

if( plotchoice == "Price & Filter" )
	{ 
		SetChartBkColor(colorBlack);
		Plot(ff,"Reconstructed Signal",colorLightBlue, styleThick);
		Plot( C, "", colorYellow, styleBar);
	}

if( plotchoice == "Spectrum" )
	{ 
		SetChartBkColor(colorBlack);
		plot(spectrum, "Amplitude Spectrum", colorYellow, styleHistogram|styleThick);
		Plot(th*maxamp/100, "", colorRed, styleDashed);

	}



http://www.amibroker.com/kb/2014/09/22/do-not-make-assumptions-on-number-of-bars/

i think we are more or less talking about the same thing. I think the way you have originally split the ffc array up in a real and imaginary part is correct.

in the code below I use a simple synthetic data array where I add 2 functions with different frequency and different amplitude, the first a cosine with a wavelength of 60 days and amplitude 10 and the second a sinus with wavelength of 10 days and a amplitude of 3.

When you switch to the frequency domain you can see 2 distinct peaks in the power spectrum, 1 at approx 60 days and 1 at approx 10 days wavelength. So that is correct. Also their relative power seems correct, which should be 10^2 versus 3^2

Now if you change the “Number of Freq bins” in the param window you can cut off the higher frequencies by putting the number at 10 or so

It all seems to function correctly imo

SetBarsRequired( sbrAll, sbrAll );

pi = 3.1415926;
bi = BarIndex();

//len = Param( "Number of Data", 256, 8, 4096, 2 );
len = 128;
bins = Param( "Number of Frequency Bins to show (Time Domain)", 20, 1, len / 2, 1 );
view = ParamToggle( "Display", "Frequency Domain| Time domain", 1 );
st = Param( "Start bin position (Power Domain)", 0, 0, len / 2, 1 );

data = 10 * cos( 2 * pi * bi / 60 ) + 3 * sin( 2 * pi * bi / 10 );

// DFT
ffc = FFT( data, len );

// Real and Imaginary Parts in Frequency Domain
fRe = fIm = 0;

for( k = 0; k < ( len / 2 ) ; k++ )
{
    fRe[k] = ffc[ 2 * k ];
    fIm[k] = ffc[2 * k + 1];
}

// Real and Imaginary Parts in Time Domain
tRe  = 2 * fRe / len;
tRe[0]  = fRe[0] / len;
tRe[len / 2]  = fRe[len / 2] / len;
tIm  = -2 * fIm / len;

// Inverse DFT
if( view )
{
    ff = 0;

    for( i = 0; i < len; i++ )
    {
        for( k = 0; k < bins; k++ )
        {
            ff[BarCount - len + i] = ff[BarCount - len + i] + tRe[k] * cos( 2 * pi * k * i / len ) + tIm[k] * sin( 2 * pi * k * i / len );
        }
    }

    ff = IIf( ff != 0 , ff, Null );
	
    SetChartBkColor( colorBlack );
    SetChartOptions( 1, chartShowDates, chartGridMiddle, 0, 0, 0 );
    Plot( data, "data", colorBlue, styleLine, Null, Null, 0, 0, 2 );
    Plot( ff, "Reconstructed Signal", colorWhite, styleLine | styleNoRescale, Null, Null, 0, 0, 8 );
}
else
{
    // power spectrum calculation
    power = tRe ^ 2 + tIm ^ 2;

    GfxSetOverlayMode( 2 );
    GfxSetCoordsMode( 0 );
    GfxSelectPen( ColorRGB( 0, 255, 255 ), 2 );

    pxwidth = Status( "pxwidth" );
    pxheight = Status( "pxheight" );
    margin = int( pxwidth / 20 );

    // chart border
    GfxMoveTo( margin, margin );
    GfxLineTo( pxwidth - margin, margin );
    GfxLineTo( pxwidth - margin, pxheight - margin );
    GfxLineTo( margin, pxheight - margin );
    GfxLineTo( margin, margin );

    GfxSelectFont( "Tahoma", 15, 700, False, False, 0 );
    GfxTextOut( "Power Spectrum" , pxwidth / 2, margin / 2 );
    GfxTextOut( "Frequency" , pxwidth / 2, pxheight - margin / 2 );
    GfxSelectFont( "Tahoma", 15, 700, False, False, 900 );
    GfxTextOut( "Power" , margin / 2, pxheight / 2 );

    maxpower = 0;

    for( i = st; i < ( len / 2 ) + 1; i++ )
    {
        if( power[i] > maxpower )
        {
            maxpower = power[i];
        }
    }

    xbin = ( pxwidth - 2 * margin ) / ( ( len / 2 ) + 1 );
    ybin = ( pxheight - 2 * margin ) / maxpower;

    GfxSelectPen( ColorRGB( 0, 255, 255 ), 6, 0 );

    for( i = st; i < ( len / 2 ) + 1; i++ )
    {
        GfxMoveTo( margin + i * xbin, pxheight - margin );
        GfxLineTo( margin + i * xbin, pxheight - margin - power[i]*ybin );
    }

    // axis
    nyquist = 1 / ( 2 * Interval() );
    RequestMouseMoveRefresh();
    // retrieve co-ordinates in pixel units
    px = GetCursorXPosition( 1 );
    py = GetCursorYPosition( 1 );
    GfxSelectFont( "Tahoma", 10, 700, False, False, 0 );

    if( px > margin AND px < pxwidth - margin )
    {
        pxr = ( px - margin ) / ( pxwidth - 2 * margin ) ;

        freq = pxr * nyquist;
        wavelength = ( 1 / freq ) / ( 24.*60 * 60 );

        GfxSetTextColor( ColorRGB( 0, 0, 0 ) );
        GfxSetBkColor( ColorRGB( 255, 255, 0 ) );

        if( px < pxwidth / 2 )
            GfxTextOut( "Frequency: " + pxr * nyquist + " Hz, Wavelength: " + wavelength + " days" , px + 15, py );

        if( px > pxwidth / 2 )
            GfxTextOut( "Frequency: " + pxr * nyquist + " Hz, Wavelength: " + wavelength + " days" , px - 400, py );
    }

    GfxSelectFont( "Tahoma", 15, 700, False, False, 0 );
    GfxSetTextColor( ColorRGB( 255, 0, 0 ) );
    GfxSetBkColor( ColorRGB( 255, 255, 255 ) );
    GfxTextOut( "Do not mouseclick on chart => will change spectrum" , margin, margin / 2 );
}

with respect to the post from Tomasz I added some stuff that avoids error messages. Also made little change in the Frequency domain chart.

I am using this as my param settings:

len = Param( "Number of Data", 512, 8, 4096, 2 );
bins = Param( "Number of Frequency Bins to show (Time Domain)", 20, 1, len / 2, 1 );
view = ParamToggle( "Display", "Frequency Domain| Time domain", 1 );
st = Param( "Start bin position (Power Domain)", 0, 0, len / 2, 1 );

I have a question if this is correct since default/min/max/step parameters have to be constant numbers. It seems to work correctly but I wonder if this is allowed.

adjusted code below

SetBarsRequired( sbrAll, sbrAll );

pi = 3.1415926;
bi = BarIndex();

len = Param( "Number of Data", 512, 8, 4096, 2 );
bins = Param( "Number of Frequency Bins to show (Time Domain)", 20, 1, len / 2, 1 );
view = ParamToggle( "Display", "Frequency Domain| Time domain", 1 );
st = Param( "Start bin position (Power Domain)", 0, 0, len / 2, 1 );

data = 10 * cos( 2 * pi * bi / 60 ) + 3 * sin( 2 * pi * bi / 10 );

// DFT
ffc = FFT( data, len );

// Real and Imaginary Parts in Frequency Domain
fRe = fIm = 0;

if( len / 2 < BarCount )
{
    for( k = 0; k < ( len / 2 ) ; k++ )
    {
        fRe[k] = ffc[ 2 * k ];
        fIm[k] = ffc[2 * k + 1];
    }

    // Real and Imaginary Parts in Time Domain
    tRe  = 2 * fRe / len;
    tRe[0]  = fRe[0] / len;
    tRe[len / 2]  = fRe[len / 2] / len;
    tIm  = -2 * fIm / len;
}

// Inverse DFT
if( view )
{
    ff = 0;

    if( len < BarCount )
    {
        for( i = 0; i < len; i++ )
        {
            for( k = 0; k < bins; k++ )
            {
                ff[BarCount - len + i] = ff[BarCount - len + i] + tRe[k] * cos( 2 * pi * k * i / len ) + tIm[k] * sin( 2 * pi * k * i / len );
            }
        }
    }

    ff = IIf( ff != 0 , ff, Null );

    SetChartBkColor( colorBlack );
    SetChartOptions( 1, chartShowDates, chartGridMiddle, 0, 0, 0 );
    Plot( data, "data", colorBlue, styleLine, Null, Null, 0, 0, 2 );
    Plot( ff, "Reconstructed Signal", colorWhite, styleLine | styleNoRescale, Null, Null, 0, 0, 8 );
}
else
{
    if( len < BarCount )
    {
        // power spectrum calculation
        power = tRe ^ 2 + tIm ^ 2;

        GfxSetOverlayMode( 2 );
        GfxSetCoordsMode( 0 );
        GfxSelectPen( ColorRGB( 0, 255, 255 ), 2 );

        pxwidth = Status( "pxwidth" );
        pxheight = Status( "pxheight" );
        margin = int( pxwidth / 20 );

        // chart border
        GfxMoveTo( margin, margin );
        GfxLineTo( pxwidth - margin, margin );
        GfxLineTo( pxwidth - margin, pxheight - margin );
        GfxLineTo( margin, pxheight - margin );
        GfxLineTo( margin, margin );

        GfxSelectFont( "Tahoma", 15, 700, False, False, 0 );
        GfxTextOut( "Power Spectrum" , pxwidth / 2, margin / 2 );
        GfxTextOut( "Frequency" , pxwidth / 2, pxheight - margin / 2 );
        GfxSelectFont( "Tahoma", 15, 700, False, False, 900 );
        GfxTextOut( "Power" , margin / 2, pxheight / 2 );

        maxpower = 0;

        for( i = st; i < ( len / 2 ); i++ )
        {
            if( power[i] > maxpower )
            {
                maxpower = power[i];
            }
        }

        xbin = ( pxwidth - 2 * margin ) / ( ( len / 2 ) - 1 );
        ybin = ( pxheight - 2 * margin ) / maxpower;

        GfxSelectPen( ColorRGB( 0, 255, 255 ), 6, 0 );

        for( i = st; i < ( len / 2 ); i++ )
        {
            GfxMoveTo( margin + i * xbin, pxheight - margin );
            GfxLineTo( margin + i * xbin, pxheight - margin - power[i]*ybin );
        }

        // axis
        nyquist = 1 / ( 2 * Interval() );
        RequestMouseMoveRefresh();
        // retrieve co-ordinates in pixel units
        px = GetCursorXPosition( 1 );
        py = GetCursorYPosition( 1 );
        GfxSelectFont( "Tahoma", 10, 700, False, False, 0 );

        if( px > margin AND px < pxwidth - margin )
        {
            pxr = ( px - margin ) / ( pxwidth - 2 * margin ) ;

            freq = pxr * nyquist;
            wavelength = ( 1 / freq ) / ( 24.*60 * 60 );

            GfxSetTextColor( ColorRGB( 0, 0, 0 ) );
            GfxSetBkColor( ColorRGB( 255, 255, 0 ) );

            if( px < pxwidth / 2 )
                GfxTextOut( "Frequency: " + pxr * nyquist + " Hz, Wavelength: " + wavelength + " days" , px + 15, py );

            if( px > pxwidth / 2 )
                GfxTextOut( "Frequency: " + pxr * nyquist + " Hz, Wavelength: " + wavelength + " days" , px - 400, py );
        }

        GfxSelectFont( "Tahoma", 15, 700, False, False, 0 );
        GfxSetTextColor( ColorRGB( 255, 0, 0 ) );
        GfxSetBkColor( ColorRGB( 255, 255, 255 ) );
        GfxTextOut( "Do not mouseclick on chart => will change spectrum" , margin, margin / 2 );
    }
}

to answer by own question, no it is not correct to use

len = Param( "Number of Data", 512, 8, 4096, 2 );
bins = Param( "Number of Frequency Bins to show (Time Domain)", 20, 1, len / 2, 1 );
view = ParamToggle( "Display", "Frequency Domain| Time domain", 1 );
st = Param( "Start bin position (Power Domain)", 0, 0, len / 2, 1 );

this should be correct but if you change the value of len one should also do a “reset all” in the param window

len = 2048;
//len = Param( "Number of Data", 512, 8, 4096, 2 );
bins = Param( "Number of Frequency Bins to show (Time Domain)", 20, 1, len / 2, 1 );
view = ParamToggle( "Display", "Frequency Domain| Time domain", 1 );
st = Param( "Start bin position (Power Domain)", 0, 0, len / 2, 1 );

still did some work on this even though I suspect it will be worthless for trading.

I found some of the calculations Waz uses in his code in the DSP guide, chapter 8 ( http://www.dspguide.com/CH8.PDF ) page 152, 153. What I still do not understand is that in that same chapter 8 on page 147 they describe the data contains N samples while the Re and Im part in the frequency domain contain N/2+1 samples. This is however not what we use. We use N Re and Im pairs to reconstruct the signal. In the DSPguide (page 148) he writes: “That is, N points in the time domain corresponds to N/2+1 points in the frequency domain (not N/2 points). Forgetting about this extra point is a common bug in DFT programs.”

We however use N/2 points for signal reconstruction. Still it seems to look good.

added high and low frequency cutoff filter although not sure if my filtering procedure is mathematically correct it seems to work. You can toggle to the Frequency domain and set the filter

/*
© AFL code by Waz and E.M.Pottasch, 11/2017
originated here: http://forum.amibroker.com/t/fast-fourier-transform-fft/3098/11

code is an attempt to explore the use of frequency analysis in trading

some reading: http://www.dspguide.com/  
http://www.dspguide.com/CH8.PDF
*/
SetBarsRequired( sbrAll, sbrAll );

pi = 3.1415926;
bi = BarIndex();

len = 256;
view = ParamToggle( "Display", "Frequency Domain| Time domain", 1 );
bins1 = Param( "Low Cutoff Frequency", 0, 0, len / 2, 1 );
bins2 = Param( "High Cutoff Frequency", len / 2, 0, len / 2, 1 );
st = Param( "Start bin position (Frequency Domain, Display purposes only)", 0, 0, len / 2, 1 );

data = 10 * cos( 2 * pi * bi / 60 ) +
       3 * sin( 2 * pi * bi / 10 ) +
       6 * sin( 2 * pi * bi / 5 + 1 / 4 * pi ) +
       4 * sin( 2 * pi * bi / 20 + 2 / 5 * pi );

// DFT
ffc = FFT( data, len );

// Real and Imaginary Parts in Frequency Domain (Amibroker manual)
fRe = fIm = 0;

if( len < BarCount )
{
    for( k = 0; k < ( len / 2 ); k++ )
    {
        fRe[k] = ffc[2 * k];
        fIm[k] = ffc[2 * k + 1];
    }

	// http://www.dspguide.com/CH8.PDF, page 153
    tRe  = 2 * fRe / len;
    tRe[0]  = fRe[0] / len;
    tRe[len / 2 - 1]  = fRe[len / 2 - 1] / len;
    tIm  = -2 * fIm / len;
}

// Inverse DFT
if( view )
{
    ff = 0;

    if( len < BarCount )
    {
        for( n = 0; n < len; n++ )
        {
            for( k = bins1; k < bins2; k++ )
            {
                // http://www.dspguide.com/CH8.PDF, page 152
                ff[BarCount - len + n] = ff[BarCount - len + n] + tRe[k] * cos( 2 * pi * k * n / len ) + tIm[k] * sin( 2 * pi * k * n / len );
            }
        }
    }

    ff = IIf( ff != 0 , ff, Null );

    SetChartBkColor( colorBlack );
    SetChartOptions( 1, chartShowDates, chartGridMiddle, 0, 0, 0 );
    Plot( data, "data", ColorRGB( 0, 0, 255 ), styleLine, Null, Null, 0, 0, 2 );
    Plot( ff, "Reconstructed Signal", ColorRGB( 0, 255, 0 ), styleLine | styleNoRescale, Null, Null, 0, 0, 8 );
}
else
{
    if( len < BarCount )
    {
        // power spectrum calculation
        power = tRe ^ 2 + tIm ^ 2;

        GfxSetOverlayMode( 2 );
        GfxSetCoordsMode( 0 );
        GfxSelectPen( ColorRGB( 90, 90, 90 ), 1, 0 );

        pxwidth = Status( "pxwidth" );
        pxheight = Status( "pxheight" );
        margin = int( pxwidth / 20 );

        // chart border
        GfxMoveTo( margin, margin );
        GfxLineTo( pxwidth - margin, margin );
        GfxLineTo( pxwidth - margin, pxheight - margin );
        GfxLineTo( margin, pxheight - margin );
        GfxLineTo( margin, margin );

        GfxSelectFont( "Tahoma", 15, 700, False, False, 0 );
        GfxTextOut( " Frequency " , pxwidth / 2, pxheight - margin / 2 );
        GfxSelectFont( "Tahoma", 15, 700, False, False, 900 );
        GfxTextOut( " Power " , margin / 2, pxheight / 2 );

        maxpower = 0;

        for( i = st; i < ( len / 2 ); i++ )
        {
            if( power[i] > maxpower )
            {
                maxpower = power[i];
            }
        }

        xbin = ( pxwidth - 2 * margin ) / ( ( len / 2 ) - 1 );
        ybin = ( pxheight - 2 * margin ) / maxpower;

        GfxSelectPen( ColorRGB( 0, 255, 255 ), 6, 0 );

        for( i = st; i < ( len / 2 ); i++ )
        {
            GfxMoveTo( margin + i * xbin, pxheight - margin );
            GfxLineTo( margin + i * xbin, pxheight - margin - power[i]*ybin );
        }

        nyquist = 1 / ( 2 * Interval() );
        RequestMouseMoveRefresh();
        px = GetCursorXPosition( 1 );
        py = GetCursorYPosition( 1 );
        GfxSelectFont( "Tahoma", 10, 700, False, False, 0 );

        if( px > margin AND px < pxwidth - margin )
        {
            pxr = ( px - margin ) / ( pxwidth - 2 * margin ) ;

            freq = pxr * nyquist;
            wavelength = ( 1 / freq ) / ( 24 * 60 * 60 );

            GfxSetTextColor( ColorRGB( 0, 0, 0 ) );
            GfxSetBkColor( ColorRGB( 255, 255, 0 ) );

            if( px < pxwidth / 2 )
                GfxTextOut( " Frequency: " + pxr * nyquist + " Hz, Wavelength: " + prec( wavelength, 2 ) + " days " , px + 15, py );

            if( px > pxwidth / 2 )
                GfxTextOut( " Frequency: " + pxr * nyquist + " Hz, Wavelength: " + Prec( wavelength, 2 ) + " days " , px - 400, py );
        }

        GfxSetTextColor( ColorRGB( 0, 0, 0 ) );
        GfxSetBkColor( ColorRGB( 255, 255, 255 ) );
        GfxSelectPen( ColorRGB( 255, 0, 0 ), 1, 0 );
        // low cutoff frequency
        GfxMoveTo( margin + bins1 * xbin, pxheight - margin );
        GfxLineTo( margin + bins1 * xbin, margin );
        pxr1 = ( ( margin + bins1 * xbin ) - margin ) / ( pxwidth - 2 * margin ) ;
        freq1 = pxr1 * nyquist;
        wavelength1 = ( 1 / freq1 ) / ( 24 * 60 * 60 );
        GfxSelectFont( "Tahoma", 9, 700, False, False, 900 );
        GfxTextOut( " Low Cutoff Frequency " + Prec( wavelength1, 2 ) + " Days ", margin + bins1 * xbin , pxheight / 3 );
        // high cutoff frequency
        GfxMoveTo( margin + ( bins2 - 1 ) * xbin, pxheight - margin );
        GfxLineTo( margin + ( bins2 - 1 ) * xbin, margin );
        pxr2 = ( ( margin + ( bins2 - 1 ) * xbin ) - margin ) / ( pxwidth - 2 * margin ) ;
        freq2 = pxr2 * nyquist;
        wavelength2 = ( 1 / freq2 ) / ( 24 * 60 * 60 );
        GfxSelectFont( "Tahoma", 9, 700, False, False, 900 );
        GfxTextOut( " High Cutoff Frequency " + Prec( wavelength2, 2 ) + " Days ", margin + ( bins2 - 1 ) * xbin , pxheight / 3 );
    }
}

Concerning the number of points, I do not think that it is a big issue. Anyway AB fft is computed using 2n points from 0 to 2n-1.
For your test data, I would be better to add trend and noise to the data.

Cheers.

not sure about these 2n points you mention since the code starts of with “len” points and the FFT returns also “len” points, which are len/2 pairs. So if you replace “len” with N then we have N data points that correspond to N/2 points in the frequency domain. And according to the theory that should be N/2+1. But all seems to look good.

That was my basic goal with this that we have a way to analyse the frequency domain and reconstruct the signal correctly. But I looked at a few power spectra of EOD data and there is very little significant peaks to be seen at the higher frequencies, just looks like noise to me. So I think it will be useless but will look if there is any information out there how others claim to use it.

With all “reconstruction” attempts you need to understand that Fourier transformation assumes that input signal is periodic. With market data you almost never have pure periodic signal and there are discontinuities at the beginning/end of analysis window, so you need to precondition the input using DC component removal, de-trending and windowing, see: https://dsp.stackexchange.com/questions/11312/why-should-one-use-windowing-functions-for-fft

2 Likes