"BB-MACD" Created by KelvinHand Using Ehlers Periodogram for the Dominant Cycle

/*
        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         {
                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         {
                //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                 {
                        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         {
                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         {
                cosinePart = 0;
                sinePart = 0;
               
                for(n = 3; 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                 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         {
                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                 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         {
                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);
        dominantCycle =(dominantCycle);
       
        if(isHeatMap)
        {
                //Plot as a Heatmap
                for(period = 10; 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(LastValue(dominantCycle,1), "Dominant Cycle", colorBlue, styleThick, Null, Null, 0, 1);

        pds=LastValue(dominantCycle,1);

//"BB-MACD")    Created by KelvinHand
       
SetChartBkColor( ParamColor("background",colorBlack) );
A1=EMA(C,pds/2)-EMA(C,pds);
BBtop=BBandTop(A1,10,1);
BBbot=BBandBot(A1,10,1);
Color=IIf(a1<0 AND a1>Ref(a1,-1), colorLime,IIf(a1>0 AND a1>Ref(a1,-1),colorBrightGreen,IIf(a1>0 AND a1
Plot(a1,"MACD",color,styleThick+styleLine);
Plot(BBtop,"BBtop",colorWhite,styleDashed);
Plot(BBbot,"BBbot",colorWhite,styleDashed);
Plot(0,"",colorWhite,1);

        AddColumn(round(dominantCycle), " DominantCycle", 2.0 ,  colorRed, colorLightYellow,50) ;

        return dominantCycle;

}