Adaptive VWAP Stdev Bands
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © simwai

//@version=5
indicator('VWAP Stdev Bands', overlay=true)

// -- Input --
src = input.source(hl2, 'Choose Source', group='General', inline='1')
resetMode = input.string('Adaptive – Homodyne Discriminator', 'Choose Reset Mode', options=['Resolution', 'Adaptive – Inphase-Quadrature Transform', 'Adaptive – Hilbert Transform', 'Adaptive – Homodyne Discriminator'], group='General', inline='2')
resetResolution = input.timeframe('D', 'Choose Reset Timeframe', group='General', inline='3', tooltip='Used for resolution reset mode')
minLength = input.int(2, minval=2, title='Min Length', group='General', inline='4', tooltip='Used for adaptive reset mode')
maxLength = input.int(200, minval=2, title='Max Length', group='General', inline='4', tooltip='Used for adaptive reset mode')
showPrevVWAP = input(false, title='Enable Previous VWAP Close Plot', group='General', inline='5')
isSmoothingEnabled  = input(false, 'Enable Smoothing', group='General', inline='6', tooltip='Hann Window is the smoothing method')
smoothingLength = input.int(10, 'Smoothing Length', minval=2, group='General', inline='6')

devUp1 = input(1.28, title='Stdev above (1)', group='Band Levels', inline='1')
devDn1 = input(1.28, title='Stdev below (1)', group='Band Levels', inline='1')

devUp2 = input(2.01, title='Stdev above (2)', group='Band Levels', inline='2')
devDn2 = input(2.01, title='Stdev below (2)', group='Band Levels', inline='2')

devUp3 = input(2.51, title='Stdev above (3)', group='Band Levels', inline='3')
devDn3 = input(2.51, title='Stdev below (3)', group='Band Levels', inline='3')

devUp4 = input(3.09, title='Stdev above (4)', group='Band Levels', inline='4')
devDn4 = input(3.09, title='Stdev below (4)', group='Band Levels', inline='4')

devUp5 = input(4.01, title='Stdev above (5)', group='Band Levels', inline='5')
devDn5 = input(4.01, title='Stdev below (5)', group='Band Levels', inline='5')

showDv1 = input(true, title='Show first group of bands?', group='Band Levels', inline='6')
showDv2 = input(true, title='Show second group of bands?', group='Band Levels', inline='7')
showDv3 = input(true, title='Show third group of bands?', group='Band Levels', inline='8')
showDv4 = input(true, title='Show fourth group of bands?', group='Band Levels', inline='9')
showDv5 = input(true, title='Show fifth group of bands?', group='Band Levels', inline='10')

// -- Colors --
color lightGray = color.rgb(214, 214, 214)

color lightCyan = color.rgb(207, 243, 235)
color aeroBlue = color.rgb(189, 239, 227)
color middleBlueGreen = color.rgb(156, 231, 214)
color middleBlueGreen2 = color.rgb(140, 227, 207)
color turquoise = color.rgb(107, 219, 193)

color mikadoYellow = color.rgb(255, 197, 51)
color selectiveYellow = color.rgb(255, 183, 3)
color orangePeel = color.rgb(253, 158, 2)
color yellowOrangeColorWheel = color.rgb(252, 146, 1)
color tangerine = color.rgb(251, 133, 0)

color transpTurquoise = color.new(turquoise, 50)
color transpTangerine = color.new(tangerine, 50)

// -- Functions --
// Get Dominant Cyber Cycle – Inphase-Quadrature -- Credits to @DasanC
getIqDc(float src, int min, int max) =>
    PI = 3.14159265359
    P = src - src[7]
    lenIQ = 0.0
    lenC = 0.0
    imult = 0.635
    qmult = 0.338
    inphase = 0.0
    quadrature = 0.0
    re = 0.0
    im = 0.0
    deltaIQ = 0.0
    instIQ = 0.0
    V = 0.0

    inphase := 1.25 * (P[4] - imult * P[2]) + imult * nz(inphase[3])
    quadrature := P[2] - qmult * P + qmult * nz(quadrature[2])
    re := 0.2 * (inphase * inphase[1] + quadrature * quadrature[1]) + 0.8 * nz(re[1])
    im := 0.2 * (inphase * quadrature[1] - inphase[1] * quadrature) + 0.8 * nz(im[1])
    if re != 0.0
        deltaIQ := math.atan(im / re)
        deltaIQ
    for i = 0 to max by 1
        V += deltaIQ[i]
        if V > 2 * PI and instIQ == 0.0
            instIQ := i
            instIQ
    if instIQ == 0.0
        instIQ := nz(instIQ[1])
        instIQ
    lenIQ := 0.25 * instIQ + 0.75 * nz(lenIQ[1], 1)
	lenIQ < min ? min : lenIQ

