Band Pass Normalized Suite (BPNS)
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © peacefulLizard50262

//@version=5 
indicator("BPNS", explicit_plot_zorder = true)

bes(float source = close, float alpha = 0.7) =>
    var float smoothed = na
    smoothed := na(smoothed) ? source : alpha * source + (1 - alpha) * nz(smoothed[1])

sma(source, length) =>
    // Initialize the moving average to the first element of the source
    float average = source[0]
    // Initialize the count to 1
    int count = 1

    // Calculate the moving average using the Welford method
    for i = 1 to length
        count := count + 1
        delta = source[i] - average
        average := average + delta / count

    // Return the moving average
    average

ssmma(src, len, smooth)=>
    smma = 0.0
    avg = sma(src, len)
    smma := na(smma[1]) ? avg : bes((smma[1] * (len - 1) + src) / len, smooth)

smma(src, len)=>
    smma = 0.0
    avg = sma(src, len)
    smma := na(smma[1]) ? avg : (smma[1] * (len - 1) + src) / len

ema(src, len, smoothing, a) => //{
    alpha = a / (len + 1)
    sum   = 0.0
    sum  := na(sum[1]) ? src : alpha * bes(src, smoothing) + (1 - alpha) * nz(bes(sum[1], smoothing))

dema(src, length, smooth, a)=>
    e1 = ema(src, length, smooth, a)
    e2 = ema(e1, length, smooth, a)
    dema = 2 * e1 - e2

tema(src, length, smooth, a) =>
    e1 = ema(src, length, smooth, a)
    e2 = ema(e1, length, smooth, a)
    e3 = ema(e2, length, smooth, a)
    tema = 3 * (e1 - e2) + e3

wma(src, len) =>
    ta.wma(src, len)
        //// dont look for too long...

hma(src, len)=>
    wma(2 * wma(src, int(len/2)) - wma(src, len), int(math.sqrt(len/2)))

thma(src, len) =>
    wma1 = wma(src, len)
    wma2 = wma(src, int(len/2))
    wma3 = 2 * wma2 - wma1
    hma = wma(wma3, int(math.sqrt(len/2)))
    3 * (wma1 - wma2) + hma

ehma(src, len, smoothing, a)=>
    ema(2 * ema(src, int(len/2), smoothing, a) - ema(src, len, smoothing, a), int(math.sqrt(len/2)), smoothing, a)
 
gma(values, length) =>
    // Define the standard deviation of the Gaussian distribution
    stddev = length / 4
    // Generate an array of indices from 1 to length
    indices = ta.range(1, length)
    // Calculate the Gaussian weights for each value in the array
    weights = math.exp(-0.5 * (math.pow((indices - length), 2) / math.pow(stddev, 2)))
    // Calculate the weighted sum of the values
    sum = math.sum(values * weights, length)
    // Return the moving average
    sum / math.sum(weights, length)

rms(source, length)=>
    math.sqrt(math.sum(math.pow(source, 2), length)/length)

max(source, outlier_level, dev_lookback)=>
    var float max = na
    src = array.new<float>()
    stdev = math.abs((source - bes(source, 0.1))/ta.stdev(source, dev_lookback))
    array.push(src, stdev < outlier_level ? source : -1.7976931348623157e+308)
    max := math.max(nz(max[1]), array.get(src, 0))

min(source, outlier_level, dev_lookback) =>
    var float min = na
    src = array.new<float>()
    stdev = math.abs((source - bes(source, 0.1))/ta.stdev(source, dev_lookback))
    array.push(src, stdev < outlier_level ? source : 1.7976931348623157e+308)
    min := math.min(nz(min[1]), array.get(src, 0))

min_max(src, outlier_level, dev_lookback) =>
    (src - min(src, outlier_level, dev_lookback))/(max(src, outlier_level, dev_lookback) - min(src, outlier_level, dev_lookback)) * 100

