Ehlers Homodyne Discriminator (HOMOD)

// Licensed under the Apache License, Version 2.0
// © mihakralj
//@version=6
indicator("Ehlers Homodyne Discriminator (HOMOD)","HOMOD",overlay=false)

//@function Quadrant-aware angle calculation using stable atan2
//@param y Imaginary component
//@param x Real component
//@returns Angle in radians from -π to π
atan2(series float y,series float x)=>
    if y==0.0 and x==0.0
        runtime.error("atan2: y and x cannot both be zero")
    float ay=math.abs(y)
    float ax=math.abs(x)
    float angle=0.0
    if ax>ay
        angle:=math.atan(ay/ax)
    else
        angle:=(math.pi/2.0)-math.atan(ax/ay)
    if x<0.0
        angle:=math.pi-angle
    if y<0.0
        angle:=-angle
    angle

//@function Measures dominant cycle period using Ehlers homodyne discriminator
//@param source Price input series
//@param minPeriod Minimum dominant cycle length
//@param maxPeriod Maximum dominant cycle length
//@returns Smoothed dominant cycle estimate
//@optimized Exponential warmup compensation for dominant cycle smoothing
//@validation wolfram:"atan2(y,x)" external:"TradingView Homodyne Discriminator","Forex-Station Homodyne Discriminator","MQL5 Adaptive Lookback Homodyne","tindicators hd.cc"
homod(series float source,simple float minPeriod,simple float maxPeriod)=>
    if minPeriod<=0
        runtime.error("Min period must be greater than 0")
    if maxPeriod<=minPeriod
        runtime.error("Max period must be greater than min period")
    var float smooth_price=0.0
    var float detrender=0.0
    var float i1=0.0
    var float q1=0.0
    var float ji=0.0
    var float jq=0.0
    var float i2=0.0
    var float q2=0.0
    var float re=0.0
    var float im=0.0
    var float period=15.0
    var float smooth_period=15.0
    var float warm_decay=1.0
    var bool warmup=true
    float price=nz(source)
    float bandwidth=0.075*smooth_period+0.54
    smooth_price:=(4.0*price+3.0*nz(price[1])+2.0*nz(price[2])+nz(price[3]))/10.0
    detrender:=(0.0962*smooth_price+0.5769*nz(smooth_price[2])-0.5769*nz(smooth_price[4])-0.0962*nz(smooth_price[6]))*bandwidth
    q1:=(0.0962*detrender+0.5769*nz(detrender[2])-0.5769*nz(detrender[4])-0.0962*nz(detrender[6]))*bandwidth
    i1:=nz(detrender[3])
    ji:=(0.0962*i1+0.5769*nz(i1[2])-0.5769*nz(i1[4])-0.0962*nz(i1[6]))*bandwidth
    jq:=(0.0962*q1+0.5769*nz(q1[2])-0.5769*nz(q1[4])-0.0962*nz(q1[6]))*bandwidth
    float i2_raw=i1-jq
    float q2_raw=q1+ji
    i2:=0.2*i2_raw+0.8*nz(i2[1])
    q2:=0.2*q2_raw+0.8*nz(q2[1])
    float re_raw=i2*nz(i2[1])+q2*nz(q2[1])
    float im_raw=i2*nz(q2[1])-q2*nz(i2[1])
    re:=0.2*re_raw+0.8*nz(re[1])
    im:=0.2*im_raw+0.8*nz(im[1])
    float magnitude=math.abs(re)+math.abs(im)
    if magnitude>1e-10
        float angle=atan2(im,re)
        if math.abs(angle)>1e-10
            float candidate=2.0*math.pi/angle
            float clamped=math.max(minPeriod,math.min(maxPeriod,math.abs(candidate)))
            period:=0.2*clamped+0.8*period
    // Rate limiter: ±50% bar-to-bar, then clamp 6..50
    float prevPeriod = nz(period[1], 15.0)
    period := math.max(period, 0.67 * prevPeriod)
    period := math.min(period, 1.5 * prevPeriod)
    period := math.max(period, 6.0)
    period := math.min(period, 50.0)
    smooth_period := 0.2 * period + 0.8 * nz(smooth_period[1], period)
    float result=smooth_period
    float alpha=0.2
    if warmup
        warm_decay*=1.0-alpha
        float denom=1.0-warm_decay
        result:=denom>1e-10?result/denom:result
        warmup:=warm_decay>1e-10
    result

// ---------- Main loop ----------
i_source=input.source(hlc3,"Source")
i_minPeriod=input.float(6,"Min Period",minval=1,maxval=5000,step=0.5)
i_maxPeriod=input.float(50,"Max Period",minval=2,maxval=5000,step=0.5)
homodPeriod=homod(i_source,i_minPeriod,i_maxPeriod)
plot(homodPeriod,"Dominant Cycle Period",color=color.yellow,linewidth=2)