// Get Dominant Cyber Cycle – Hilbert Transform -- Credits to @DasanC
getHtDc(float src) =>
    Imult = .635
    Qmult = .338
    PI = 3.14159
    InPhase = 0.0
    Quadrature = 0.0
    Phase = 0.0
    DeltaPhase = 0.0
    InstPeriod = 0.0
    Period = 0.0
    Value4 = 0.0

    if bar_index > 5
        // Detrend src
        Value3 = src - src[7]
        // Compute InPhase and Quadrature components
        InPhase := 1.25 * (Value3[4] - Imult * Value3[2]) + Imult * nz(InPhase[3])
        Quadrature := Value3[2] - Qmult * Value3 + Qmult * nz(Quadrature[2])
        // Use ArcTangent to compute the current phase
        if math.abs(InPhase + InPhase[1]) > 0
            Phase := 180 / PI * math.atan(math.abs((Quadrature + Quadrature[1]) / (InPhase + InPhase[1])))
            Phase
        // Resolve the ArcTangent ambiguity
        if InPhase < 0 and Quadrature > 0
            Phase := 180 - Phase
            Phase
        if InPhase < 0 and Quadrature < 0
            Phase := 180 + Phase
            Phase
        if InPhase > 0 and Quadrature < 0
            Phase := 360 - Phase
            Phase
        // Compute a differential phase, resolve phase wraparound, and limit delta phase errors
        DeltaPhase := Phase[1] - Phase
        if Phase[1] < 90 and Phase > 270
            DeltaPhase := 360 + Phase[1] - Phase
            DeltaPhase
        if DeltaPhase < 1
            DeltaPhase := 1
            DeltaPhase
        if DeltaPhase > 60
            DeltaPhase := 60
            DeltaPhase
        // Sum DeltaPhases to reach 360 degrees. The sum is the instantaneous period.
        for i = 0 to 50 by 1
            Value4 += DeltaPhase[i]
            if Value4 > 360 and InstPeriod == 0
                InstPeriod := i
                InstPeriod
        // Resolve Instantaneous Period errors and smooth
        if InstPeriod == 0
			nz(InstPeriod[1])
        Period := .25 * InstPeriod + .75 * Period[1]
        Period
    Period

// Get Dominant Cyber Cycle - Homodyne Discriminator With Hilbert Dominant Cycle -- Credits to @blackcat1402
getHdDc(float Price, int min, int max, simple float CycPart = 0.5) =>
    Smooth = 0.00
    Detrender = 0.00
    I1 = 0.00
    Q1 = 0.00
    jI = 0.00
    jQ = 0.00
    I2 = 0.00
    Q2 = 0.00
    Re = 0.00
    Im = 0.00
    Period = 0.00
    SmoothPeriod = 0.00
    pi = 2 * math.asin(1)
    DomCycle = 0.0

    //Hilbert Transform
    Smooth := bar_index > 7 ? (4 * Price + 3 * nz(Price[1]) + 2 * nz(Price[2]) + nz(Price[3])) / 10 : Smooth
    Detrender := bar_index > 7 ? (.0962 * Smooth + .5769 * nz(Smooth[2]) - .5769 * nz(Smooth[4]) - .0962 * nz(Smooth[6])) * (.075 * nz(Period[1]) + .54) : Detrender
    //Compute InPhase and Quadrature components
    Q1 := bar_index > 7 ? (.0962 * Detrender + .5769 * nz(Detrender[2]) - .5769 * nz(Detrender[4]) - .0962 * nz(Detrender[6])) * (.075 * nz(Period[1]) + .54) : Q1
    I1 := bar_index > 7 ? nz(Detrender[3]) : I1

    //Advance the phase of I1 and Q1 by 90 degrees
    jI := (.0962 * I1 + .5769 * nz(I1[2]) - .5769 * nz(I1[4]) - .0962 * nz(I1[6])) * (.075 * nz(Period[1]) + .54)
    jQ := (.0962 * Q1 + .5769 * nz(Q1[2]) - .5769 * nz(Q1[4]) - .0962 * nz(Q1[6])) * (.075 * nz(Period[1]) + .54)

    //Phasor addition for 3 bar averaging
    I2 := I1 - jQ
    Q2 := Q1 + jI

    //Smooth the I and Q components before applying the discriminator
    I2 := .2 * I2 + .8 * nz(I2[1])
    Q2 := .2 * Q2 + .8 * nz(Q2[1])

    //Homodyne Discriminator
    Re := I2 * nz(I2[1]) + Q2 * nz(Q2[1])
    Im := I2 * nz(Q2[1]) - Q2 * nz(I2[1])
    Re := .2 * Re + .8 * nz(Re[1])
    Im := .2 * Im + .8 * nz(Im[1])

    Period := Im != 0 and Re != 0 ? 2 * pi / math.atan(Im / Re) : Period
    Period := Period > 1.5 * nz(Period[1]) ? 1.5 * nz(Period[1]) : Period
    Period := Period < .67 * nz(Period[1]) ? .67 * nz(Period[1]) : Period

    Period := Period < min ? min : Period
    Period := Period > max ? max : Period
    Period := .2 * Period + .8 * nz(Period[1])
    SmoothPeriod := .33 * Period + .67 * nz(SmoothPeriod[1])
    //it can add filter to Period here
    DomCycle := math.ceil(CycPart * SmoothPeriod) > 34 ? 34 : math.ceil(CycPart * SmoothPeriod) < 1 ? 1 : math.ceil(CycPart * SmoothPeriod)
    // Limit dominant cycle
    DomCycle := DomCycle < min ? min : DomCycle
    DomCycle := DomCycle > max ? max : DomCycle