source = input.source(close, "Source", group = "Source")
moving = input.string("BES", "MA Style", ["BES","SSMMA","SMMA","SMA","RMS","EMA","DEMA","TEMA","WMA","HMA","EHMA","THMA","GMA"], group = "Source")
blur   = input.float(1, "Global Smoothing", 0.005, 1, 0.005, group = "Source")
bound  = input.float(70, "Bounds", 50, 100, group = "Source")

lf_brown      = input.float(0.15, "BES Low Pass Filter Cutoff", 0.001, 1, 0.001, tooltip = "Adjusts BES", group = "Browns")
hf_brown      = input.float(0.01, "BES High Pass Filter Cutoff", 0.001, 1, 0.001, tooltip = "Adjusts BES", group = "Browns")

lf_len        = input.int(200, "Low Pass Filter Length", 1, tooltip = "Inputs for the MA's.", group = "MA")
hf_len        = input.int(12, "High Pass Filter Length", 1, tooltip = "Inputs for the MA's.", group = "MA")
lf_alpha      = input.float(2, "Low Pass Filter Alpha", 0.05, 10, 0.05, tooltip = "Inputs for the MA's.", group = "MA")
hf_alpha      = input.float(2, "High Pass Filter Alpha", 0.05, 10, 0.05, tooltip = "Inputs for the MA's.", group = "MA")
lf_smooth     = input.float(0.7, "Low Pass Filter Smoothing", 0.005, 1, 0.005, tooltip = "Inputs for the MA's.", group = "MA")
hf_smooth     = input.float(0.7, "High Pass Filter Smoothing", 0.005, 1, 0.005, tooltip = "Inputs for the MA's.", group = "MA")
outlier_level = input.float(2, "Outlier Level", 0.01, 100, 0.01, group = "Outlier Detection")
dev_lookback  = input.int(30, "Outlier Length", 2, group = "Outlier Detection")

enable   = input.bool(false, "Enable Global Mean", group = "Mean")
mean     = input.string("BES", "MA Style", ["BES","SSMMA","SMMA","SMA","RMS","EMA","DEMA","TEMA","WMA","HMA","EHMA","THMA","GMA"], group = "Mean")
mean_brn = input.float(0.01, "Global Smoothing Browns Alpha", group = "Mean")
mean_len = input.int(14, "Global Smoothing MA Length", group = "Mean")
mean_alp = input.float(2, "Low Pass Filter Alpha", 0.05, 10, 0.05, group = "Mean")
mean_smo = input.float(0.7, "Low Pass Filter Smoothing", 0.005, 1, 0.005), group = "Mean"

ma(src, len, smooth, alpha, brown, moving)=>
    switch moving
        "BES"   => bes(src, brown)
        "SSMMA" => ssmma(src, len, smooth)
        "SMMA"  => smma(src, len)
        "SMA"   => sma(src, len)
        "RMS"   => rms(src, len)
        "EMA"   => ema(src, len, smooth, alpha)
        "DEMA"  => dema(src, len, smooth, alpha)
        "TEMA"  => tema(src, len, smooth, alpha)
        "WMA"   => wma(src, len)
        "HMA"   => hma(src, len)
        "EHMA"  => ehma(src, len, smooth, alpha)
        "THMA"  => thma(src, len)
        "GMA"   => gma(src, len)
        "RMS"   => rms(src, len)

band_pass = bes(min_max(ma(source - ma(source, hf_len, hf_smooth, hf_alpha, hf_brown, moving), lf_len, lf_smooth, lf_alpha, lf_brown, moving), outlier_level, dev_lookback), blur)
g_mean    = ma(band_pass, mean_len, mean_smo, mean_alp, mean_brn, mean)

plot(band_pass > 100 ? 100 : band_pass < 0 ? 0 : band_pass, "Band Pass", color.blue)
plot(enable ? g_mean : na, "Global Mean", color.orange)

_100 = hline(bound, "Upper Bound", linestyle = hline.style_solid)
_50  = hline(50, "50", color.new(color.gray, 50))
_0   = hline(100 - bound, "Lower Bound", linestyle = hline.style_solid)

fill(_100, _0, color.new(color.navy, 90))