diff --git a/_sidebar.md b/_sidebar.md index dbfd22d4..9bbb7486 100644 --- a/_sidebar.md +++ b/_sidebar.md @@ -329,6 +329,7 @@ * [HT_DCPHASE - Ehlers Hilbert Transform Dominant Cycle Phase](/lib/cycles/ht_dcphase/HtDcphase.md) * [HT_PHASOR - Ehlers Hilbert Transform Phasor Components](/lib/cycles/ht_phasor/HtPhasor.md) * [HT_SINE - Ehlers Hilbert Transform SineWave (also known as SINE)](/lib/cycles/ht_sine/HtSine.md) + * [LPF - Ehlers Linear Predictive Filter](/lib/cycles/lpf/Lpf.md) * [LUNAR - Lunar Phase](/lib/cycles/lunar/Lunar.md) * [SOLAR - Solar Activity Cycle](/lib/cycles/solar/Solar.md) * [SSFDSP - Ehlers SSF Detrended Synthetic Price](/lib/cycles/ssfdsp/Ssfdsp.md) diff --git a/docs/indicators.md b/docs/indicators.md index 8605bb9a..19c63a36 100644 --- a/docs/indicators.md +++ b/docs/indicators.md @@ -447,6 +447,7 @@ Periodic pattern detection and dominant frequency extraction. Markets exhibit cy | [**HT_DCPHASE**](../lib/cycles/ht_dcphase/HtDcphase.md) | Ehlers HT Dominant Cycle Phase | Hilbert Transform phase angle | | [**HT_PHASOR**](../lib/cycles/ht_phasor/HtPhasor.md) | Ehlers HT Phasor Components | In-phase and quadrature components | | [**HT_SINE**](../lib/cycles/ht_sine/HtSine.md) | Ehlers HT SineWave (also known as SINE) | Dominant cycle phase with lead signal | +| [**LPF**](../lib/cycles/lpf/Lpf.md) | Ehlers Linear Predictive Filter | Griffiths LMS predictor ? DFT spectrum ? dominant cycle | | [**LUNAR**](../lib/cycles/lunar/Lunar.md) | Lunar Phase | Moon phase cycle | | [**SOLAR**](../lib/cycles/solar/Solar.md) | Solar Activity Cycle | Solar activity periodicity | | [**SSFDSP**](../lib/cycles/ssfdsp/Ssfdsp.md) | Ehlers SSF Detrended Synthetic Price | Dual Super Smoother oscillator | diff --git a/docs/pinescript.md b/docs/pinescript.md index 3502e9e7..2dee8743 100644 --- a/docs/pinescript.md +++ b/docs/pinescript.md @@ -435,6 +435,7 @@ Markets oscillate. These indicators try to measure the oscillation itself β€” th | HT_DCPHASE | Ehlers Hilbert Transform Dominant Cycle Phase | [ht_dcphase.pine](../lib/cycles/ht_dcphase/ht_dcphase.pine) | | HT_PHASOR | Ehlers Hilbert Transform Phasor Components | [phasor.pine](../lib/cycles/ht_phasor/phasor.pine) | | HT_SINE | Ehlers Hilbert Transform SineWave | [ht_sine.pine](../lib/cycles/ht_sine/ht_sine.pine) | +| LPF | Ehlers Linear Predictive Filter | [lpf.pine](../lib/cycles/lpf/lpf.pine) | | LUNAR | Lunar Phase | [lunar.pine](../lib/cycles/lunar/lunar.pine) | | SOLAR | Solar Activity Cycle | [solar.pine](../lib/cycles/solar/solar.pine) | | SSFDSP | Ehlers SSF Detrended Synthetic Price | [ssfdsp.pine](../lib/cycles/ssfdsp/ssfdsp.pine) | diff --git a/lib/_index.md b/lib/_index.md index ea297956..2e9c8c4c 100644 --- a/lib/_index.md +++ b/lib/_index.md @@ -1,4 +1,4 @@ -ο»Ώ# QuanTAlib Indicators +# QuanTAlib Indicators | Indicator | Full Name | Category | | :--------------------------------------------------------- | :------------------------------------------------------------------------ | :----------- | @@ -201,6 +201,7 @@ | [LOGNORMDIST](numerics/lognormdist/Lognormdist.md) | Log-Normal Distribution | Numerics | | [LOGTRANS](numerics/logtrans/Logtrans.md) | Logarithmic Transform | Numerics | | [LOWEST](numerics/lowest/Lowest.md) | Rolling Minimum | Numerics | +| [LPF](cycles/lpf/Lpf.md) | Ehlers Linear Predictive Filter | Cycles | | [LRSI](oscillators/lrsi/Lrsi.md) | Ehlers Laguerre RSI | Oscillators | | [LSMA](trends_FIR/lsma/Lsma.md) | Least Squares MA | Trends (FIR) | | [LTMA](trends_IIR/ltma/Ltma.md) | Linear Trend MA | Trends (IIR) | @@ -312,7 +313,7 @@ | [RSE](errors/rse/Rse.md) | Relative Squared Error | Errors | | [RSI](momentum/rsi/Rsi.md) | Relative Strength Index | Momentum | | [RSIH](oscillators/rsih/Rsih.md) | Ehlers Hann-Windowed RSI | Oscillators | -| [RSQUARED](errors/rsquared/Rsquared.md) | RΒ² (Coefficient of Determination) | Errors | +| [RSQUARED](errors/rsquared/Rsquared.md) | R² (Coefficient of Determination) | Errors | | [RSV](volatility/rsv/Rsv.md) | Rogers-Satchell Volatility | Volatility | | [RSX](momentum/rsx/Rsx.md) | Relative Strength Quality Index | Momentum | | [RV](volatility/rv/Rv.md) | Realized Volatility | Volatility | diff --git a/lib/cycles/_index.md b/lib/cycles/_index.md index 5cdfe326..f035dc1a 100644 --- a/lib/cycles/_index.md +++ b/lib/cycles/_index.md @@ -18,6 +18,7 @@ Cycle analysis identifies repeating patterns in price data. John Ehlers pioneere | [HT_DCPHASE](ht_dcphase/Htdcphase.md) | Ehlers Hilbert Transform Dominant Cycle Phase | Ehlers Hilbert Transform. Measures current position in cycle. | | [HT_PHASOR](ht_phasor/HtPhasor.md) | Ehlers Hilbert Transform Phasor Components | Ehlers. In-phase and quadrature components. | | [HT_SINE](ht_sine/HtSine.md) | Ehlers Hilbert Transform SineWave (also known as SINE) | Ehlers Hilbert Transform. Sine and lead sine for cycle timing. | +| [LPF](lpf/Lpf.md) | Ehlers Linear Predictive Filter | Ehlers. Griffiths LMS predictor coefficients β†’ DFT spectrum β†’ dominant cycle.| | [LUNAR](lunar/Lunar.md) | Lunar Phase | 29.5-day lunar cycle. Studied for market correlations. | | [SOLAR](solar/Solar.md) | Solar Activity Cycle | ~11-year sunspot cycle. Long-term research indicator. | | [SSFDSP](ssfdsp/Ssfdsp.md) | Ehlers SSF Detrended Synthetic Price | Super Smoother Filter based DSP. Cleaner cycle extraction. | diff --git a/lib/cycles/lpf/Lpf.Quantower.cs b/lib/cycles/lpf/Lpf.Quantower.cs new file mode 100644 index 00000000..bc4ec916 --- /dev/null +++ b/lib/cycles/lpf/Lpf.Quantower.cs @@ -0,0 +1,79 @@ +using System.Drawing; +using System.Runtime.CompilerServices; +using TradingPlatform.BusinessLayer; + +namespace QuanTAlib; + +[SkipLocalsInit] +public sealed class LpfIndicator : Indicator, IWatchlistIndicator +{ + [InputParameter("Lower Bound", sortIndex: 1, 8, 200, 1, 0)] + public int LowerBound { get; set; } = 18; + + [InputParameter("Upper Bound", sortIndex: 2, 10, 500, 1, 0)] + public int UpperBound { get; set; } = 40; + + [InputParameter("Data Length", sortIndex: 3, 4, 200, 1, 0)] + public int DataLength { get; set; } = 40; + + [IndicatorExtensions.DataSourceInput] + public SourceType Source { get; set; } = SourceType.Close; + + [InputParameter("Show cold values", sortIndex: 21)] + public bool ShowColdValues { get; set; } = true; + + private Lpf _lpf = null!; + private readonly LineSeries _cycleSeries; + private readonly LineSeries _signalSeries; + private readonly LineSeries _predictSeries; + private Func _priceSelector = null!; + + public static int MinHistoryDepths => 0; + int IWatchlistIndicator.MinHistoryDepths => MinHistoryDepths; + + public override string ShortName => $"LPF ({LowerBound},{UpperBound},{DataLength})"; + public override string SourceCodeLink => "https://github.com/mihakralj/QuanTAlib/blob/main/lib/cycles/lpf/Lpf.Quantower.cs"; + + public LpfIndicator() + { + OnBackGround = true; + SeparateWindow = true; + Name = "LPF - Ehlers Linear Predictive Filter"; + Description = "Ehlers' Linear Predictive Filter estimates the dominant cycle period using Griffiths adaptive coefficients and spectral analysis"; + + _cycleSeries = new LineSeries(name: "Cycle", color: IndicatorExtensions.Oscillators, width: 2, style: LineStyle.Solid); + _signalSeries = new LineSeries(name: "Signal", color: Color.Lime, width: 1, style: LineStyle.Solid); + _predictSeries = new LineSeries(name: "Predict", color: Color.Red, width: 1, style: LineStyle.Dot); + AddLineSeries(_cycleSeries); + AddLineSeries(_signalSeries); + AddLineSeries(_predictSeries); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + protected override void OnInit() + { + _lpf = new Lpf(LowerBound, UpperBound, DataLength); + _priceSelector = Source.GetPriceSelector(); + base.OnInit(); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + protected override void OnUpdate(UpdateArgs args) + { + if (args.Reason != UpdateReason.NewBar && args.Reason != UpdateReason.HistoricalBar) + { + return; + } + + var item = this.HistoricalData[this.Count - 1, SeekOriginHistory.Begin]; + double value = _priceSelector(item); + var time = this.HistoricalData.Time(); + + var input = new TValue(time, value); + TValue result = _lpf.Update(input, args.IsNewBar()); + + _cycleSeries.SetValue(result.Value, _lpf.IsHot, ShowColdValues); + _signalSeries.SetValue(_lpf.Signal * UpperBound * 0.5 + (LowerBound + UpperBound) * 0.5, _lpf.IsHot, ShowColdValues); + _predictSeries.SetValue(_lpf.Predict * UpperBound * 0.5 + (LowerBound + UpperBound) * 0.5, _lpf.IsHot, ShowColdValues); + } +} diff --git a/lib/cycles/lpf/Lpf.cs b/lib/cycles/lpf/Lpf.cs new file mode 100644 index 00000000..c6aa5ea7 --- /dev/null +++ b/lib/cycles/lpf/Lpf.cs @@ -0,0 +1,441 @@ +using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; + +namespace QuanTAlib; + +/// +/// LPF: Ehlers Linear Predictive Filter +/// Estimates the dominant cycle period using a Griffiths linear predictive filter +/// combined with spectral analysis of the adapted filter coefficients. +/// +/// +/// The algorithm is based on John F. Ehlers' "Linear Predictive Filters And +/// Instantaneous Frequency" from Technical Analysis of Stocks & Commodities, +/// January 2025. +/// +/// Algorithm stages: +/// 1. Roofing filter (HP + SuperSmoother) band-limits the input +/// 2. AGC normalizes the band-limited signal to Β±1 +/// 3. Griffiths LMS predictor adapts coefficients to minimize prediction error +/// 4. DFT of coefficients yields power spectrum +/// 5. Center of gravity of spectrum identifies dominant cycle +/// +/// Key properties: +/// - Self-calibrating via adaptive coefficient updates +/// - Dominant cycle constrained to change by at most 2 bars per update +/// - Also provides predictive signal and trigger outputs +/// - Based on Lloyd Griffiths' "Rapid Measurement of Digital Instantaneous +/// Frequency" (IEEE Trans. ASSP-23, 1975) +/// +/// Pine Script reference: +/// https://github.com/mihakralj/pinescript/blob/main/indicators/cycles/lpf.md +/// +/// Complexity: O(LengthΒ² + Length Γ— (UpperBound βˆ’ LowerBound)) per bar +/// +[SkipLocalsInit] +public sealed class Lpf : AbstractBase +{ + private readonly int _lowerBound; + private readonly int _upperBound; + private readonly int _dataLength; + + // Roofing filter coefficients (precomputed) + private readonly double _hpAlpha; + private readonly double _ssC1, _ssC2, _ssC3; + + // Griffiths coefficient arrays (need snapshot/restore) + private readonly double[] _coef; + private readonly double[] _p_coef; + + // Signal history buffer + private readonly double[] _xx; + private readonly double[] _p_xx; + + // Power spectrum buffer (reused each bar) + private readonly double[] _pwr; + + // Filter state + [StructLayout(LayoutKind.Auto)] + private record struct State( + double Src1, double Src2, + double Hp1, double Hp2, + double Lp1, double Lp2, + double Peak, + double Signal, + double Dom, double PrevDom, + double Predict, + double DomPower, + int BarCount, + double LastValid + ) + { + public static State New() => new() + { + Peak = 0.1, + Dom = 0.0, + PrevDom = 0.0, + LastValid = double.NaN + }; + } + + private State _state; + private State _p_state; + + private ITValuePublisher? _publisher; + private TValuePublishedHandler? _handler; + private bool _isNew; + + /// Gets the current dominant cycle period. + public double DominantCycle => _state.Dom; + + /// Gets the AGC-normalized band-limited signal (Β±1 range). + public double Signal => _state.Signal; + + /// Gets the predicted value (2-bar ahead prediction of the normalized signal). + public double Predict => _state.Predict; + + /// Lower bandpass boundary (minimum period). + public int LowerBound => _lowerBound; + + /// Upper bandpass boundary (maximum period). + public int UpperBound => _upperBound; + + /// Data length for Griffiths predictor. + public int DataLength => _dataLength; + + public bool IsNew => _isNew; + public override bool IsHot => _state.BarCount >= WarmupPeriod; + + /// + /// Creates a new Ehlers Linear Predictive Filter indicator. + /// + /// Lower bandpass boundary - minimum period to detect (β‰₯ 8). + /// Upper bandpass boundary - maximum period to detect (> lowerBound). + /// Data length for adaptive predictor (β‰₯ 4). Ehlers recommends matching upperBound. + public Lpf(int lowerBound = 18, int upperBound = 40, int dataLength = 40) + { + if (lowerBound < 8) + { + throw new ArgumentOutOfRangeException(nameof(lowerBound), "Lower bound must be at least 8."); + } + if (upperBound <= lowerBound) + { + throw new ArgumentOutOfRangeException(nameof(upperBound), "Upper bound must be greater than lower bound."); + } + if (dataLength < 4) + { + throw new ArgumentOutOfRangeException(nameof(dataLength), "Data length must be at least 4."); + } + + _lowerBound = lowerBound; + _upperBound = upperBound; + _dataLength = dataLength; + + Name = $"LPF({lowerBound},{upperBound},{dataLength})"; + WarmupPeriod = upperBound * 2; + + // Precompute HP Butterworth coefficient + double angle707 = 0.707 * 2.0 * Math.PI / upperBound; + _hpAlpha = (Math.Cos(angle707) + Math.Sin(angle707) - 1.0) / Math.Cos(angle707); + + // Precompute SuperSmoother (LP) coefficients + double sqrt2Pi = Math.Sqrt(2.0) * Math.PI / lowerBound; + double a1 = Math.Exp(-sqrt2Pi); + double b1 = 2.0 * a1 * Math.Cos(sqrt2Pi); + _ssC2 = b1; + _ssC3 = -(a1 * a1); + _ssC1 = 1.0 - _ssC2 - _ssC3; + + // Allocate arrays + _coef = new double[dataLength + 1]; + _p_coef = new double[dataLength + 1]; + _xx = new double[dataLength + 1]; + _p_xx = new double[dataLength + 1]; + _pwr = new double[upperBound + 2]; + + // Initialize state + _state = State.New(); + _state = _state with { Dom = (lowerBound + upperBound) * 0.5, PrevDom = (lowerBound + upperBound) * 0.5 }; + _p_state = _state; + } + + /// + /// Creates a chained LPF indicator subscribed to a source publisher. + /// + public Lpf(ITValuePublisher source, int lowerBound = 18, int upperBound = 40, int dataLength = 40) + : this(lowerBound, upperBound, dataLength) + { + ArgumentNullException.ThrowIfNull(source); + _publisher = source; + _handler = Handle; + source.Pub += _handler; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private void Handle(object? sender, in TValueEventArgs args) + { + Update(args.Value, args.IsNew); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public override TValue Update(TValue input, bool isNew = true) + { + _isNew = isNew; + if (isNew) + { + _p_state = _state; + Array.Copy(_coef, _p_coef, _coef.Length); + Array.Copy(_xx, _p_xx, _xx.Length); + } + else + { + _state = _p_state; + Array.Copy(_p_coef, _coef, _coef.Length); + Array.Copy(_p_xx, _xx, _xx.Length); + } + + // Handle non-finite input + double val = input.Value; + if (!double.IsFinite(val)) + { + val = double.IsFinite(_state.LastValid) ? _state.LastValid : 0.0; + } + else + { + _state = _state with { LastValid = val }; + } + + int barCount = isNew ? _state.BarCount + 1 : _state.BarCount; + + // ── Stage 1: Roofing Filter ────────────────────────────────────── + // Highpass: 2nd-order Butterworth + double hpCoef = (1.0 - _hpAlpha * 0.5); + double hpInput = hpCoef * hpCoef * (val - 2.0 * _state.Src1 + _state.Src2); + double oneMinusAlpha = 1.0 - _hpAlpha; + double hp = hpInput + 2.0 * oneMinusAlpha * _state.Hp1 - oneMinusAlpha * oneMinusAlpha * _state.Hp2; + + // SuperSmoother lowpass + double lp = Math.FusedMultiplyAdd(_ssC1, (hp + _state.Hp1) * 0.5, + Math.FusedMultiplyAdd(_ssC2, _state.Lp1, _ssC3 * _state.Lp2)); + + // ── Stage 2: AGC Normalization ─────────────────────────────────── + double peak = 0.991 * _state.Peak; + double absLp = Math.Abs(lp); + if (absLp > peak) + { + peak = absLp; + } + double signal = peak > 0.0 ? lp / peak : 0.0; + + // ── Stage 3: Griffiths Adaptive Predictor ──────────────────────── + // Shift data buffer (shift right by 1) + for (int i = _dataLength; i >= 1; i--) + { + _xx[i] = _xx[i - 1]; + } + _xx[0] = signal; + + // Compute signal power for normalization + double sigPower = 0.0; + for (int i = 0; i < _dataLength; i++) + { + sigPower += _xx[i] * _xx[i]; + } + sigPower /= _dataLength; + + // Convergence factor (ΞΌ) + double mu = sigPower > 0.0 ? 0.25 / (sigPower * _dataLength) : 0.0; + + // Predict current value from past + double xBar = 0.0; + for (int i = 1; i <= _dataLength; i++) + { + xBar += _coef[i] * _xx[i]; + } + + // Error = actual - predicted + double error = _xx[0] - xBar; + + // Update coefficients via LMS rule + for (int i = 1; i <= _dataLength; i++) + { + _coef[i] += mu * error * _xx[i]; + } + + // ── Stage 4: Spectrum from Coefficients (DFT) ──────────────────── + double maxPwrLocal = 0.0; + for (int period = _lowerBound; period <= _upperBound; period++) + { + double realPart = 0.0; + double imagPart = 0.0; + double angleStep = 2.0 * Math.PI / period; + for (int i = 1; i <= _dataLength; i++) + { + double angle = angleStep * i; + realPart += _coef[i] * Math.Cos(angle); + imagPart += _coef[i] * Math.Sin(angle); + } + double p = realPart * realPart + imagPart * imagPart; + _pwr[period] = p; + if (p > maxPwrLocal) + { + maxPwrLocal = p; + } + } + + // ── Stage 5: Dominant Cycle via Center of Gravity ──────────────── + double spx = 0.0; + double sp = 0.0; + for (int period = _lowerBound; period <= _upperBound; period++) + { + double normPwr = maxPwrLocal > 0.0 ? _pwr[period] / maxPwrLocal : 0.0; + if (normPwr >= 0.5) + { + spx += period * normPwr; + sp += normPwr; + } + } + + double rawDom = sp > 0.0 ? spx / sp : _state.Dom; + + // Constrain change to Β±2 bars per update (prevents bouncing) + double prevDom = _state.PrevDom; + if (rawDom - prevDom > 2.0) + { + rawDom = prevDom + 2.0; + } + if (prevDom - rawDom > 2.0) + { + rawDom = prevDom - 2.0; + } + + double dom = Math.Clamp(rawDom, _lowerBound, _upperBound); + + // Two-bar prediction for trigger line + double xPred = 0.0; + if (_dataLength > 2) + { + for (int i = 1; i <= _dataLength - 2; i++) + { + xPred += _coef[i] * _xx[i]; + } + } + + // Compute power at dominant cycle + int domIdx = Math.Clamp((int)Math.Round(dom), _lowerBound, _upperBound); + double domPower = maxPwrLocal > 0.0 ? Math.Clamp(_pwr[domIdx] / maxPwrLocal, 0.0, 1.0) : 0.0; + + // Update state + _state = new State( + Src1: val, Src2: _state.Src1, + Hp1: hp, Hp2: _state.Hp1, + Lp1: lp, Lp2: _state.Lp1, + Peak: peak, + Signal: signal, + Dom: dom, PrevDom: dom, + Predict: xPred, + DomPower: domPower, + BarCount: barCount, + LastValid: _state.LastValid + ); + + Last = new TValue(input.Time, dom); + PubEvent(Last, isNew); + return Last; + } + + public override TSeries Update(TSeries source) + { + if (source.Count == 0) + { + return []; + } + + int len = source.Count; + var t = new List(len); + var v = new List(len); + CollectionsMarshal.SetCount(t, len); + CollectionsMarshal.SetCount(v, len); + + var tSpan = CollectionsMarshal.AsSpan(t); + var vSpan = CollectionsMarshal.AsSpan(v); + + for (int i = 0; i < len; i++) + { + var result = Update(source[i]); + vSpan[i] = result.Value; + } + source.Times.CopyTo(tSpan); + + return new TSeries(t, v); + } + + public override void Prime(ReadOnlySpan source, TimeSpan? step = null) + { + foreach (double value in source) + { + Update(new TValue(DateTime.MinValue, value)); + } + } + + /// Computes LPF for a given time series. + public static TSeries Batch(TSeries source, int lowerBound = 18, int upperBound = 40, int dataLength = 40) + { + var indicator = new Lpf(lowerBound, upperBound, dataLength); + return indicator.Update(source); + } + + /// Computes LPF in-place using pre-allocated output span. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static void Batch(ReadOnlySpan source, Span output, + int lowerBound = 18, int upperBound = 40, int dataLength = 40) + { + if (source.Length != output.Length) + { + throw new ArgumentException("Source and output must have the same length.", nameof(output)); + } + + var indicator = new Lpf(lowerBound, upperBound, dataLength); + for (int i = 0; i < source.Length; i++) + { + var result = indicator.Update(new TValue(DateTime.MinValue, source[i])); + output[i] = result.Value; + } + } + + /// Calculates LPF for a time series and returns both results and indicator state. + public static (TSeries Results, Lpf Indicator) Calculate(TSeries source, + int lowerBound = 18, int upperBound = 40, int dataLength = 40) + { + var indicator = new Lpf(lowerBound, upperBound, dataLength); + TSeries results = indicator.Update(source); + return (results, indicator); + } + + public override void Reset() + { + _state = State.New(); + _state = _state with { Dom = (_lowerBound + _upperBound) * 0.5, PrevDom = (_lowerBound + _upperBound) * 0.5 }; + _p_state = _state; + Array.Clear(_coef); + Array.Clear(_p_coef); + Array.Clear(_xx); + Array.Clear(_p_xx); + Array.Clear(_pwr); + Last = default; + } + + /// + /// Unsubscribes from the source publisher if one was provided during construction. + /// + protected override void Dispose(bool disposing) + { + if (disposing && _publisher != null && _handler != null) + { + _publisher.Pub -= _handler; + _publisher = null; + _handler = null; + } + base.Dispose(disposing); + } +} diff --git a/lib/cycles/lpf/Lpf.md b/lib/cycles/lpf/Lpf.md new file mode 100644 index 00000000..3ad12cd8 --- /dev/null +++ b/lib/cycles/lpf/Lpf.md @@ -0,0 +1,134 @@ +# LPF: Ehlers Linear Predictive Filter + +Griffiths-adapted dominant cycle estimator β€” uses LMS-predicted filter coefficients as a spectral window. + +| Property | Value | +|:-------------- |:---------------------------------------------------- | +| **Category** | Cycles | +| **Inputs** | Single series (close) | +| **Parameters** | `lowerBound` (int, 18), `upperBound` (int, 40), `dataLength` (int, 40) | +| **Outputs** | Dominant Cycle (period), Signal (Β±1 AGC), Predict | +| **Output range**| \[lowerBound, upperBound\] | +| **Warmup** | `2 Γ— upperBound` | +| **PineScript** | [lpf.pine](lpf.pine) | + +- Applies a roofing filter (HP + SuperSmoother) to produce band-limited data, then adapts Griffiths LMS coefficients to minimize one-step prediction error. +- Transforms the adapted coefficients into a frequency-domain power spectrum via DFT, identifying the dominant cycle as the spectral center of gravity. +- Constrains the dominant cycle to change by at most 2 bars per update, preventing erratic mode-switching in noisy data. + +## Historical Context + +John F. Ehlers introduced the Linear Predictive Filter in "Linear Predictive Filters And Instantaneous Frequency" (*Technical Analysis of Stocks & Commodities*, January 2025). The algorithm adapts Lloyd Griffiths' "Rapid Measurement of Digital Instantaneous Frequency" (IEEE Trans. ASSP-23, 1975) β€” a time-domain gradient method for adaptive spectral estimation originally developed for radar and sonar signal processing. Ehlers' innovation was combining this with his roofing filter and AGC normalization to create a self-calibrating cycle detector for financial data. Unlike his earlier Autocorrelation Periodogram (ACP), LPF estimates the spectrum from the *predictor coefficients* rather than from autocorrelation lags, yielding sharper spectral resolution with fewer data points. + +## Mathematical Foundation + +### Stage 1: Roofing Filter (Band-Limiting) + +A 2nd-order Butterworth highpass removes trend (periods > `upperBound`), followed by a SuperSmoother lowpass that removes noise (periods < `lowerBound`): + +$$\alpha_{HP} = \frac{\cos(0.707 \cdot 2\pi / U) + \sin(0.707 \cdot 2\pi / U) - 1}{\cos(0.707 \cdot 2\pi / U)}$$ + +$$HP_n = (1-\tfrac{\alpha}{2})^2 (x_n - 2x_{n-1} + x_{n-2}) + 2(1-\alpha)\,HP_{n-1} - (1-\alpha)^2\,HP_{n-2}$$ + +$$a_1 = e^{-\sqrt{2}\pi/L}, \quad b_1 = 2a_1\cos(\sqrt{2}\pi/L)$$ + +$$LP_n = (1-b_1+a_1^2)\tfrac{HP_n+HP_{n-1}}{2} + b_1\,LP_{n-1} - a_1^2\,LP_{n-2}$$ + +### Stage 2: AGC Normalization + +$$\text{Peak}_n = \max(0.991 \cdot \text{Peak}_{n-1},\; |LP_n|)$$ + +$$\text{Signal}_n = LP_n / \text{Peak}_n$$ + +### Stage 3: Griffiths LMS Predictor + +The heart of the algorithm β€” adaptive coefficient update minimizing prediction error: + +$$P_{\text{sig}} = \frac{1}{N}\sum_{i=0}^{N-1} x_i^2, \qquad \mu = \frac{0.25}{P_{\text{sig}} \cdot N}$$ + +$$\hat{x}_0 = \sum_{i=1}^{N} c_i \cdot x_i, \qquad \varepsilon = x_0 - \hat{x}_0$$ + +$$c_i \leftarrow c_i + \mu \cdot \varepsilon \cdot x_i \quad \forall\, i \in [1, N]$$ + +### Stage 4: Spectral Estimation via DFT of Coefficients + +$$\text{Pwr}(P) = \left(\sum_{i=1}^{N} c_i \cos\tfrac{2\pi i}{P}\right)^2 + \left(\sum_{i=1}^{N} c_i \sin\tfrac{2\pi i}{P}\right)^2$$ + +### Stage 5: Dominant Cycle (Center of Gravity) + +$$DC = \frac{\sum_{P: \text{Pwr}(P) \geq 0.5} P \cdot \text{Pwr}(P)}{\sum_{P: \text{Pwr}(P) \geq 0.5} \text{Pwr}(P)}, \qquad |\Delta DC| \leq 2$$ + +### Parameter Mapping + +| Parameter | Effect | Recommended | +|:------------- |:----------------------------------- |:--------------------- | +| `lowerBound` | Shortest cycle detected | 18 swing, 8 minimum | +| `upperBound` | Longest cycle detected | 40 swing, 125+ position | +| `dataLength` | Predictor adaptation window | Match `upperBound` | + +## Performance Profile + +### Operation Count (Streaming Mode, Scalar) + +| Operation | Count per bar | +|:---------------- |:------------------------------------ | +| HP filter | 5 mul, 4 add | +| SuperSmoother | 3 mul, 3 add | +| AGC | 2 mul, 1 cmp | +| Buffer shift | N copies | +| Signal power | N mul, N add | +| LMS predict | N mul, N add | +| Coef update | 2N mul, N add | +| DFT spectrum | 2N(Uβˆ’L) trig, 2N(Uβˆ’L) mul | +| CoG | (Uβˆ’L) mul, (Uβˆ’L) add | +| **Total** | **O(NΒ² + N(Uβˆ’L))** | + +### Batch Mode (SIMD Analysis) + +The DFT spectrum loop (Stage 4) is the dominant cost. For typical parameters (N=40, Uβˆ’L=22), each bar requires ~1,760 trig evaluations. Due to the adaptive coefficient state, vectorization is limited to within-period parallelization. + +### Quality Metrics + +| Metric | Value / Estimate | +|:----------------------|:-----------------------------| +| Spectral resolution | Higher than ACP for same N | +| Adaptation speed | ~N bars to converge | +| Phase lag | ≀2 bars (constrained) | +| Noise sensitivity | Low (roofing + AGC) | + +## Validation + +| Test | Input | Expected | +|:---------------------- |:---------------------------------------- |:--------------------------- | +| Default parameters | Close series, 500 bars | Dominant cycle in [18, 40] | +| Pure sine (30-bar) | sin(2π·n/30) for 500 bars | Converges near 30 | +| Constant input | All values = 100.0 | Stable, no NaN/Inf | +| Short series ( + if lowerBound<8 + runtime.error("Lower bound must be at least 8") + if upperBound<=lowerBound + runtime.error("Upper bound must be greater than lower bound") + if dataLength<4 + runtime.error("Data length must be at least 4") + var array xx=array.new_float(0) + var array coef=array.new_float(0) + var array pwr=array.new_float(0) + var int storedLen=0 + var int storedUpper=0 + var bool configured=false + var float hp=0.0 + var float lp=0.0 + var float peak=0.1 + var float signal=0.0 + var float dom=0.0 + var float prevDom=0.0 + if not configured or storedLen!=dataLength or storedUpper!=upperBound + xx:=array.new_float(dataLength+1,0.0) + coef:=array.new_float(dataLength+1,0.0) + pwr:=array.new_float(upperBound+2,0.0) + storedLen:=dataLength + storedUpper:=upperBound + configured:=true + hp:=0.0 + lp:=0.0 + peak:=0.1 + signal:=0.0 + dom:=(lowerBound+upperBound)*0.5 + prevDom:=dom + float price=nz(source) + // Stage 1: Roofing filter β€” Highpass (Butterworth 2nd order) + float alphaHP=(math.cos(0.707*2.0*math.pi/float(upperBound))+math.sin(0.707*2.0*math.pi/float(upperBound))-1.0)/math.cos(0.707*2.0*math.pi/float(upperBound)) + hp:=math.pow(1.0-alphaHP/2.0,2.0)*(price-2.0*nz(price[1])+nz(price[2]))+2.0*(1.0-alphaHP)*nz(hp[1])-math.pow(1.0-alphaHP,2.0)*nz(hp[2]) + // Stage 1b: SuperSmoother lowpass + float a1=math.exp(-math.sqrt(2.0)*math.pi/float(lowerBound)) + float b1=2.0*a1*math.cos(math.sqrt(2.0)*math.pi/float(lowerBound)) + float ssC2=b1 + float ssC3=-(a1*a1) + float ssC1=1.0-ssC2-ssC3 + lp:=ssC1*(hp+nz(hp[1]))*0.5+ssC2*nz(lp[1])+ssC3*nz(lp[2]) + // Stage 2: AGC normalization + peak:=0.991*peak + if math.abs(lp)>peak + peak:=math.abs(lp) + signal:=peak>0.0?lp/peak:0.0 + // Stage 3: Griffiths adaptive predictor + // Shift data buffer + for count=dataLength to 1 + array.set(xx,count,count>1?array.get(xx,count-1):0.0) + array.set(xx,0,signal) + // Compute signal power + float sigPower=0.0 + for count=0 to dataLength-1 + float xVal=array.get(xx,count) + sigPower+=xVal*xVal + sigPower/=float(dataLength) + // Convergence factor + float mu=sigPower>0.0?0.25/(sigPower*float(dataLength)):0.0 + // Predict + float xBar=0.0 + for count=1 to dataLength + xBar+=array.get(coef,count)*array.get(xx,count) + // Error + coefficient update + float err=array.get(xx,0)-xBar + for count=1 to dataLength + float c=array.get(coef,count)+mu*err*array.get(xx,count) + array.set(coef,count,c) + // Stage 4: Spectrum from coefficients + float maxPwrLocal=0.0 + for period=lowerBound to upperBound + float realPart=0.0 + float imagPart=0.0 + for count=1 to dataLength + float angle=2.0*math.pi*float(count)/float(period) + realPart+=array.get(coef,count)*math.cos(angle) + imagPart+=array.get(coef,count)*math.sin(angle) + float p=realPart*realPart+imagPart*imagPart + array.set(pwr,period,p) + if p>maxPwrLocal + maxPwrLocal:=p + // Stage 5: Dominant cycle via center of gravity + float spx=0.0 + float sp=0.0 + for period=lowerBound to upperBound + float p=maxPwrLocal>0.0?array.get(pwr,period)/maxPwrLocal:0.0 + if p>=0.5 + spx+=float(period)*p + sp+=p + float rawDom=sp>0.0?spx/sp:dom + // Constrain change to Β±2 bars per update + float maxDelta=2.0 + if rawDom-prevDom>maxDelta + rawDom:=prevDom+maxDelta + if prevDom-rawDom>maxDelta + rawDom:=prevDom-maxDelta + dom:=math.max(float(lowerBound),math.min(float(upperBound),rawDom)) + prevDom:=dom + // Two-bar prediction for trigger + float xPred=0.0 + if dataLength>2 + for count=1 to dataLength-2 + xPred+=array.get(coef,count)*array.get(xx,count) + [dom,signal,xPred] + +// ---------- Main loop ---------- +i_source=input.source(close,"Source") +i_lower=input.int(18,"Lower Bound",minval=8,maxval=200) +i_upper=input.int(40,"Upper Bound",minval=10,maxval=500) +i_length=input.int(40,"Data Length",minval=4,maxval=200) +[dominantCycle,sig,pred]=lpf(i_source,i_lower,i_upper,i_length) +plot(dominantCycle,"Dominant Cycle",color=color.yellow,linewidth=2) +plot(sig*20+30,"Signal",color=color.lime,linewidth=1) +plot(pred*20+30,"Predict",color=color.red,linewidth=1) diff --git a/lib/cycles/lpf/tests/Lpf.Quantower.Tests.cs b/lib/cycles/lpf/tests/Lpf.Quantower.Tests.cs new file mode 100644 index 00000000..361f2e16 --- /dev/null +++ b/lib/cycles/lpf/tests/Lpf.Quantower.Tests.cs @@ -0,0 +1,289 @@ +using TradingPlatform.BusinessLayer; + +namespace QuanTAlib.Quantower.Tests; + +public class LpfIndicatorTests +{ + [Fact] + public void LpfIndicator_Constructor_SetsDefaults() + { + var indicator = new LpfIndicator(); + + Assert.Equal(18, indicator.LowerBound); + Assert.Equal(40, indicator.UpperBound); + Assert.Equal(40, indicator.DataLength); + Assert.Equal(SourceType.Close, indicator.Source); + Assert.True(indicator.ShowColdValues); + Assert.Equal("LPF - Ehlers Linear Predictive Filter", indicator.Name); + Assert.True(indicator.SeparateWindow); + Assert.True(indicator.OnBackGround); + } + + [Fact] + public void LpfIndicator_MinHistoryDepths_EqualsZero() + { + var indicator = new LpfIndicator(); + + Assert.Equal(0, LpfIndicator.MinHistoryDepths); + Assert.Equal(0, ((IWatchlistIndicator)indicator).MinHistoryDepths); + } + + [Fact] + public void LpfIndicator_ShortName_IncludesParameters() + { + var indicator = new LpfIndicator { LowerBound = 18, UpperBound = 40, DataLength = 40 }; + + Assert.True(indicator.ShortName.Contains("LPF", StringComparison.Ordinal)); + Assert.True(indicator.ShortName.Contains("18", StringComparison.Ordinal)); + Assert.True(indicator.ShortName.Contains("40", StringComparison.Ordinal)); + } + + [Fact] + public void LpfIndicator_Initialize_CreatesInternalLpf() + { + var indicator = new LpfIndicator { LowerBound = 18, UpperBound = 40, DataLength = 40 }; + + indicator.Initialize(); + + // After init, line series should exist (Cycle + Signal + Predict) + Assert.Equal(3, indicator.LinesSeries.Count); + } + + [Fact] + public void LpfIndicator_ProcessUpdate_HistoricalBar_ComputesValue() + { + var indicator = new LpfIndicator { LowerBound = 18, UpperBound = 40, DataLength = 40 }; + indicator.Initialize(); + + var now = DateTime.UtcNow; + indicator.HistoricalData.AddBar(now, 100, 105, 95, 102); + + var args = new UpdateArgs(UpdateReason.HistoricalBar); + indicator.ProcessUpdate(args); + + Assert.Equal(1, indicator.LinesSeries[0].Count); + Assert.True(double.IsFinite(indicator.LinesSeries[0].GetValue(0))); + } + + [Fact] + public void LpfIndicator_ProcessUpdate_NewBar_ComputesValue() + { + var indicator = new LpfIndicator { LowerBound = 18, UpperBound = 40, DataLength = 40 }; + indicator.Initialize(); + + var now = DateTime.UtcNow; + indicator.HistoricalData.AddBar(now, 100, 105, 95, 102); + indicator.HistoricalData.AddBar(now.AddMinutes(1), 102, 108, 100, 106); + + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar)); + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.NewBar)); + + Assert.Equal(2, indicator.LinesSeries[0].Count); + } + + [Fact] + public void LpfIndicator_ProcessUpdate_NewTick_ProcessesWithoutError() + { + var indicator = new LpfIndicator { LowerBound = 18, UpperBound = 40, DataLength = 40 }; + indicator.Initialize(); + + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.NewTick)); + + Assert.NotNull(indicator); + } + + [Fact] + public void LpfIndicator_MultipleUpdates_ProducesCorrectSequence() + { + var indicator = new LpfIndicator { LowerBound = 18, UpperBound = 40, DataLength = 40 }; + indicator.Initialize(); + + var now = DateTime.UtcNow; + double[] closes = { 100, 102, 105, 103, 107, 110, 108, 112, 115, 113 }; + + foreach (var close in closes) + { + indicator.HistoricalData.AddBar(now, close, close + 2, close - 2, close); + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar)); + now = now.AddMinutes(1); + } + + for (int i = 0; i < closes.Length; i++) + { + Assert.True(double.IsFinite(indicator.LinesSeries[0].GetValue(closes.Length - 1 - i))); + } + } + + [Fact] + public void LpfIndicator_DifferentSourceTypes_Work() + { + var sources = new[] { SourceType.Open, SourceType.High, SourceType.Low, SourceType.Close, SourceType.HL2, SourceType.HLC3 }; + + foreach (var source in sources) + { + var indicator = new LpfIndicator { LowerBound = 18, UpperBound = 40, DataLength = 40, Source = source }; + indicator.Initialize(); + + var now = DateTime.UtcNow; + indicator.HistoricalData.AddBar(now, 100, 110, 90, 105); + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar)); + + Assert.True(double.IsFinite(indicator.LinesSeries[0].GetValue(0)), + $"Source {source} should produce finite value"); + } + } + + [Fact] + public void LpfIndicator_LowerBound_CanBeChanged() + { + var indicator = new LpfIndicator { LowerBound = 18 }; + + Assert.Equal(18, indicator.LowerBound); + + indicator.LowerBound = 10; + Assert.Equal(10, indicator.LowerBound); + } + + [Fact] + public void LpfIndicator_UpperBound_CanBeChanged() + { + var indicator = new LpfIndicator { UpperBound = 40 }; + + Assert.Equal(40, indicator.UpperBound); + + indicator.UpperBound = 100; + Assert.Equal(100, indicator.UpperBound); + } + + [Fact] + public void LpfIndicator_DataLength_CanBeChanged() + { + var indicator = new LpfIndicator { DataLength = 40 }; + + Assert.Equal(40, indicator.DataLength); + + indicator.DataLength = 60; + Assert.Equal(60, indicator.DataLength); + } + + [Fact] + public void LpfIndicator_Source_CanBeChanged() + { + var indicator = new LpfIndicator { Source = SourceType.Close }; + + Assert.Equal(SourceType.Close, indicator.Source); + + indicator.Source = SourceType.Open; + Assert.Equal(SourceType.Open, indicator.Source); + } + + [Fact] + public void LpfIndicator_ShowColdValues_CanBeChanged() + { + var indicator = new LpfIndicator { ShowColdValues = true }; + + Assert.True(indicator.ShowColdValues); + + indicator.ShowColdValues = false; + Assert.False(indicator.ShowColdValues); + } + + [Fact] + public void LpfIndicator_ShortName_UpdatesWhenParametersChange() + { + var indicator = new LpfIndicator { LowerBound = 18, UpperBound = 40, DataLength = 40 }; + string initialName = indicator.ShortName; + + Assert.True(initialName.Contains("18", StringComparison.Ordinal)); + Assert.True(initialName.Contains("40", StringComparison.Ordinal)); + + indicator.LowerBound = 10; + indicator.UpperBound = 60; + string updatedName = indicator.ShortName; + + Assert.True(updatedName.Contains("10", StringComparison.Ordinal)); + Assert.True(updatedName.Contains("60", StringComparison.Ordinal)); + } + + [Fact] + public void LpfIndicator_ProcessUpdate_IgnoresNonBarUpdates() + { + var indicator = new LpfIndicator { LowerBound = 18, UpperBound = 40, DataLength = 40 }; + indicator.Initialize(); + + var now = DateTime.UtcNow; + indicator.HistoricalData.AddBar(now, 100, 105, 95, 102); + + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar)); + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.NewTick)); + + Assert.NotNull(indicator); + } + + [Fact] + public void LpfIndicator_CycleSeries_HasCorrectProperties() + { + var indicator = new LpfIndicator(); + indicator.Initialize(); + + var lineSeries = indicator.LinesSeries[0]; + + Assert.Equal("Cycle", lineSeries.Name); + Assert.Equal(2, lineSeries.Width); + Assert.Equal(LineStyle.Solid, lineSeries.Style); + } + + [Fact] + public void LpfIndicator_SignalSeries_HasCorrectProperties() + { + var indicator = new LpfIndicator(); + indicator.Initialize(); + + var signalSeries = indicator.LinesSeries[1]; + + Assert.Equal("Signal", signalSeries.Name); + Assert.Equal(1, signalSeries.Width); + Assert.Equal(LineStyle.Solid, signalSeries.Style); + } + + [Fact] + public void LpfIndicator_PredictSeries_HasCorrectProperties() + { + var indicator = new LpfIndicator(); + indicator.Initialize(); + + var predictSeries = indicator.LinesSeries[2]; + + Assert.Equal("Predict", predictSeries.Name); + Assert.Equal(1, predictSeries.Width); + Assert.Equal(LineStyle.Dot, predictSeries.Style); + } + + [Fact] + public void LpfIndicator_SineWave_ProducesFiniteValues() + { + var indicator = new LpfIndicator { LowerBound = 10, UpperBound = 50, DataLength = 50 }; + indicator.Initialize(); + + var now = DateTime.UtcNow; + const int knownPeriod = 30; + + for (int i = 0; i < 200; i++) + { + double price = 100.0 + 10.0 * Math.Sin(2.0 * Math.PI * i / knownPeriod); + indicator.HistoricalData.AddBar(now.AddMinutes(i), price, price + 1, price - 1, price); + indicator.ProcessUpdate(new UpdateArgs(UpdateReason.HistoricalBar)); + } + + double cycleValue = indicator.LinesSeries[0].GetValue(0); + Assert.InRange(cycleValue, 10, 50); + } + + [Fact] + public void LpfIndicator_SourceCodeLink_IsValid() + { + var indicator = new LpfIndicator(); + Assert.Contains("github.com", indicator.SourceCodeLink); + Assert.Contains("Lpf.Quantower.cs", indicator.SourceCodeLink); + } +} diff --git a/lib/cycles/lpf/tests/Lpf.Tests.cs b/lib/cycles/lpf/tests/Lpf.Tests.cs new file mode 100644 index 00000000..4499d549 --- /dev/null +++ b/lib/cycles/lpf/tests/Lpf.Tests.cs @@ -0,0 +1,507 @@ +using Xunit; + +namespace QuanTAlib.Tests; + +public class LpfTests +{ + private const double Tolerance = 1e-9; + + #region Constructor Tests + + [Fact] + public void Constructor_DefaultParameters_SetsProperties() + { + var lpf = new Lpf(); + + Assert.Equal("LPF(18,40,40)", lpf.Name); + Assert.False(lpf.IsHot); + Assert.Equal(18, lpf.LowerBound); + Assert.Equal(40, lpf.UpperBound); + Assert.Equal(40, lpf.DataLength); + } + + [Fact] + public void Constructor_CustomParameters_SetsProperties() + { + var lpf = new Lpf(lowerBound: 10, upperBound: 60, dataLength: 50); + + Assert.Equal("LPF(10,60,50)", lpf.Name); + Assert.Equal(10, lpf.LowerBound); + Assert.Equal(60, lpf.UpperBound); + Assert.Equal(50, lpf.DataLength); + } + + [Theory] + [InlineData(7)] + [InlineData(0)] + [InlineData(-1)] + public void Constructor_InvalidLowerBound_ThrowsArgumentOutOfRange(int lower) + { + var ex = Assert.Throws(() => new Lpf(lower, 40)); + Assert.Equal("lowerBound", ex.ParamName); + } + + [Theory] + [InlineData(18, 18)] + [InlineData(18, 10)] + [InlineData(20, 20)] + public void Constructor_UpperBoundNotGreaterThanLower_ThrowsArgumentOutOfRange(int lower, int upper) + { + var ex = Assert.Throws(() => new Lpf(lower, upper)); + Assert.Equal("upperBound", ex.ParamName); + } + + [Theory] + [InlineData(3)] + [InlineData(0)] + [InlineData(-1)] + public void Constructor_InvalidDataLength_ThrowsArgumentOutOfRange(int len) + { + var ex = Assert.Throws(() => new Lpf(18, 40, len)); + Assert.Equal("dataLength", ex.ParamName); + } + + [Fact] + public void Constructor_WithNullSource_ThrowsArgumentNullException() + { + Assert.Throws(() => new Lpf(null!, 18, 40)); + } + + [Fact] + public void Constructor_WithValidSource_Subscribes() + { + var source = new TSeries(); + var lpf = new Lpf(source, 18, 40); + + source.Add(new TValue(DateTime.UtcNow, 100.0)); + Assert.NotEqual(default, lpf.Last); + } + + #endregion + + #region Basic Calculation Tests + + [Fact] + public void Update_ReturnsValidTValue() + { + var lpf = new Lpf(); + var result = lpf.Update(new TValue(DateTime.UtcNow, 100.0)); + + Assert.True(double.IsFinite(result.Value)); + } + + [Fact] + public void Update_AfterWarmup_IsHotTrue() + { + var lpf = new Lpf(); + + var gbm = new GBM(seed: 42); + var bars = gbm.Fetch(200, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1)); + + foreach (var bar in bars) + { + lpf.Update(new TValue(bar.Time, bar.Close)); + } + + Assert.True(lpf.IsHot); + } + + [Fact] + public void Update_DominantCycle_WithinRange() + { + var lpf = new Lpf(18, 40, 40); + + var gbm = new GBM(seed: 42); + var bars = gbm.Fetch(500, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1)); + + foreach (var bar in bars) + { + lpf.Update(new TValue(bar.Time, bar.Close)); + } + + Assert.InRange(lpf.DominantCycle, 18, 40); + } + + [Fact] + public void Update_Signal_WithinUnitRange() + { + var lpf = new Lpf(); + + var gbm = new GBM(seed: 42); + var bars = gbm.Fetch(500, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1)); + + foreach (var bar in bars) + { + lpf.Update(new TValue(bar.Time, bar.Close)); + } + + // AGC normalization should keep signal within [-1, 1] + Assert.InRange(lpf.Signal, -1.0, 1.0); + } + + [Fact] + public void Update_InitialValue_WithinBounds() + { + var lpf = new Lpf(18, 40, 40); + + var result = lpf.Update(new TValue(DateTime.UtcNow, 100.0)); + + Assert.True(result.Value >= 18 && result.Value <= 40); + } + + [Fact] + public void Update_PureSine_ConvergesNearTruePeriod() + { + int truePeriod = 30; + var lpf = new Lpf(lowerBound: 10, upperBound: 50, dataLength: 50); + + // Feed a pure sine wave with known period + for (int i = 0; i < 500; i++) + { + double val = Math.Sin(2.0 * Math.PI * i / truePeriod); + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(i), val)); + } + + // Should converge reasonably close to true period + // Allow generous tolerance since LPF needs time to adapt + Assert.InRange(lpf.DominantCycle, truePeriod - 10, truePeriod + 10); + } + + #endregion + + #region Bar Correction Tests + + [Fact] + public void Update_IsNewTrue_AdvancesState() + { + var lpf = new Lpf(); + + lpf.Update(new TValue(DateTime.UtcNow, 100.0), isNew: true); + var first = lpf.Last.Value; + + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(1), 110.0), isNew: true); + var second = lpf.Last.Value; + + Assert.True(double.IsFinite(first) && double.IsFinite(second)); + } + + [Fact] + public void Update_IsNewFalse_ReplacesCurrentBar() + { + var lpf = new Lpf(); + + for (int i = 0; i < 100; i++) + { + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(i), 100.0 + Math.Sin(i * 0.1) * 10), isNew: true); + } + + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(100), 110.0), isNew: true); + var beforeCorrection = lpf.Last.Value; + + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(100), 90.0), isNew: false); + var afterCorrection = lpf.Last.Value; + + Assert.True(double.IsFinite(beforeCorrection) && double.IsFinite(afterCorrection)); + } + + [Fact] + public void Update_MultipleCorrections_RestoresToSnapshot() + { + var lpf = new Lpf(); + + for (int i = 0; i < 100; i++) + { + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(i), 100.0 + i), isNew: true); + } + + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(100), 150.0), isNew: true); + var originalValue = lpf.Last.Value; + + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(100), 160.0), isNew: false); + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(100), 140.0), isNew: false); + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(100), 150.0), isNew: false); + var restoredValue = lpf.Last.Value; + + Assert.Equal(originalValue, restoredValue, Tolerance); + } + + #endregion + + #region Reset Tests + + [Fact] + public void Reset_ClearsState() + { + var lpf = new Lpf(); + + for (int i = 0; i < 200; i++) + { + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(i), 100.0 + i)); + } + + Assert.True(lpf.IsHot); + + lpf.Reset(); + + Assert.False(lpf.IsHot); + Assert.Equal(default, lpf.Last); + } + + [Fact] + public void Reset_AllowsReuse() + { + var lpf = new Lpf(); + + for (int i = 0; i < 200; i++) + { + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(i), 100.0 + Math.Sin(i * 0.1) * 10)); + } + var firstResult = lpf.Last.Value; + + lpf.Reset(); + + for (int i = 0; i < 200; i++) + { + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(i), 100.0 + Math.Sin(i * 0.1) * 10)); + } + var secondResult = lpf.Last.Value; + + Assert.Equal(firstResult, secondResult, Tolerance); + } + + #endregion + + #region NaN/Infinity Handling Tests + + [Fact] + public void Update_NaN_UsesLastValidValue() + { + var lpf = new Lpf(); + + lpf.Update(new TValue(DateTime.UtcNow, 100.0)); + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(1), double.NaN)); + + Assert.True(double.IsFinite(lpf.Last.Value)); + } + + [Fact] + public void Update_Infinity_UsesLastValidValue() + { + var lpf = new Lpf(); + + lpf.Update(new TValue(DateTime.UtcNow, 100.0)); + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(1), double.PositiveInfinity)); + + Assert.True(double.IsFinite(lpf.Last.Value)); + } + + [Fact] + public void Update_NegativeInfinity_UsesLastValidValue() + { + var lpf = new Lpf(); + + lpf.Update(new TValue(DateTime.UtcNow, 100.0)); + lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(1), double.NegativeInfinity)); + + Assert.True(double.IsFinite(lpf.Last.Value)); + } + + #endregion + + #region Consistency Tests + + [Theory] + [InlineData(42)] + [InlineData(123)] + [InlineData(456)] + public void Update_Deterministic_AcrossSeeds(int seed) + { + var lpf1 = new Lpf(); + var lpf2 = new Lpf(); + + var gbm = new GBM(seed: seed); + var bars = gbm.Fetch(300, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1)); + + foreach (var bar in bars) + { + var input = new TValue(bar.Time, bar.Close); + lpf1.Update(input); + lpf2.Update(input); + } + + Assert.Equal(lpf1.Last.Value, lpf2.Last.Value, Tolerance); + Assert.Equal(lpf1.DominantCycle, lpf2.DominantCycle, Tolerance); + Assert.Equal(lpf1.Signal, lpf2.Signal, Tolerance); + Assert.Equal(lpf1.Predict, lpf2.Predict, Tolerance); + } + + [Fact] + public void Batch_MatchesStreaming() + { + var gbm = new GBM(seed: 42); + var bars = gbm.Fetch(200, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1)); + + // TSeries from bars + var source = new TSeries(); + foreach (var bar in bars) + { + source.Add(new TValue(bar.Time, bar.Close)); + } + + // Batch + var batchResult = Lpf.Batch(source, 18, 40, 40); + + // Streaming + var lpf = new Lpf(18, 40, 40); + var streamResults = new List(); + foreach (var bar in bars) + { + var r = lpf.Update(new TValue(bar.Time, bar.Close)); + streamResults.Add(r.Value); + } + + Assert.Equal(streamResults.Count, batchResult.Count); + for (int i = 0; i < streamResults.Count; i++) + { + Assert.Equal(streamResults[i], batchResult[i].Value, Tolerance); + } + } + + [Fact] + public void BatchSpan_MatchesStreaming() + { + var gbm = new GBM(seed: 42); + var bars = gbm.Fetch(200, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1)); + + double[] input = bars.Select(b => b.Close).ToArray(); + double[] output = new double[input.Length]; + + Lpf.Batch(input, output, 18, 40, 40); + + var lpf = new Lpf(18, 40, 40); + for (int i = 0; i < input.Length; i++) + { + var r = lpf.Update(new TValue(DateTime.MinValue, input[i])); + Assert.Equal(r.Value, output[i], Tolerance); + } + } + + #endregion + + #region Constant Input Tests + + [Fact] + public void Update_ConstantInput_NoNaNOrInf() + { + var lpf = new Lpf(); + + for (int i = 0; i < 200; i++) + { + var result = lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(i), 100.0)); + Assert.True(double.IsFinite(result.Value), $"Non-finite result at bar {i}: {result.Value}"); + } + } + + [Fact] + public void Update_ZeroInput_NoNaNOrInf() + { + var lpf = new Lpf(); + + for (int i = 0; i < 200; i++) + { + var result = lpf.Update(new TValue(DateTime.UtcNow.AddSeconds(i), 0.0)); + Assert.True(double.IsFinite(result.Value), $"Non-finite result at bar {i}: {result.Value}"); + } + } + + #endregion + + #region Calculate + Dispose Tests + + [Fact] + public void Calculate_ReturnsResultsAndIndicator() + { + var gbm = new GBM(seed: 42); + var bars = gbm.Fetch(200, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1)); + + var source = new TSeries(); + foreach (var bar in bars) + { + source.Add(new TValue(bar.Time, bar.Close)); + } + + var (results, indicator) = Lpf.Calculate(source, 18, 40, 40); + + Assert.Equal(200, results.Count); + Assert.True(indicator.IsHot); + } + + [Fact] + public void Dispose_UnsubscribesFromSource() + { + var source = new TSeries(); + var lpf = new Lpf(source, 18, 40, 40); + + source.Add(new TValue(DateTime.UtcNow, 100.0)); + Assert.NotEqual(default, lpf.Last); + + lpf.Dispose(); + + // After dispose, adding to source should not update lpf + var lastBefore = lpf.Last; + source.Add(new TValue(DateTime.UtcNow.AddSeconds(1), 200.0)); + Assert.Equal(lastBefore, lpf.Last); + } + + #endregion + + #region DominantCycle Rate Constraint Tests + + [Fact] + public void Update_DominantCycle_ChangeConstrainedToTwo() + { + var lpf = new Lpf(10, 50, 40); + + // Feed data and track DC changes + var gbm = new GBM(seed: 42); + var bars = gbm.Fetch(500, DateTime.UtcNow.Ticks, TimeSpan.FromMinutes(1)); + + double prevDC = 0; + bool first = true; + foreach (var bar in bars) + { + lpf.Update(new TValue(bar.Time, bar.Close)); + double dc = lpf.DominantCycle; + if (!first) + { + double delta = Math.Abs(dc - prevDC); + Assert.True(delta <= 2.0 + 1e-10, $"DC changed by {delta} > 2.0"); + } + prevDC = dc; + first = false; + } + } + + #endregion + + #region Prime Tests + + [Fact] + public void Prime_SetsState() + { + var lpf = new Lpf(); + + double[] data = new double[200]; + for (int i = 0; i < 200; i++) + { + data[i] = 100.0 + Math.Sin(i * 0.1) * 10; + } + + lpf.Prime(data); + + Assert.True(lpf.IsHot); + Assert.True(double.IsFinite(lpf.Last.Value)); + } + + #endregion +} diff --git a/python/quantalib/_bridge.py b/python/quantalib/_bridge.py index 2371da0c..c8893529 100644 --- a/python/quantalib/_bridge.py +++ b/python/quantalib/_bridge.py @@ -573,6 +573,7 @@ HAS_DSP = _bind("qtl_dsp", [_dp, _ci, _dp, _ci]) HAS_CCOR = _bind("qtl_ccor", [_dp, _ci, _dp, _ci, _cd]) HAS_EBSW = _bind("qtl_ebsw", [_dp, _ci, _dp, _ci, _ci]) HAS_ACP = _bind("qtl_acp", [_dp, _ci, _dp, _ci, _ci, _ci, _ci]) +HAS_LPF = _bind("qtl_lpf", [_dp, _ci, _dp, _ci, _ci, _ci]) HAS_AMFM = _bind("qtl_amfm", [_dp, _dp, _ci, _dp, _dp, _ci]) # ── Numerics (Exports.cs β€” manual) ── diff --git a/python/quantalib/cycles.py b/python/quantalib/cycles.py index 46b21f38..51ccf136 100644 --- a/python/quantalib/cycles.py +++ b/python/quantalib/cycles.py @@ -13,6 +13,7 @@ __all__ = [ "ht_dcphase", "ht_phasor", "ht_sine", + "lpf", "lunar", "solar", "ssfdsp", @@ -79,6 +80,16 @@ def ht_sine(close: object, offset: int = 0, **kwargs) -> object: return _wrap_multi({"sine": sine, "leadSine": leadSine}, idx, "cycles", offset) +def lpf(close: object, lower_bound: int = 18, upper_bound: int = 40, + data_length: int = 40, offset: int = 0, **kwargs) -> object: + """Ehlers Linear Predictive Filter (dominant cycle).""" + lower_bound = int(lower_bound); upper_bound = int(upper_bound) + data_length = int(data_length); offset = int(offset) + src, idx = _arr(close); n = len(src); dst = _out(n) + _check(_lib.qtl_lpf(_ptr(src), n, _ptr(dst), lower_bound, upper_bound, data_length)) + return _wrap(dst, idx, f"LPF_{lower_bound}_{upper_bound}", "cycles", offset) + + def lunar(close: object, offset: int = 0, **kwargs) -> object: """Lunar Cycle.""" offset = int(offset) diff --git a/python/src/Exports.cs b/python/src/Exports.cs index 3fe4220a..4f3ba2e5 100644 --- a/python/src/Exports.cs +++ b/python/src/Exports.cs @@ -1546,6 +1546,16 @@ public static unsafe partial class Exports catch { return StatusCodes.QTL_ERR_INTERNAL; } } + // Lpf: Pattern A (source β†’ output, int lowerBound, int upperBound, int dataLength) + [UnmanagedCallersOnly(EntryPoint = "qtl_lpf")] + public static int QtlLpf(double* src, int n, double* dst, int lowerBound, int upperBound, int dataLength) + { + int v = Chk1(src, dst, n); if (v != 0) return v; + v = ChkPeriod(lowerBound); if (v != 0) return v; + try { Lpf.Batch(Src(src, n), Dst(dst, n), lowerBound, upperBound, dataLength); return StatusCodes.QTL_OK; } + catch { return StatusCodes.QTL_ERR_INTERNAL; } + } + // ═══════════════════════════════════════════════════════════════════════ // Β§8.14 Numerics / transforms // ═══════════════════════════════════════════════════════════════════════