// Zero-lag Least Squared Moving Average (ZLSMA) -- Credits to @veryfid
zlsma(float _src, int _period, int _offset=0) => 
    lsma = ta.linreg(_src, _period, _offset)
    lsma2 = ta.linreg(lsma, _period, _offset)
    eq = lsma - lsma2
    zlsma = lsma + eq

// -- Calculation --
float _src = src[1]
float _volume = volume[1]
float _close = close[1]
float _high = high[1] 
float _low = low[1]

float start = request.security(syminfo.tickerid, resetResolution, _src, barmerge.gaps_off, barmerge.lookahead_on)
int newSession = ta.change(start) ? 1 : 0

bool isReset = switch resetMode
    'Resolution' => newSession
    'Adaptive – Inphase-Quadrature Transform' => bar_index % math.round(getIqDc(_src, minLength, maxLength) / 2) == 0
    'Adaptive – Hilbert Transform' => bar_index % math.round(getHtDc(_src) / 2) == 0
    'Adaptive – Homodyne Discriminator' => bar_index % math.round(getHdDc(_src, minLength, maxLength)) == 0

vwapsum = 0.0
volumesum = 0.0
v2sum = 0.0
vwapsum := isReset ? _src * _volume : vwapsum[1] + _src * _volume
volumesum := isReset ? _volume : volumesum[1] + _volume
v2sum := isReset ? _volume * _src * _src : v2sum[1] + _volume * _src * _src
myvwap = vwapsum / volumesum
dev = math.sqrt(math.max(v2sum / volumesum - myvwap * myvwap, 0))

myvwap := isSmoothingEnabled ? zlsma(myvwap, smoothingLength) : myvwap

plot(myvwap, title='VWAP', color=lightGray)

getUpperPlot(float stdev, float vwap, float dev) => myvwap + stdev * dev
getLowerPlot(float stdev, float vwap, float dev) => myvwap - stdev * dev

upperPlot1 = plot(showDv1 ? getUpperPlot(devUp1, myvwap, dev) : na, title='VWAP Upper (1)', color=mikadoYellow)
lowerPlot1 = plot(showDv1 ? getLowerPlot(devDn1, myvwap, dev) : na, title='VWAP Lower (1)', color=lightCyan)

upperPlot2 = plot(showDv2 ? getUpperPlot(devUp2, myvwap, dev) : na, title='VWAP Upper (2)', color=selectiveYellow)
lowerPlot2 = plot(showDv2 ? getLowerPlot(devDn2, myvwap, dev) : na, title='VWAP Lower (2)', color=aeroBlue)

upperPlot3 = plot(showDv3 ? getUpperPlot(devUp3, myvwap, dev) : na, title='VWAP Upper (3)', color=orangePeel)
lowerPlot3 = plot(showDv3 ? getLowerPlot(devDn3, myvwap, dev) : na, title='VWAP Lower (3)', color=middleBlueGreen)

upperPlot4 = plot(showDv4 ? getUpperPlot(devUp4, myvwap, dev) : na, title='VWAP Upper (4)', color=yellowOrangeColorWheel)
lowerPlot4 = plot(showDv4 ? getLowerPlot(devDn4, myvwap, dev) : na, title='VWAP Lower (4)', color=middleBlueGreen2)

upperPlot5 = plot(showDv5 ? getUpperPlot(devUp5, myvwap, dev) : na, title='VWAP Upper (5)', color=tangerine)
lowerPlot5 = plot(showDv5 ? getLowerPlot(devDn5, myvwap, dev) : na, title='VWAP Lower (5)', color=turquoise)

fill(upperPlot5, upperPlot4, color=showDv5 ? transpTangerine : na)
fill(lowerPlot5, lowerPlot4, color=showDv5 ? transpTurquoise : na)

fill(upperPlot4, upperPlot3, color=((not showDv5) and showDv4) ? transpTangerine : na)
fill(lowerPlot4, lowerPlot3, color=((not showDv5) and showDv4) ? transpTurquoise : na)

fill(upperPlot3, upperPlot2, color=((not showDv5) and (not showDv4) and showDv3) ? transpTangerine : na)
fill(lowerPlot3, lowerPlot2, color=((not showDv5) and (not showDv4) and showDv3) ? transpTurquoise : na)

fill(upperPlot2, upperPlot1, color=((not showDv5) and (not showDv4) and (not showDv3) and showDv2) ? transpTangerine : na)
fill(lowerPlot2, lowerPlot1, color=((not showDv5) and (not showDv4) and (not showDv3) and showDv2) ? transpTurquoise : na)

prevwap = 0.0
prevwap := newSession ? myvwap[1] : prevwap[1]
plot(showPrevVWAP ? prevwap : na, style=plot.style_circles, color=_src > prevwap ? middleBlueGreen : orangePeel)