quantwave_core/indicators/
fourier_transform.rs

1use crate::indicators::metadata::{IndicatorMetadata, ParamDef};
2use crate::traits::Next;
3use crate::indicators::high_pass::HighPass;
4use std::collections::VecDeque;
5
6/// Fourier Dominant Cycle
7///
8/// Based on John Ehlers' "Fourier Transform for Traders".
9/// Estimates the dominant cycle period using a Discrete Fourier Transform (DFT)
10/// with a Kay & Demeure resolution-enhancing transformation, followed by a
11/// center-of-gravity calculation.
12#[derive(Debug, Clone)]
13pub struct FourierDominantCycle {
14    window_len: usize,
15    hp: HighPass,
16    hp_history: [f64; 6],
17    cleaned_window: VecDeque<f64>,
18    count: usize,
19}
20
21impl FourierDominantCycle {
22    pub fn new(window_len: usize) -> Self {
23        Self {
24            window_len,
25            hp: HighPass::new(40),
26            hp_history: [0.0; 6],
27            cleaned_window: VecDeque::with_capacity(window_len),
28            count: 0,
29        }
30    }
31}
32
33impl Default for FourierDominantCycle {
34    fn default() -> Self {
35        Self::new(50)
36    }
37}
38
39impl Next<f64> for FourierDominantCycle {
40    type Output = f64;
41
42    fn next(&mut self, input: f64) -> Self::Output {
43        self.count += 1;
44        let hp_val = self.hp.next(input);
45
46        // FIR Smoothing: CleanedData = (HP + 2*HP[1] + 3*HP[2] + 3*HP[3] + 2*HP[4] + HP[5])/12
47        let cleaned = (hp_val 
48            + 2.0 * self.hp_history[0] 
49            + 3.0 * self.hp_history[1] 
50            + 3.0 * self.hp_history[2] 
51            + 2.0 * self.hp_history[3] 
52            + self.hp_history[4]) / 12.0;
53
54        // Shift HP history
55        for i in (1..6).rev() {
56            self.hp_history[i] = self.hp_history[i-1];
57        }
58        self.hp_history[0] = hp_val;
59
60        self.cleaned_window.push_front(cleaned);
61        if self.cleaned_window.len() > self.window_len {
62            self.cleaned_window.pop_back();
63        }
64
65        if self.cleaned_window.len() < self.window_len {
66            return 0.0;
67        }
68
69        // DFT
70        let mut pwr = vec![0.0; 51]; // Periods 8 to 50
71        for (period_idx, p) in pwr.iter_mut().enumerate().take(51).skip(8) {
72            let period = period_idx as f64;
73            let mut cos_part = 0.0;
74            let mut sin_part = 0.0;
75            for n in 0..self.window_len {
76                let angle = std::f64::consts::TAU * n as f64 / period;
77                cos_part += self.cleaned_window[n] * angle.cos();
78                sin_part += self.cleaned_window[n] * angle.sin();
79            }
80            *p = cos_part * cos_part + sin_part * sin_part;
81        }
82
83        // Max Power
84        let mut max_pwr = 0.0;
85        for &p in &pwr[8..=50] {
86            if p > max_pwr { max_pwr = p; }
87        }
88
89        // dB and CG
90        let mut num = 0.0;
91        let mut denom = 0.0;
92        for (period_idx, &p) in pwr.iter().enumerate().take(51).skip(8) {
93            let db = if max_pwr > 0.0 && p > 0.0 {
94                let val = 0.01 / (1.0 - 0.99 * p / max_pwr);
95                -10.0 * val.log10()
96            } else {
97                20.0
98            }.min(20.0);
99
100            if db < 3.0 {
101                num += period_idx as f64 * (3.0 - db);
102                denom += 3.0 - db;
103            }
104        }
105
106        if denom != 0.0 {
107            num / denom
108        } else {
109            0.0
110        }
111    }
112}
113
114pub const FOURIER_DOMINANT_CYCLE_METADATA: IndicatorMetadata = IndicatorMetadata {
115    name: "FourierDominantCycle",
116    description: "Dominant cycle period estimation using resolution-enhanced DFT and center of gravity.",
117    usage: "Use to compute the dominant market cycle period via DFT. Feed the output period into adaptive indicators like DSMA or Ehlers Stochastic to make them cycle-synchronized.",
118    keywords: &["cycle", "spectral", "ehlers", "dsp", "dominant-cycle", "fourier"],
119    ehlers_summary: "Ehlers implements a Discrete Fourier Transform cycle measurement in Cybernetic Analysis using a Hann-windowed data segment. The DFT computes power across periods from 6 to 50 bars, and the peak power identifies the dominant cycle period driving price movement.",
120    params: &[
121        ParamDef {
122            name: "window_len",
123            default: "50",
124            description: "DFT window length",
125        },
126    ],
127    formula_source: "https://github.com/lavs9/quantwave/blob/main/references/Ehlers%20Papers/FourierTransformForTraders.pdf",
128    formula_latex: r#"
129\[
130HP = \text{HighPass}(Price, 40)
131\]
132\[
133Cleaned = \frac{HP + 2HP_{t-1} + 3HP_{t-2} + 3HP_{t-3} + 2HP_{t-4} + HP_{t-5}}{12}
134\]
135\[
136Pwr(P) = \left(\sum_{n=0}^{W-1} Cleaned_{t-n} \cos\left(\frac{2\pi n}{P}\right)\right)^2 + \left(\sum_{n=0}^{W-1} Cleaned_{t-n} \sin\left(\frac{2\pi n}{P}\right)\right)^2
137\]
138\[
139DB(P) = \min\left(20, -10 \log_{10}\left(\frac{0.01}{1 - 0.99 \frac{Pwr(P)}{\max(Pwr)}}\right)\right)
140\]
141\[
142DC = \frac{\sum_{P=8}^{50} P \cdot (3 - DB(P)) \text{ where } DB(P) < 3}{\sum (3 - DB(P))}
143\]
144"#,
145    gold_standard_file: "fourier_dominant_cycle.json",
146    category: "Ehlers DSP",
147};
148
149#[cfg(test)]
150mod tests {
151    use std::f64::consts::PI;
152    use super::*;
153    use crate::traits::Next;
154    use proptest::prelude::*;
155
156    #[test]
157    fn test_fourier_dc_basic() {
158        let mut fdc = FourierDominantCycle::new(50);
159        for i in 0..200 {
160            // Sine wave with period 20
161            let val = fdc.next((2.0 * PI * i as f64 / 20.0).sin());
162            if i > 150 {
163                // Should be around 20
164                assert!(val > 15.0 && val < 25.0);
165            }
166        }
167    }
168
169    proptest! {
170        #[test]
171        fn test_fourier_dc_parity(
172            inputs in prop::collection::vec(1.0..100.0, 100..150),
173        ) {
174            let window_len = 50;
175            let mut fdc = FourierDominantCycle::new(window_len);
176            let streaming_results: Vec<f64> = inputs.iter().map(|&x| fdc.next(x)).collect();
177
178            // Batch implementation
179            let mut batch_results = Vec::with_capacity(inputs.len());
180            let mut hp = HighPass::new(40);
181            let hp_vals: Vec<f64> = inputs.iter().map(|&x| hp.next(x)).collect();
182            let mut cleaned_vals = Vec::new();
183
184            for i in 0..hp_vals.len() {
185                let c = (hp_vals[i] 
186                    + 2.0 * (if i > 0 { hp_vals[i-1] } else { 0.0 })
187                    + 3.0 * (if i > 1 { hp_vals[i-2] } else { 0.0 })
188                    + 3.0 * (if i > 2 { hp_vals[i-3] } else { 0.0 })
189                    + 2.0 * (if i > 3 { hp_vals[i-4] } else { 0.0 })
190                    + (if i > 4 { hp_vals[i-5] } else { 0.0 })) / 12.0;
191                cleaned_vals.push(c);
192
193                if i < window_len + 5 { // account for HP and FIR delay
194                    batch_results.push(0.0);
195                    continue;
196                }
197
198                let mut pwr = vec![0.0; 51];
199                for period_idx in 8..=50 {
200                    let period = period_idx as f64;
201                    let mut cos_part = 0.0;
202                    let mut sin_part = 0.0;
203                    for n in 0..window_len {
204                        let angle = 2.0 * PI * n as f64 / period;
205                        cos_part += cleaned_vals[i-n] * angle.cos();
206                        sin_part += cleaned_vals[i-n] * angle.sin();
207                    }
208                    pwr[period_idx] = cos_part * cos_part + sin_part * sin_part;
209                }
210
211                let mut max_p = 0.0;
212                for &p in &pwr[8..=50] {
213                    if p > max_p { max_p = p; }
214                }
215
216                let mut num = 0.0;
217                let mut den = 0.0;
218                for period_idx in 8..=50 {
219                    let p = pwr[period_idx];
220                    let db = if max_p > 0.0 && p > 0.0 {
221                        let val = 0.01 / (1.0 - 0.99 * p / max_p);
222                        -10.0 * val.log10()
223                    } else {
224                        20.0
225                    }.min(20.0);
226
227                    if db < 3.0 {
228                        num += period_idx as f64 * (3.0 - db);
229                        den += 3.0 - db;
230                    }
231                }
232
233                batch_results.push(if den != 0.0 { num / den } else { 0.0 });
234            }
235
236            // Note: streaming uses front-pushing window, batch uses indexing.
237            // There might be minor startup differences due to history initialization.
238            for i in 60..inputs.len() {
239                approx::assert_relative_eq!(streaming_results[i], batch_results[i], epsilon = 1e-10);
240            }
241        }
242    }
243}
quantwave_core/indicators/fourier_transform.rs

quantwave_core/indicators/
fourier_transform.rs
