sim_lib_numbers_signal/
reference.rs

1//! Definition-level O(N^2) reference transforms.
2//!
3//! These routines intentionally remain direct definitions. They provide the
4//! oracle for fast implementations and are not described as FFT algorithms.
5
6use std::f64::consts::{PI, SQRT_2, TAU};
7
8use crate::{
9    DctType, Direction, DstType, Normalization, SignConvention, SignalError, fft::Complex,
10};
11
12/// Computes the direct, unnormalized complex DFT definition.
13///
14/// The selected [`SignConvention`] and [`Direction`] determine the exponential
15/// sign. No scaling is applied.
16pub fn reference_dft(
17    input: &[(f64, f64)],
18    direction: Direction,
19    sign: SignConvention,
20) -> Result<Vec<(f64, f64)>, SignalError> {
21    if input.is_empty() {
22        return Err(SignalError::InvalidLength {
23            len: 0,
24            reason: "DFT requires at least one value",
25        });
26    }
27    validate_complex(input)?;
28    let sign = sign.angle_sign(direction);
29    let len = input.len() as f64;
30    let output = (0..input.len())
31        .map(|frequency| {
32            input
33                .iter()
34                .copied()
35                .enumerate()
36                .fold(Complex::ZERO, |sum, (sample, value)| {
37                    let angle = sign * TAU * frequency as f64 * sample as f64 / len;
38                    sum + Complex::from(value) * Complex::cis(angle)
39                })
40                .into()
41        })
42        .collect::<Vec<_>>();
43    validate_complex(&output)?;
44    Ok(output)
45}
46
47/// Computes DCT-I, II, III, or IV directly from its cosine definition.
48pub fn reference_dct(
49    input: &[f64],
50    kind: DctType,
51    direction: Direction,
52    normalization: Normalization,
53) -> Result<Vec<f64>, SignalError> {
54    validate_real(input)?;
55    if kind == DctType::I && input.len() < 2 {
56        return Err(SignalError::InvalidLength {
57            len: input.len(),
58            reason: "DCT-I requires at least two values",
59        });
60    }
61    if normalization == Normalization::Orthonormal {
62        let output = orthonormal_dct(input, kind, direction);
63        validate_real(&output)?;
64        return Ok(output);
65    }
66    let (kernel, factor) = match (kind, direction) {
67        (DctType::I, _) => (DctType::I, 2.0 * (input.len() - 1) as f64),
68        (DctType::II, Direction::Forward) | (DctType::III, Direction::Inverse) => {
69            (DctType::II, 2.0 * input.len() as f64)
70        }
71        (DctType::III, Direction::Forward) | (DctType::II, Direction::Inverse) => {
72            (DctType::III, 2.0 * input.len() as f64)
73        }
74        (DctType::IV, _) => (DctType::IV, 2.0 * input.len() as f64),
75    };
76    let mut output = raw_dct(input, kernel);
77    apply_pair_normalization(&mut output, factor, direction, normalization);
78    validate_real(&output)?;
79    Ok(output)
80}
81
82/// Computes DST-I, II, III, or IV directly from its sine definition.
83pub fn reference_dst(
84    input: &[f64],
85    kind: DstType,
86    direction: Direction,
87    normalization: Normalization,
88) -> Result<Vec<f64>, SignalError> {
89    validate_real(input)?;
90    if normalization == Normalization::Orthonormal {
91        let output = orthonormal_dst(input, kind, direction);
92        validate_real(&output)?;
93        return Ok(output);
94    }
95    let (kernel, factor) = match (kind, direction) {
96        (DstType::I, _) => (DstType::I, 2.0 * (input.len() + 1) as f64),
97        (DstType::II, Direction::Forward) | (DstType::III, Direction::Inverse) => {
98            (DstType::II, 2.0 * input.len() as f64)
99        }
100        (DstType::III, Direction::Forward) | (DstType::II, Direction::Inverse) => {
101            (DstType::III, 2.0 * input.len() as f64)
102        }
103        (DstType::IV, _) => (DstType::IV, 2.0 * input.len() as f64),
104    };
105    let mut output = raw_dst(input, kernel);
106    apply_pair_normalization(&mut output, factor, direction, normalization);
107    validate_real(&output)?;
108    Ok(output)
109}
110
111fn raw_dct(input: &[f64], kind: DctType) -> Vec<f64> {
112    let len = input.len();
113    match kind {
114        DctType::I => (0..len)
115            .map(|frequency| {
116                input[0]
117                    + if frequency % 2 == 0 {
118                        input[len - 1]
119                    } else {
120                        -input[len - 1]
121                    }
122                    + 2.0
123                        * input[1..len - 1]
124                            .iter()
125                            .enumerate()
126                            .map(|(offset, value)| {
127                                let sample = offset + 1;
128                                value
129                                    * (PI * frequency as f64 * sample as f64 / (len - 1) as f64)
130                                        .cos()
131                            })
132                            .sum::<f64>()
133            })
134            .collect(),
135        DctType::II => (0..len)
136            .map(|frequency| {
137                2.0 * input
138                    .iter()
139                    .enumerate()
140                    .map(|(sample, value)| {
141                        value
142                            * (PI * frequency as f64 * (2 * sample + 1) as f64 / (2 * len) as f64)
143                                .cos()
144                    })
145                    .sum::<f64>()
146            })
147            .collect(),
148        DctType::III => (0..len)
149            .map(|frequency| {
150                input[0]
151                    + 2.0
152                        * input[1..]
153                            .iter()
154                            .enumerate()
155                            .map(|(offset, value)| {
156                                let sample = offset + 1;
157                                value
158                                    * (PI * sample as f64 * (2 * frequency + 1) as f64
159                                        / (2 * len) as f64)
160                                        .cos()
161                            })
162                            .sum::<f64>()
163            })
164            .collect(),
165        DctType::IV => (0..len)
166            .map(|frequency| {
167                2.0 * input
168                    .iter()
169                    .enumerate()
170                    .map(|(sample, value)| {
171                        value
172                            * (PI * (2 * frequency + 1) as f64 * (2 * sample + 1) as f64
173                                / (4 * len) as f64)
174                                .cos()
175                    })
176                    .sum::<f64>()
177            })
178            .collect(),
179    }
180}
181
182fn raw_dst(input: &[f64], kind: DstType) -> Vec<f64> {
183    let len = input.len();
184    match kind {
185        DstType::I => (0..len)
186            .map(|frequency| {
187                2.0 * input
188                    .iter()
189                    .enumerate()
190                    .map(|(sample, value)| {
191                        value
192                            * (PI * (frequency + 1) as f64 * (sample + 1) as f64 / (len + 1) as f64)
193                                .sin()
194                    })
195                    .sum::<f64>()
196            })
197            .collect(),
198        DstType::II => (0..len)
199            .map(|frequency| {
200                2.0 * input
201                    .iter()
202                    .enumerate()
203                    .map(|(sample, value)| {
204                        value
205                            * (PI * (frequency + 1) as f64 * (2 * sample + 1) as f64
206                                / (2 * len) as f64)
207                                .sin()
208                    })
209                    .sum::<f64>()
210            })
211            .collect(),
212        DstType::III => (0..len)
213            .map(|frequency| {
214                let endpoint = if frequency % 2 == 0 {
215                    input[len - 1]
216                } else {
217                    -input[len - 1]
218                };
219                endpoint
220                    + 2.0
221                        * input[..len - 1]
222                            .iter()
223                            .enumerate()
224                            .map(|(sample, value)| {
225                                value
226                                    * (PI * (sample + 1) as f64 * (2 * frequency + 1) as f64
227                                        / (2 * len) as f64)
228                                        .sin()
229                            })
230                            .sum::<f64>()
231            })
232            .collect(),
233        DstType::IV => (0..len)
234            .map(|frequency| {
235                2.0 * input
236                    .iter()
237                    .enumerate()
238                    .map(|(sample, value)| {
239                        value
240                            * (PI * (2 * frequency + 1) as f64 * (2 * sample + 1) as f64
241                                / (4 * len) as f64)
242                                .sin()
243                    })
244                    .sum::<f64>()
245            })
246            .collect(),
247    }
248}
249
250fn orthonormal_dct(input: &[f64], kind: DctType, direction: Direction) -> Vec<f64> {
251    let len = input.len();
252    match kind {
253        DctType::I => matrix_apply(input, |row, column| {
254            let row_weight = endpoint_weight(row, len);
255            let column_weight = endpoint_weight(column, len);
256            (2.0 / (len - 1) as f64).sqrt()
257                * row_weight
258                * column_weight
259                * (PI * row as f64 * column as f64 / (len - 1) as f64).cos()
260        }),
261        DctType::II | DctType::III => {
262            let transpose = (kind == DctType::II) == (direction == Direction::Inverse);
263            matrix_apply(input, |row, column| {
264                let (frequency, sample) = if transpose {
265                    (column, row)
266                } else {
267                    (row, column)
268                };
269                (2.0 / len as f64).sqrt()
270                    * if frequency == 0 { 1.0 / SQRT_2 } else { 1.0 }
271                    * (PI * frequency as f64 * (2 * sample + 1) as f64 / (2 * len) as f64).cos()
272            })
273        }
274        DctType::IV => matrix_apply(input, |row, column| {
275            (2.0 / len as f64).sqrt()
276                * (PI * (2 * row + 1) as f64 * (2 * column + 1) as f64 / (4 * len) as f64).cos()
277        }),
278    }
279}
280
281fn orthonormal_dst(input: &[f64], kind: DstType, direction: Direction) -> Vec<f64> {
282    let len = input.len();
283    match kind {
284        DstType::I => matrix_apply(input, |row, column| {
285            (2.0 / (len + 1) as f64).sqrt()
286                * (PI * (row + 1) as f64 * (column + 1) as f64 / (len + 1) as f64).sin()
287        }),
288        DstType::II | DstType::III => {
289            let transpose = (kind == DstType::II) == (direction == Direction::Inverse);
290            matrix_apply(input, |row, column| {
291                let (frequency, sample) = if transpose {
292                    (column, row)
293                } else {
294                    (row, column)
295                };
296                (2.0 / len as f64).sqrt()
297                    * if frequency + 1 == len {
298                        1.0 / SQRT_2
299                    } else {
300                        1.0
301                    }
302                    * (PI * (frequency + 1) as f64 * (2 * sample + 1) as f64 / (2 * len) as f64)
303                        .sin()
304            })
305        }
306        DstType::IV => matrix_apply(input, |row, column| {
307            (2.0 / len as f64).sqrt()
308                * (PI * (2 * row + 1) as f64 * (2 * column + 1) as f64 / (4 * len) as f64).sin()
309        }),
310    }
311}
312
313fn endpoint_weight(index: usize, len: usize) -> f64 {
314    if index == 0 || index + 1 == len {
315        1.0 / SQRT_2
316    } else {
317        1.0
318    }
319}
320
321fn matrix_apply(input: &[f64], coefficient: impl Fn(usize, usize) -> f64) -> Vec<f64> {
322    (0..input.len())
323        .map(|row| {
324            input
325                .iter()
326                .enumerate()
327                .map(|(column, value)| coefficient(row, column) * value)
328                .sum()
329        })
330        .collect()
331}
332
333fn apply_pair_normalization(
334    output: &mut [f64],
335    factor: f64,
336    direction: Direction,
337    normalization: Normalization,
338) {
339    let scale = match (normalization, direction) {
340        (Normalization::None, _) => 1.0,
341        (Normalization::Forward, Direction::Forward)
342        | (Normalization::Inverse, Direction::Inverse) => 1.0 / factor,
343        (Normalization::Forward | Normalization::Inverse, _) => 1.0,
344        (Normalization::Orthonormal, _) => unreachable!("handled by orthonormal definitions"),
345    };
346    output.iter_mut().for_each(|value| *value *= scale);
347}
348
349fn validate_real(input: &[f64]) -> Result<(), SignalError> {
350    if input.is_empty() {
351        return Err(SignalError::InvalidLength {
352            len: 0,
353            reason: "real transforms require at least one value",
354        });
355    }
356    for (index, value) in input.iter().enumerate() {
357        if !value.is_finite() {
358            return Err(SignalError::NonFinite {
359                index,
360                component: "value",
361            });
362        }
363    }
364    Ok(())
365}
366
367fn validate_complex(input: &[(f64, f64)]) -> Result<(), SignalError> {
368    for (index, (real, imag)) in input.iter().copied().enumerate() {
369        if !real.is_finite() {
370            return Err(SignalError::NonFinite {
371                index,
372                component: "real",
373            });
374        }
375        if !imag.is_finite() {
376            return Err(SignalError::NonFinite {
377                index,
378                component: "imag",
379            });
380        }
381    }
382    Ok(())
383}
sim_lib_numbers_signal/reference.rs

sim_lib_numbers_signal/
reference.rs
