/*
	Autocorrelation Periodogram
	2013 John F. Ehlers
*/

SetBarsRequired(sbrAll);

LPeriod = Param("Low-pass Period", 10, 3, 20);
HPeriod = Param("High-pass Period", 48, 22, 80);
IsPlotHeatMap = ParamToggle("Show HeatMap?", "No|Yes", 1);
IsPlotDominantCycle = ParamToggle("Show Dom. Cycle?", "No|Yes");

pi=3.1415926;

function RoofingFilter(lpPeriod, hpPeriod)
{
	alpha1 = (cos(0.707*2*pi / hpPeriod) + sin(0.707*2*pi / hpPeriod) - 1) / cos(0.707*2*pi / hpPeriod);
	a1 = exp(-1.414*pi / lpPeriod);
	b1 = 2*a1*cos(1.414*pi / lpPeriod);
	c2 = b1;
	c3 = -a1*a1;
	c1 = 1 - c2 - c3;
	
	HP = Close;
	Filt = HP;
	
	for(i = 2; i < BarCount; i++)
	{
		HP[i] = ((1 - alpha1 / 2)^2)*(Close[i] - 2*Close[i-1] + Close[i-2]) + 2*(1 - alpha1)*HP[i-1] - ((1 - alpha1)^2)*HP[i-2];
		Filt[i] = c1*(HP[i] + HP[i-1]) / 2 + c2*Filt[i-1] + c3*Filt[i-2];
	}
	
	return Filt;
}

function AGC(lowerCutoff, higherCutoff, acceptableSlope)
{	
	factor = 0;
	accSlope = -acceptableSlope;	//acceptableSlope = 1.5 dB
	halfLC = lowerCutoff / 2;
	halfHC = higherCutoff / 2;
	ratio = 10^(accSlope/20);
	if(halfHC - halfLC > 0)
		factor = ratio^(1/(halfHC - halfLC));
	return factor;
}

function AutocorrelationPeriodogram(data, isHeatMap, isDomCyc)
{
	avgLength = 3;
	dominantCycle = 0;

	//Pearson correlation for each value of lag
	for(lag = 0; lag <= 48; lag++)
	{
		//Set the average length as M
		M = avgLength;
		if(avgLength == 0)
			M = lag;
		//Initialize correlation sums
		Sx = 0;
		Sy = 0;
		Sxx = 0;
		Syy = 0;
		Sxy = 0;
		//Advance samples of both data streams and sum Pearson components
		for(count = 0; count <= M-1; count++)
		{
			X = Ref(data, -count);
			Y = Ref(data, -(lag + count));
			Sx += X;
			Sy += Y;
			Sxx += X^2;
			Syy += Y^2;
			Sxy += X*Y;
		}
		var1 = (M*Sxx - Sx^2)*(M*Syy - Sy^2);
		VarSet("corr" + lag, IIf(var1 > 0, (M*Sxy - Sx*Sy)/sqrt(var1), 0));	//Compute correlation for each value of lag
		//VarSet("corrScale" + lag, IIf(var1 > 0, 0.5*((M*Sxy - Sx*Sy)/sqrt(var1) + 1), 0));	//Scale each correlation to range between 0 and 1
	}
	
	/*//Plot as a Heatmap (for scale each correlation to range between 0 and 1)
	for(period = 3; period <= 48; period++)
	{
		corr = VarGet("corrScale" + period);
		Red = IIf(corr > 0.5, 255*(2 - 2*corr), 255);
		Green = IIf(corr > 0.5, 255, 2*255*corr);
		PlotOHLC( period-1, period-1, period, period, "", ColorRGB( Red, Green, 0 ), styleCloud | styleNoLabel);
	}*/
	
	/*
		The DFT 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 the relative power at each period.
	*/
	for(period = 10; period <= 48; period++)
	{
		cosinePart = 0;
		sinePart = 0;
		
		for(n = 3; n <= 48; n++)
		{
			cosinePart += VarGet("corr" + n)*cos(2*pi*n / period);
			sinePart += VarGet("corr" + n)*sin(2*pi*n / period);
		}
		VarSet("sqSum" + period, cosinePart^2 + sinePart^2);
	}
	
	//EMA is used to smooth the power measurement at each period
	for(period = 10; period <= 48; period++)
		VarSet("r" + period, AMA((VarGet("sqSum" + period))^2, 0.2));

	//Find Maximum Power Level for Normalization
	K = AGC(10, 48, 1.5);
	for(period = 10; period <= 48; period++)
	{
		if(period == 10)
			VarSet("maxPwr", 0);		
		VarSet("maxPwr", IIf(VarGet("r" + period) > VarGet("maxPwr"), K*VarGet("r" + period), VarGet("maxPwr")));
	}
	
	//Normalization power
	for(period = 10; period <= 48; period++)
		VarSet("pwr" + period, VarGet("r" + period)/VarGet("maxPwr"));
		
	//Compute the dominant cycle using the CG of the spectrum
	Spx = 0;
	Sp = 0;
	for(period = 10; period <= 48; period++)
	{
		Spx += IIf(VarGet("pwr" + period) >= 0.5,  period*VarGet("pwr" + period), 0);
		Sp += IIf(VarGet("pwr" + period) >= 0.5, VarGet("pwr" + period), 0);
	}
	dominantCycle = IIf(Sp != 0, Spx / Sp, 0);
	
	if(isHeatMap)
	{
		//Plot as a Heatmap
		for(period = 10; period <= 48; period++)
		{
			pwr = VarGet("pwr" + period);
			Red = IIf(pwr > 0.5, 255, 2*255*pwr);
			Green = IIf(pwr > 0.5, 255*(2*pwr - 1), 0);
			PlotOHLC( period-1, period-1, period, period, "", ColorRGB( Red, Green, 0 ), styleCloud | styleNoLabel, Null, Null, 0, 0);
		}
	}
	
	if(isDomCyc)
		Plot(dominantCycle, "Dominant Cycle", colorBlue, styleThick, Null, Null, 0, 1);
	
	return dominantCycle;
}

filtData = RoofingFilter(LPeriod, HPeriod);
AutocorrelationPeriodogram(filtData, IsPlotHeatMap, IsPlotDominantCycle);