Hello jcl and others

I was wondering if you could help me translate John Ehlers' Autocorrelation Periodogram into Lite-C. The indicator is found in the book 'Cycle Analytics for Traders'. The indicator uses an autocorrleation process to calculate the dominant period, and I believe will do so with less lag than Zorro's DominantPeriod function.

I will firstly describe the indicator's calculation, then I will post the TradeStation code that appears in the book, and finally I will post my attempt and describe where I am stuck.

The indicator is calculated by firstly pre-filtering the data using the roofing filter. I chose low and high cutoff periods of 10 and 50. Then, we correlate the output of the filter with the output delayed by each of a series of lag values, using Person correlation. Ehlers autocorrelates using 10 to 48 bars of lag, and I will try to do the same. Each value of lag is assigned a correlation coefficient between -1 (perfect anticorellation) and 1 (perfect correlation).

After calculating the autocorrelation for each value of lag, Ehlers extracts the cyclic information using a discrete Fourier transform of the autocorrelation results. This is accomplished by correlating the autocorrelation at each value of lag with the cosine and sine of each period of interest. The sum of the squares of each of these values represents relative power at each period from trigonometry. An EMA is used to smooth the power measurement at each period.

Next, the automatic gain control function is used to normalize the spectral components. Finally, the dominant cycle is extracted using a centre of gravity algorithm.

Ehlers uses a heat map display to show the strength of each period at each bar. I don't think this is possible in Zorro, but that shouldn't stop us from extracting the dominant cycle using this method.

I've attached the pages from 'Cycle Analytics for Traders' that describe the indicator, including the TradeStation code. The TradeStation code is also below:

Code:
{
Autocorrelation Periodogram
© 2013	John F. Ehlers
}

Vars:
AvgLength(3), M(0),
N(0),	 
X(0),
Y(0),
alpha1(0), HP(0),
a1(0),
b1(0),
c1(0),
c2(0),
c3(0),
Filt(0),
Lag(0),
count(0),
Sx(0),
Sy(0),
Sxx(0),
Syy(0),
Sxy(0),
Period(0), Sp(0),
Spx(0),
MaxPwr(0), DominantCycle(0), Color1(0), Color2(0), Color3(0);

Arrays:
Corr[48](0),
CosinePart[48](0), SinePart[48](0),
SqSum[48](0),
R[48, 2](0),
Pwr[48](0);

//Highpass filter cyclic components whose periods are shorter than 48 bars
alpha1 =  (Cosine(.707*360 / 48) +  Sine (.707*360 / 48) - 1) / Cosine(.707*360 / 48);
HP =  (1 - alpha1 / 2)*(1 - alpha1 / 2)*(Close - 2*Close[1] +
Close[2]) +  2*(1 - alpha1)*HP[1] - (1 - alpha1)* (1 - alpha1)*HP[2];
//Smooth with a Super Smoother Filter from equation 3-3
 

a1 =  expvalue(-1.414*3.14159 / 10); b1 =  2*a1*Cosine(1.414*180 / 10); c2 =  b1;
c3 =  -a1*a1;
c1 =  1 - c2 - c3;
Filt =  c1*(HP +  HP[1]) / 2 +  c2*Filt[1] +  c3*Filt[2];
 

//Pearson correlation for each value of lag For Lag =  0 to 48 Begin
//Set the averaging length as M M =  AvgLength;
If AvgLength =  0 Then M =  Lag;
Sx =  0;
Sy =  0;
Sxx =  0;
Syy =  0;
Sxy =  0;
For count =  0 to M - 1 Begin X =  Filt[count];
Y =  Filt[Lag +  count]; Sx =  Sx +  X;
Sy =  Sy +  Y;
Sxx = Sxx + X*X; Sxy = Sxy + X*Y; Syy =  Syy +  Y*Y;
End;
If (M*Sxx - Sx*Sx)*(M*Syy - Sy*Sy) > 0 Then Corr[Lag] =
(M*Sxy - Sx*Sy)/SquareRoot((M*Sxx - Sx*Sx)*(M*Syy - Sy*Sy)); End;
For Period =  10 to 48 Begin CosinePart[Period] =  0;
SinePart[Period] =  0; For N =  3 to 48 Begin
CosinePart[Period] =   CosinePart[Period] +
Corr[N]*Cosine(370*N / Period);
SinePart[Period] = SinePart[Period] + Corr[N]*Sine(370*N / Period);
End;
SqSum[Period] =    CosinePart[Period]*CosinePart[Period] +
SinePart[Period]*SinePart[Period];
End;	(Continued )
 
For Period =  10 to 48 Begin R[Period, 2] =  R[Period, 1];
R[Period, 1] =  .2*SqSum[Period]*SqSum[Period] +
.8*R[Period, 2]; End;

//Find Maximum Power Level for Normalization MaxPwr =  .995*MaxPwr;
For Period =  10 to 48 Begin
If R[Period, 1] > MaxPwr Then MaxPwr =  R[Period, 1]; End;
For Period =  3 to 48 Begin
Pwr[Period] =  R[Period, 1] / MaxPwr; End;
 
//Compute the dominant cycle using the CG of the spectrum Spx =  0;
Sp =  0;
For Period =  10 to 48 Begin
If Pwr[Period] >=  .5 Then Begin Spx =  Spx +  Period*Pwr[Period]; Sp =  Sp +  Pwr[Period];
End; End;
If Sp <> 0 Then DominantCycle =  Spx / Sp;
Plot2(DominantCycle, “DC”, RGB(0, 0, 255), 0, 2);

{
//Increase Display Resolution by raising the NormPwr to a higher mathematical power (optional)
For Period =  10 to 48 Begin Pwr[Period] =  Power(Pwr[Period], 2);
End;
}


And here is my so far incomplete attempt:
Code:
function run()
{
//=====John Ehlers' Autocorrelation Periodogram=====

//Roofing Filter used on prices
vars RoofPrice = series(Roof(Price, 10, 40)); //standard roofing filter
vars RoofAGC = series(AGC(RoofPrice, 1)); //output normalized using AGC algorithm

//Ehlers Autocorrelation Periodogram

var corr[48]; //an array that will store the correlation coefficient for each value of lag from 1 to 48
for (i = 1; i < 49; i++)
	corr[i] = Correlation(RoofAGC, RoofAGC+i, i); //assign a correlation coefficent to each value of lag. 

//discrete Fourier transform of autcorrelation results
}


As you can see, I'm stuck on how to do the discrete Fourier transform. I just don't understand how to translate the TradeStation code for this part into Lite-C. I'm also unsure of how to best handle the arrays. Any advice or tips would be really appreciated, and hopefully enable me to translate the remaining TradeStation code.

Thanks in advance!

Attached Files
AutocorrelationPeriodogram.rar (8 downloads)
www.robotwealth.com