quantwave_core/indicators/
fourier_transform.rs

1use crate::indicators::high_pass::HighPass;
2use crate::indicators::metadata::{IndicatorMetadata, ParamDef};
3use crate::traits::Next;
4use crate::utils::RingBuffer as 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])
53            / 12.0;
54
55        // Shift HP history
56        for i in (1..6).rev() {
57            self.hp_history[i] = self.hp_history[i - 1];
58        }
59        self.hp_history[0] = hp_val;
60
61        self.cleaned_window.push_front(cleaned);
62        if self.cleaned_window.len() > self.window_len {
63            self.cleaned_window.pop_back();
64        }
65
66        if self.cleaned_window.len() < self.window_len {
67            return 0.0;
68        }
69
70        // DFT
71        let mut pwr = vec![0.0; 51]; // Periods 8 to 50
72        for (period_idx, p) in pwr.iter_mut().enumerate().take(51).skip(8) {
73            let period = period_idx as f64;
74            let mut cos_part = 0.0;
75            let mut sin_part = 0.0;
76            for n in 0..self.window_len {
77                let angle = std::f64::consts::TAU * n as f64 / period;
78                cos_part += self.cleaned_window[n] * angle.cos();
79                sin_part += self.cleaned_window[n] * angle.sin();
80            }
81            *p = cos_part * cos_part + sin_part * sin_part;
82        }
83
84        // Max Power
85        let mut max_pwr = 0.0;
86        for &p in &pwr[8..=50] {
87            if p > max_pwr {
88                max_pwr = p;
89            }
90        }
91
92        // dB and CG
93        let mut num = 0.0;
94        let mut denom = 0.0;
95        for (period_idx, &p) in pwr.iter().enumerate().take(51).skip(8) {
96            let db = if max_pwr > 0.0 && p > 0.0 {
97                let val = 0.01 / (1.0 - 0.99 * p / max_pwr);
98                -10.0 * val.log10()
99            } else {
100                20.0
101            }
102            .min(20.0);
103
104            if db < 3.0 {
105                num += period_idx as f64 * (3.0 - db);
106                denom += 3.0 - db;
107            }
108        }
109
110        if denom != 0.0 { num / denom } else { 0.0 }
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: &[
119        "cycle",
120        "spectral",
121        "ehlers",
122        "dsp",
123        "dominant-cycle",
124        "fourier",
125    ],
126    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.",
127    params: &[ParamDef {
128        name: "window_len",
129        default: "50",
130        description: "DFT window length",
131    }],
132    formula_source: "https://github.com/lavs9/quantwave/blob/main/references/Ehlers%20Papers/FourierTransformForTraders.pdf",
133    formula_latex: r#"
134\[
135HP = \text{HighPass}(Price, 40)
136\]
137\[
138Cleaned = \frac{HP + 2HP_{t-1} + 3HP_{t-2} + 3HP_{t-3} + 2HP_{t-4} + HP_{t-5}}{12}
139\]
140\[
141Pwr(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
142\]
143\[
144DB(P) = \min\left(20, -10 \log_{10}\left(\frac{0.01}{1 - 0.99 \frac{Pwr(P)}{\max(Pwr)}}\right)\right)
145\]
146\[
147DC = \frac{\sum_{P=8}^{50} P \cdot (3 - DB(P)) \text{ where } DB(P) < 3}{\sum (3 - DB(P))}
148\]
149"#,
150    gold_standard_file: "fourier_dominant_cycle.json",
151    category: "Ehlers DSP",
152};
153
154#[cfg(test)]
155mod tests {
156    use super::*;
157    use crate::traits::Next;
158    use proptest::prelude::*;
159    use std::f64::consts::PI;
160
161    #[test]
162    fn test_fourier_dc_basic() {
163        let mut fdc = FourierDominantCycle::new(50);
164        for i in 0..200 {
165            // Sine wave with period 20
166            let val = fdc.next((2.0 * PI * i as f64 / 20.0).sin());
167            if i > 150 {
168                // Should be around 20
169                assert!(val > 15.0 && val < 25.0);
170            }
171        }
172    }
173
174    proptest! {
175        #[test]
176        fn test_fourier_dc_parity(
177            inputs in prop::collection::vec(1.0..100.0, 100..150),
178        ) {
179            let window_len = 50;
180            let mut fdc = FourierDominantCycle::new(window_len);
181            let streaming_results: Vec<f64> = inputs.iter().map(|&x| fdc.next(x)).collect();
182
183            // Batch implementation
184            let mut batch_results = Vec::with_capacity(inputs.len());
185            let mut hp = HighPass::new(40);
186            let hp_vals: Vec<f64> = inputs.iter().map(|&x| hp.next(x)).collect();
187            let mut cleaned_vals = Vec::new();
188
189            for i in 0..hp_vals.len() {
190                let c = (hp_vals[i]
191                    + 2.0 * (if i > 0 { hp_vals[i-1] } else { 0.0 })
192                    + 3.0 * (if i > 1 { hp_vals[i-2] } else { 0.0 })
193                    + 3.0 * (if i > 2 { hp_vals[i-3] } else { 0.0 })
194                    + 2.0 * (if i > 3 { hp_vals[i-4] } else { 0.0 })
195                    + (if i > 4 { hp_vals[i-5] } else { 0.0 })) / 12.0;
196                cleaned_vals.push(c);
197
198                if i < window_len + 5 { // account for HP and FIR delay
199                    batch_results.push(0.0);
200                    continue;
201                }
202
203                let mut pwr = vec![0.0; 51];
204                for period_idx in 8..=50 {
205                    let period = period_idx as f64;
206                    let mut cos_part = 0.0;
207                    let mut sin_part = 0.0;
208                    for n in 0..window_len {
209                        let angle = 2.0 * PI * n as f64 / period;
210                        cos_part += cleaned_vals[i-n] * angle.cos();
211                        sin_part += cleaned_vals[i-n] * angle.sin();
212                    }
213                    pwr[period_idx] = cos_part * cos_part + sin_part * sin_part;
214                }
215
216                let mut max_p = 0.0;
217                for &p in &pwr[8..=50] {
218                    if p > max_p { max_p = p; }
219                }
220
221                let mut num = 0.0;
222                let mut den = 0.0;
223                for period_idx in 8..=50 {
224                    let p = pwr[period_idx];
225                    let db = if max_p > 0.0 && p > 0.0 {
226                        let val = 0.01 / (1.0 - 0.99 * p / max_p);
227                        -10.0 * val.log10()
228                    } else {
229                        20.0
230                    }.min(20.0);
231
232                    if db < 3.0 {
233                        num += period_idx as f64 * (3.0 - db);
234                        den += 3.0 - db;
235                    }
236                }
237
238                batch_results.push(if den != 0.0 { num / den } else { 0.0 });
239            }
240
241            // Note: streaming uses front-pushing window, batch uses indexing.
242            // There might be minor startup differences due to history initialization.
243            for i in 60..inputs.len() {
244                approx::assert_relative_eq!(streaming_results[i], batch_results[i], epsilon = 1e-10);
245            }
246        }
247    }
248}
quantwave_core/indicators/fourier_transform.rs

quantwave_core/indicators/
fourier_transform.rs
