quantwave_core/indicators/
griffiths_spectrum.rs

1use crate::indicators::high_pass::HighPass;
2use crate::indicators::metadata::{IndicatorMetadata, ParamDef};
3use crate::indicators::super_smoother::SuperSmoother;
4use crate::traits::Next;
5use crate::utils::RingBuffer as VecDeque;
6use std::f64::consts::PI;
7
8/// Griffiths Spectrum
9///
10/// Based on John Ehlers' "Linear Predictive Filters And Instantaneous Frequency" (TASC January 2025).
11/// It computes the normalized power spectrum using Griffiths adaptive filter coefficients.
12#[derive(Debug, Clone)]
13pub struct GriffithsSpectrum {
14    lb: usize,
15    ub: usize,
16    length: usize,
17    mu: f64,
18    hp: HighPass,
19    ss: SuperSmoother,
20    peak: f64,
21    signal_window: VecDeque<f64>,
22    coef: Vec<f64>,
23}
24
25impl GriffithsSpectrum {
26    pub fn new(lower_bound: usize, upper_bound: usize, length: usize) -> Self {
27        Self {
28            lb: lower_bound,
29            ub: upper_bound,
30            length,
31            mu: 1.0 / (length as f64),
32            hp: HighPass::new(upper_bound),
33            ss: SuperSmoother::new(lower_bound),
34            peak: 0.1,
35            signal_window: VecDeque::with_capacity(length + 1),
36            coef: vec![0.0; length + 1],
37        }
38    }
39}
40
41impl Default for GriffithsSpectrum {
42    fn default() -> Self {
43        Self::new(18, 40, 40)
44    }
45}
46
47impl Next<f64> for GriffithsSpectrum {
48    type Output = Vec<f64>; // Power for each period from lb to ub
49
50    fn next(&mut self, input: f64) -> Self::Output {
51        let hp_val = self.hp.next(input);
52        let lp_val = self.ss.next(hp_val);
53
54        self.peak *= 0.991;
55        if lp_val.abs() > self.peak {
56            self.peak = lp_val.abs();
57        }
58
59        let signal = if self.peak != 0.0 {
60            lp_val / self.peak
61        } else {
62            0.0
63        };
64
65        self.signal_window.push_front(signal);
66        if self.signal_window.len() > self.length {
67            self.signal_window.pop_back();
68        }
69
70        let mut results = vec![0.0; self.ub - self.lb + 1];
71
72        if self.signal_window.len() < self.length {
73            return results;
74        }
75
76        let mut xx = vec![0.0; self.length + 1];
77        for (i, val) in xx.iter_mut().enumerate().skip(1).take(self.length) {
78            *val = self.signal_window[self.length - i];
79        }
80
81        let mut x_bar = 0.0;
82        for count in 1..=self.length {
83            x_bar += xx[self.length - count] * self.coef[count];
84        }
85
86        for count in 1..=self.length {
87            self.coef[count] += self.mu * (xx[self.length] - x_bar) * xx[self.length - count];
88        }
89
90        let mut max_pwr = 0.0;
91        let mut powers = Vec::with_capacity(self.ub - self.lb + 1);
92
93        for period_idx in self.lb..=self.ub {
94            let period = period_idx as f64;
95            let mut real = 0.0;
96            let mut imag = 0.0;
97
98            for count in 1..=self.length {
99                let angle = 2.0 * PI * (count as f64) / period;
100                real += self.coef[count] * angle.cos();
101                imag += self.coef[count] * angle.sin();
102            }
103
104            let denom = (1.0 - real).powi(2) + imag.powi(2);
105            let pwr = 0.1 / denom;
106
107            if pwr > max_pwr {
108                max_pwr = pwr;
109            }
110            powers.push(pwr);
111        }
112
113        if max_pwr != 0.0 {
114            for (i, pwr) in powers.into_iter().enumerate() {
115                results[i] = pwr / max_pwr;
116            }
117        }
118
119        results
120    }
121}
122
123pub const GRIFFITHS_SPECTRUM_METADATA: IndicatorMetadata = IndicatorMetadata {
124    name: "GriffithsSpectrum",
125    description: "Normalized power spectrum estimation using Griffiths adaptive filters.",
126    usage: "Use to generate a high-resolution periodogram for cycle analysis. Best visualized as a heatmap to identify and track multiple market cycles simultaneously.",
127    keywords: &["spectrum", "cycle", "ehlers", "dsp", "periodogram"],
128    ehlers_summary: "The Griffiths Spectrum is an adaptive spectral estimation method that provides higher resolution than a standard DFT for short data segments. It fits an all-pole model to the signal using an LMS algorithm, allowing for instantaneous frequency measurement without the windowing artifacts of FFT-based methods.",
129    params: &[
130        ParamDef {
131            name: "lower_bound",
132            default: "18",
133            description: "Lower period bound",
134        },
135        ParamDef {
136            name: "upper_bound",
137            default: "40",
138            description: "Upper period bound",
139        },
140        ParamDef {
141            name: "length",
142            default: "40",
143            description: "LMS filter length",
144        },
145    ],
146    formula_source: "https://github.com/lavs9/quantwave/blob/main/references/traderstipsreference/TRADERS’%20TIPS%20-%20JANUARY%202025.html",
147    formula_latex: r#"
148\[
149Pwr(P) = \frac{0.1}{(1 - \sum coef_i \cos(2\pi i/P))^2 + (\sum coef_i \sin(2\pi i/P))^2}
150\]
151\[
152Pwr_{norm}(P) = \frac{Pwr(P)}{\max(Pwr)}
153\]
154"#,
155    gold_standard_file: "griffiths_spectrum.json",
156    category: "Ehlers DSP",
157};
158
159#[cfg(test)]
160mod tests {
161    use super::*;
162    use crate::traits::Next;
163    // use crate::test_utils;
164    // use crate::test_utils::{load_gold_standard_vec, assert_indicator_parity_vec};
165    use proptest::prelude::*;
166
167    /*
168    #[test]
169    fn test_griffiths_spectrum_gold_standard() {
170        let case = load_gold_standard_vec("griffiths_spectrum");
171        let gs = GriffithsSpectrum::new(18, 40, 40);
172        assert_indicator_parity_vec(gs, &case.input, &case.expected);
173    }
174    */
175    // TODO: Restore test once griffiths_spectrum.json is recovered.
176
177    #[test]
178    fn test_griffiths_spectrum_basic() {
179        let mut gs = GriffithsSpectrum::new(18, 40, 40);
180        let inputs = vec![10.0, 11.0, 12.0, 13.0, 14.0, 15.0];
181        for input in inputs {
182            let res = gs.next(input);
183            assert_eq!(res.len(), 40 - 18 + 1);
184        }
185    }
186
187    proptest! {
188        #[test]
189        fn test_griffiths_spectrum_parity(
190            inputs in prop::collection::vec(1.0..100.0, 100..200),
191        ) {
192            let lb = 18;
193            let ub = 40;
194            let length = 40;
195            let mut gs = GriffithsSpectrum::new(lb, ub, length);
196            let streaming_results: Vec<Vec<f64>> = inputs.iter().map(|&x| gs.next(x)).collect();
197
198            // Batch implementation
199            let mut batch_results = Vec::with_capacity(inputs.len());
200            let mut hp = HighPass::new(ub);
201            let mut ss = SuperSmoother::new(lb);
202            let lp_vals: Vec<f64> = inputs.iter().map(|&x| ss.next(hp.next(x))).collect();
203
204            let mut peak = 0.1;
205            let mut signals = Vec::new();
206            let mut coef = vec![0.0; length + 1];
207            let mu = 1.0 / length as f64;
208
209            for (i, &lp_val) in lp_vals.iter().enumerate() {
210                peak *= 0.991;
211                if lp_val.abs() > peak {
212                    peak = lp_val.abs();
213                }
214                let signal = if peak != 0.0 { lp_val / peak } else { 0.0 };
215                signals.push(signal);
216
217                if signals.len() < length {
218                    batch_results.push(vec![0.0; ub - lb + 1]);
219                    continue;
220                }
221
222                let mut xx = vec![0.0; length + 1];
223                for j in 1..=length {
224                    xx[j] = signals[i - (length - j)];
225                }
226
227                let mut x_bar = 0.0;
228                for count in 1..=length {
229                    x_bar += xx[length - count] * coef[count];
230                }
231
232                for count in 1..=length {
233                    coef[count] += mu * (xx[length] - x_bar) * xx[length - count];
234                }
235
236                let mut powers = Vec::new();
237                let mut max_pwr = 0.0;
238                for period_idx in lb..=ub {
239                    let period = period_idx as f64;
240                    let mut real = 0.0;
241                    let mut imag = 0.0;
242                    for count in 1..=length {
243                        let angle = 2.0 * PI * (count as f64) / period;
244                        real += coef[count] * angle.cos();
245                        imag += coef[count] * angle.sin();
246                    }
247                    let denom = (1.0 - real).powi(2) + imag.powi(2);
248                    let pwr = 0.1 / denom;
249                    if pwr > max_pwr { max_pwr = pwr; }
250                    powers.push(pwr);
251                }
252
253                let norm_powers = if max_pwr != 0.0 {
254                    powers.into_iter().map(|p| p / max_pwr).collect()
255                } else {
256                    vec![0.0; ub - lb + 1]
257                };
258                batch_results.push(norm_powers);
259            }
260
261            for (s, b) in streaming_results.iter().zip(batch_results.iter()) {
262                for (sv, bv) in s.iter().zip(b.iter()) {
263                    approx::assert_relative_eq!(sv, bv, epsilon = 1e-10);
264                }
265            }
266        }
267    }
268}
quantwave_core/indicators/griffiths_spectrum.rs

quantwave_core/indicators/
griffiths_spectrum.rs
