wickra_core/indicators/
hilbert_dominant_cycle.rs

1//! Ehlers Hilbert Transform Dominant Cycle period estimator.
2#![allow(clippy::manual_clamp)]
3
4use std::f64::consts::PI;
5
6use crate::traits::Indicator;
7
8/// Ehlers' Hilbert Transform–based Dominant Cycle period estimator.
9///
10/// Decomposes price into in-phase and quadrature components via Ehlers'
11/// truncated Hilbert transform, then derives the instantaneous phase. The
12/// dominant cycle period is recovered from the phase rate of change,
13/// rate-limited, clamped and EMA-smoothed (0.2/0.8, then 0.33/0.67). From *Rocket Science for Traders* (Ehlers 2001, ch. 7),
14/// implementation aligned with the formulation used in TA-Lib's `HT_DCPERIOD`.
15///
16/// The output is clamped to the band `[6, 50]` bars, which Ehlers identifies
17/// as the meaningful tradable cycle range. The estimator emits its first
18/// value after ~50 inputs as the moving-average chain fills.
19///
20/// # Example
21///
22/// ```
23/// use wickra_core::{Indicator, HilbertDominantCycle};
24///
25/// let mut ht = HilbertDominantCycle::new();
26/// let mut last = None;
27/// for i in 0..200 {
28///     last = ht.update(100.0 + (f64::from(i) * 0.4).sin() * 5.0);
29/// }
30/// assert!(last.is_some());
31/// ```
32#[derive(Debug, Clone, Default)]
33pub struct HilbertDominantCycle {
34    // Raw input window for the 4-bar WMA.
35    price_buf: Vec<f64>,
36    // WMA-smoothed price history feeding the Hilbert detrender taps.
37    smooth_buf: Vec<f64>,
38    // Detrender / Q1 / I1 ring history (need 6 prior).
39    detrender_buf: Vec<f64>,
40    q1_buf: Vec<f64>,
41    i1_buf: Vec<f64>,
42    // Smoothed I/Q lines for phase computation.
43    prev_i2: f64,
44    prev_q2: f64,
45    prev_re: f64,
46    prev_im: f64,
47    prev_period: f64,
48    prev_smooth_period: f64,
49    count: usize,
50    last_value: Option<f64>,
51}
52
53impl HilbertDominantCycle {
54    /// Construct a new dominant cycle estimator.
55    pub fn new() -> Self {
56        Self::default()
57    }
58
59    /// Current period estimate if available.
60    pub const fn value(&self) -> Option<f64> {
61        self.last_value
62    }
63}
64
65impl Indicator for HilbertDominantCycle {
66    type Input = f64;
67    type Output = f64;
68
69    fn update(&mut self, input: f64) -> Option<f64> {
70        if !input.is_finite() {
71            return None;
72        }
73        self.count += 1;
74
75        // 4-bar weighted moving average of the input (smoothed price).
76        // Ehlers: (4*x[0] + 3*x[1] + 2*x[2] + x[3]) / 10.
77        Self::push_front(&mut self.price_buf, input, 4);
78        if self.price_buf.len() < 4 {
79            return None;
80        }
81        let smooth = (4.0 * self.price_buf[0]
82            + 3.0 * self.price_buf[1]
83            + 2.0 * self.price_buf[2]
84            + self.price_buf[3])
85            / 10.0;
86        Self::push_front(&mut self.smooth_buf, smooth, 7);
87
88        // Adaptive coefficient based on the previous period estimate.
89        let period = self.prev_period.max(6.0).min(50.0);
90        let adj = 0.075 * period + 0.54;
91
92        // We need the smooth buffer to hold ≥ 7 samples for the Hilbert taps.
93        if self.smooth_buf.len() < 7 {
94            return None;
95        }
96
97        // Ehlers' Hilbert transform of `smooth` (using current + 2/4/6 lags).
98        let s0 = smooth;
99        let s2 = self.smooth_buf[2];
100        let s4 = self.smooth_buf[4];
101        let s6 = self.smooth_buf[6];
102        let detrender = (0.0962 * s0 + 0.5769 * s2 - 0.5769 * s4 - 0.0962 * s6) * adj;
103        Self::push_front(&mut self.detrender_buf, detrender, 7);
104
105        if self.detrender_buf.len() < 7 {
106            return None;
107        }
108        // In-phase and quadrature components.
109        let q1 = (0.0962 * self.detrender_buf[0] + 0.5769 * self.detrender_buf[2]
110            - 0.5769 * self.detrender_buf[4]
111            - 0.0962 * self.detrender_buf[6])
112            * adj;
113        let i1 = self.detrender_buf[3];
114
115        Self::push_front(&mut self.q1_buf, q1, 7);
116        Self::push_front(&mut self.i1_buf, i1, 7);
117        if self.q1_buf.len() < 7 || self.i1_buf.len() < 7 {
118            return None;
119        }
120
121        // Advance the phase 90 deg via a second Hilbert pass.
122        let ji = (0.0962 * self.i1_buf[0] + 0.5769 * self.i1_buf[2]
123            - 0.5769 * self.i1_buf[4]
124            - 0.0962 * self.i1_buf[6])
125            * adj;
126        let jq = (0.0962 * self.q1_buf[0] + 0.5769 * self.q1_buf[2]
127            - 0.5769 * self.q1_buf[4]
128            - 0.0962 * self.q1_buf[6])
129            * adj;
130
131        // Phasor smoothing.
132        let mut i2 = i1 - jq;
133        let mut q2 = q1 + ji;
134        i2 = 0.2 * i2 + 0.8 * self.prev_i2;
135        q2 = 0.2 * q2 + 0.8 * self.prev_q2;
136
137        // Homodyne discriminator.
138        let mut re = i2 * self.prev_i2 + q2 * self.prev_q2;
139        let mut im = i2 * self.prev_q2 - q2 * self.prev_i2;
140        re = 0.2 * re + 0.8 * self.prev_re;
141        im = 0.2 * im + 0.8 * self.prev_im;
142
143        self.prev_i2 = i2;
144        self.prev_q2 = q2;
145        self.prev_re = re;
146        self.prev_im = im;
147
148        let mut new_period = if im.abs() > f64::EPSILON && re.abs() > f64::EPSILON {
149            2.0 * PI / im.atan2(re)
150        } else {
151            self.prev_period
152        };
153        // Rate-of-change clamp per Ehlers.
154        new_period = new_period.min(1.5 * self.prev_period);
155        new_period = new_period.max(0.67 * self.prev_period);
156        new_period = new_period.clamp(6.0, 50.0);
157
158        // EMA smoothing of the period.
159        self.prev_period = 0.2 * new_period + 0.8 * self.prev_period;
160        // Second smoothing step (TA-Lib uses 0.33/0.67).
161        self.prev_smooth_period = 0.33 * self.prev_period + 0.67 * self.prev_smooth_period;
162
163        if self.count < 50 {
164            return None;
165        }
166        self.last_value = Some(self.prev_smooth_period);
167        Some(self.prev_smooth_period)
168    }
169
170    fn reset(&mut self) {
171        self.price_buf.clear();
172        self.smooth_buf.clear();
173        self.detrender_buf.clear();
174        self.q1_buf.clear();
175        self.i1_buf.clear();
176        self.prev_i2 = 0.0;
177        self.prev_q2 = 0.0;
178        self.prev_re = 0.0;
179        self.prev_im = 0.0;
180        self.prev_period = 0.0;
181        self.prev_smooth_period = 0.0;
182        self.count = 0;
183        self.last_value = None;
184    }
185
186    #[inline]
187    fn warmup_period(&self) -> usize {
188        50
189    }
190
191    #[inline]
192    fn is_ready(&self) -> bool {
193        self.last_value.is_some()
194    }
195
196    #[inline]
197    fn name(&self) -> &'static str {
198        "HilbertDominantCycle"
199    }
200}
201
202impl HilbertDominantCycle {
203    /// Push `v` at the front of `buf`, capping the length at `cap`.
204    fn push_front(buf: &mut Vec<f64>, v: f64, cap: usize) {
205        buf.insert(0, v);
206        if buf.len() > cap {
207            buf.truncate(cap);
208        }
209    }
210}
211
212#[cfg(test)]
213mod tests {
214    use super::*;
215    use crate::traits::BatchExt;
216
217    #[test]
218    fn accessors_and_metadata() {
219        let mut ht = HilbertDominantCycle::new();
220        assert_eq!(ht.warmup_period(), 50);
221        assert_eq!(ht.name(), "HilbertDominantCycle");
222        assert!(!ht.is_ready());
223        assert!(ht.value().is_none());
224        for i in 0..120 {
225            ht.update(100.0 + (f64::from(i) * 0.3).sin() * 5.0);
226        }
227        assert!(ht.is_ready());
228        assert!(ht.value().is_some());
229    }
230
231    #[test]
232    fn output_within_clamp_band() {
233        let mut ht = HilbertDominantCycle::new();
234        let prices: Vec<f64> = (0..200)
235            .map(|i| 100.0 + (f64::from(i) * 0.4).sin() * 5.0)
236            .collect();
237        let out = ht.batch(&prices);
238        for v in out.iter().flatten() {
239            assert!((6.0..=50.0).contains(v), "period {v} outside [6, 50]");
240        }
241    }
242
243    #[test]
244    fn batch_equals_streaming() {
245        let prices: Vec<f64> = (0..200)
246            .map(|i| 100.0 + (f64::from(i) * 0.3).sin() * 5.0)
247            .collect();
248        let mut a = HilbertDominantCycle::new();
249        let mut b = HilbertDominantCycle::new();
250        let batch = a.batch(&prices);
251        let streamed: Vec<_> = prices.iter().map(|p| b.update(*p)).collect();
252        assert_eq!(batch, streamed);
253    }
254
255    #[test]
256    fn ignores_non_finite_input() {
257        let mut ht = HilbertDominantCycle::new();
258        let prices: Vec<f64> = (0..120)
259            .map(|i| 100.0 + (f64::from(i) * 0.4).sin() * 5.0)
260            .collect();
261        ht.batch(&prices);
262        let before = ht.value();
263        assert!(before.is_some());
264        assert_eq!(ht.update(f64::NAN), None);
265    }
266
267    #[test]
268    fn reset_clears_state() {
269        let mut ht = HilbertDominantCycle::new();
270        let prices: Vec<f64> = (0..120)
271            .map(|i| 100.0 + (f64::from(i) * 0.4).sin() * 5.0)
272            .collect();
273        ht.batch(&prices);
274        assert!(ht.is_ready());
275        ht.reset();
276        assert!(!ht.is_ready());
277        assert!(ht.value().is_none());
278    }
279
280    use crate::traits::BatchNanExt;
281    use approx::assert_relative_eq;
282
283    fn sine_prices(n: u32) -> Vec<f64> {
284        (0..n)
285            .map(|i| 100.0 + (f64::from(i) * 0.4).sin() * 5.0)
286            .collect()
287    }
288
289    #[test]
290    fn first_value_lands_exactly_at_warmup() {
291        let mut ht = HilbertDominantCycle::new();
292        let out = ht.batch(&sine_prices(120));
293        let warmup = ht.warmup_period();
294        assert!(out[..warmup - 1].iter().all(Option::is_none));
295        assert!(out[warmup - 1].is_some());
296    }
297
298    #[test]
299    fn reset_replays_identically() {
300        let prices = sine_prices(150);
301        let fresh = HilbertDominantCycle::new().batch(&prices);
302        let mut ht = HilbertDominantCycle::new();
303        let first = ht.batch(&prices);
304        ht.reset();
305        let second = ht.batch(&prices);
306        assert_eq!(first, fresh);
307        assert_eq!(second, fresh);
308    }
309
310    #[test]
311    fn batch_nan_paths_match_streaming_bitwise() {
312        let prices = sine_prices(150);
313        let mut out = vec![0.0; prices.len()];
314        HilbertDominantCycle::new().batch_nan_into(&prices, &mut out);
315        let nan = HilbertDominantCycle::new().batch_nan(&prices);
316        let fast = HilbertDominantCycle::new().batch_fast(&prices);
317        let mut stream = HilbertDominantCycle::new();
318        let expected: Vec<u64> = prices
319            .iter()
320            .map(|&p| stream.update(p).unwrap_or(f64::NAN).to_bits())
321            .collect();
322        assert!(out.iter().zip(&expected).all(|(v, e)| v.to_bits() == *e));
323        assert!(nan.iter().zip(&expected).all(|(v, e)| v.to_bits() == *e));
324        assert!(fast.iter().zip(&expected).all(|(v, e)| v.to_bits() == *e));
325    }
326
327    #[test]
328    fn wma_of_raw_inputs_feeds_detrender_taps() {
329        let mut ht = HilbertDominantCycle::new();
330        // After exactly 4 inputs the WMA is (4*40 + 3*30 + 2*20 + 10) / 10 = 30.
331        for p in [10.0, 20.0, 30.0, 40.0] {
332            assert_eq!(ht.update(p), None);
333        }
334        assert_eq!(ht.smooth_buf, vec![30.0]);
335        assert_eq!(ht.detrender_buf.len(), 0);
336
337        // A single spike of 10 at index 7 in an all-zero series gives smoothed
338        // values 4, 3, 2 at indices 7, 8, 9, so at index 9 the smooth history is
339        // [2, 3, 4, 0, 0, 0, 0]. With the period seed of 6, adj = 0.075*6 + 0.54
340        // = 0.99, and the detrender reads the SMOOTHED taps s0, s2, s4, s6:
341        //   (0.0962*2 + 0.5769*4 - 0.5769*0 - 0.0962*0) * 0.99 = 2.5 * 0.99 = 2.475.
342        // Raw-price taps would instead give 0.5769*10*0.99 = 5.711_31.
343        let mut ht = HilbertDominantCycle::new();
344        let mut series = [0.0; 10];
345        series[7] = 10.0;
346        let _ = ht.batch(&series);
347        assert_eq!(ht.smooth_buf, vec![2.0, 3.0, 4.0, 0.0, 0.0, 0.0, 0.0]);
348        assert_eq!(ht.detrender_buf.len(), 1);
349        assert_relative_eq!(ht.detrender_buf[0], 2.475, epsilon = 1e-12);
350    }
351
352    #[test]
353    fn constant_input_period_settles_at_lower_clamp() {
354        // A flat zero series makes re == im == 0, so the period falls back to
355        // its previous value, is clamped up to 6, and the EMA chain converges
356        // on 6 from below.
357        let mut ht = HilbertDominantCycle::new();
358        let out = ht.batch(&[0.0; 400]);
359        assert!(out
360            .iter()
361            .flatten()
362            .all(|v| v.is_finite() && *v <= 6.0 + 1e-9));
363        assert_relative_eq!(ht.value().unwrap(), 6.0, epsilon = 1e-9);
364
365        // A non-zero flat series stays finite and inside the clamp band.
366        let mut ht = HilbertDominantCycle::new();
367        let out = ht.batch(&[100.0; 400]);
368        assert!(out.iter().flatten().all(|v| v.is_finite() && *v <= 50.0));
369        assert_eq!(out.iter().flatten().count(), 400 - 49);
370    }
371
372    #[test]
373    fn non_finite_input_during_warmup_does_not_advance() {
374        let prices = sine_prices(120);
375        let mut ht = HilbertDominantCycle::new();
376        let mut out = Vec::new();
377        for (i, &p) in prices.iter().enumerate() {
378            if i == 10 {
379                assert_eq!(ht.update(f64::INFINITY), None);
380            }
381            out.push(ht.update(p));
382        }
383        assert_eq!(out, HilbertDominantCycle::new().batch(&prices));
384    }
385}
wickra_core/indicators/hilbert_dominant_cycle.rs

wickra_core/indicators/
hilbert_dominant_cycle.rs
