symphonia_core/dsp/mdct/
no_simd.rs

1// Symphonia
2// Copyright (c) 2019-2023 The Project Symphonia Developers.
3//
4// This Source Code Form is subject to the terms of the Mozilla Public
5// License, v. 2.0. If a copy of the MPL was not distributed with this
6// file, You can obtain one at https://mozilla.org/MPL/2.0/.
7
8//! The Modified Discrete Cosine Transform (MDCT) implemented without SIMD optimizations.
9
10use crate::dsp::complex::Complex;
11use crate::dsp::fft::*;
12
13/// The Inverse Modified Discrete Transform (IMDCT).
14pub struct Imdct {
15    fft: Fft,
16    fft_in: Box<[Complex]>,
17    fft_out: Box<[Complex]>,
18    twiddle: Box<[Complex]>,
19}
20
21impl Imdct {
22    /// Instantiate a N-point IMDCT with no scaling.
23    ///
24    /// The value of `n` is the number of spectral samples and must be a power-of-2 and less-than or
25    /// equal to `2 * Fft::MAX_SIZE`.
26    pub fn new(n: usize) -> Self {
27        Imdct::new_scaled(n, 1.0)
28    }
29
30    /// Instantiate a N-point IMDCT with scaling.
31    ///
32    /// The value of `n` is the number of spectral samples and must be a power-of-2 and less-than or
33    /// equal to `2 * Fft::MAX_SIZE`.
34    pub fn new_scaled(n: usize, scale: f64) -> Self {
35        // The FFT requires a power-of-two N.
36        assert!(n.is_power_of_two(), "n must be a power of two");
37        // A complex FFT of size N/2 is used to compute the IMDCT. Therefore, the maximum value of N
38        // is 2 * Fft::MAX_SIZE.
39        assert!(n <= 2 * Fft::MAX_SIZE, "maximum size exceeded");
40
41        let n2 = n / 2;
42        let mut twiddle = Vec::with_capacity(n2);
43
44        let alpha = 1.0 / 8.0 + if scale.is_sign_positive() { 0.0 } else { n2 as f64 };
45        let pi_n = std::f64::consts::PI / n as f64;
46        let sqrt_scale = scale.abs().sqrt();
47
48        for k in 0..n2 {
49            let theta = pi_n * (alpha + k as f64);
50            let re = sqrt_scale * theta.cos();
51            let im = sqrt_scale * theta.sin();
52            twiddle.push(Complex::new(re as f32, im as f32));
53        }
54
55        let fft_in = vec![Default::default(); n2].into_boxed_slice();
56        let fft_out = vec![Default::default(); n2].into_boxed_slice();
57
58        Imdct { fft: Fft::new(n2), fft_in, fft_out, twiddle: twiddle.into_boxed_slice() }
59    }
60
61    /// Performs the the N-point Inverse Modified Discrete Cosine Transform.
62    ///
63    /// The number of input spectral samples provided by the slice `spec` must equal the value of N
64    /// that the IMDCT was instantiated with. The length of the output slice, `out`, must be of
65    /// length 2N. Failing to meet these requirements will throw an assertion.
66    pub fn imdct(&mut self, spec: &[f32], out: &mut [f32]) {
67        // Spectral length: 2x FFT size, 0.5x output length.
68        let n = self.fft.size() << 1;
69        // 1x FFT size, 0.25x output length.
70        let n2 = n >> 1;
71        // 0.5x FFT size.
72        let n4 = n >> 2;
73
74        // The spectrum length must be the same as N.
75        assert_eq!(spec.len(), n);
76        // The output length must be 2x the spectrum length.
77        assert_eq!(out.len(), 2 * n);
78
79        // Pre-FFT twiddling and packing of the real input signal values into complex signal values.
80        for (i, (&w, t)) in self.twiddle.iter().zip(self.fft_in.iter_mut()).enumerate() {
81            let even = spec[i * 2];
82            let odd = -spec[n - 1 - i * 2];
83
84            let re = odd * w.im - even * w.re;
85            let im = odd * w.re + even * w.im;
86            *t = Complex::new(re, im);
87        }
88
89        // Do the FFT.
90        self.fft.fft(&self.fft_in, &mut self.fft_out);
91
92        // Split the output vector (2N samples) into 4 vectors (N/2 samples each).
93        let (vec0, vec1) = out.split_at_mut(n2);
94        let (vec1, vec2) = vec1.split_at_mut(n2);
95        let (vec2, vec3) = vec2.split_at_mut(n2);
96
97        // Post-FFT twiddling and processing to expand the N/2 complex output values into 2N real
98        // output samples.
99        for (i, (x, &w)) in self.fft_out[..n4].iter().zip(self.twiddle[..n4].iter()).enumerate() {
100            // The real and imaginary components of the post-twiddled FFT samples are used to
101            // generate 4 reak output samples. Using the first half of the complex FFT output,
102            // populate each of the 4 output vectors.
103            let val = w * x.conj();
104
105            // Forward and reverse order indicies that will be populated.
106            let fi = 2 * i;
107            let ri = n2 - 1 - 2 * i;
108
109            // The odd indicies in vec0 are populated reverse order.
110            vec0[ri] = -val.im;
111            // The even indicies in vec1 are populated forward order.
112            vec1[fi] = val.im;
113            // The odd indicies in vec2 are populated reverse order.
114            vec2[ri] = val.re;
115            // The even indicies in vec3 are populated forward order.
116            vec3[fi] = val.re;
117        }
118
119        for (i, (x, &w)) in self.fft_out[n4..].iter().zip(self.twiddle[n4..].iter()).enumerate() {
120            // Using the second half of the FFT output samples, finish populating each of the 4
121            // output vectors.
122            let val = w * x.conj();
123
124            // Forward and reverse order indicies that will be populated.
125            let fi = 2 * i;
126            let ri = n2 - 1 - 2 * i;
127
128            // The even indicies in vec0 are populated in forward order.
129            vec0[fi] = -val.re;
130            // The odd indicies in vec1 are populated in reverse order.
131            vec1[ri] = val.re;
132            // The even indicies in vec2 are populated in forward order.
133            vec2[fi] = val.im;
134            // The odd indicies in vec3 are populated in reverse order.
135            vec3[ri] = val.im;
136        }
137
138        // Note: As of Rust 1.58, there doesn't appear to be any measurable difference between using
139        // iterators or indexing like above. Either the bounds checks are elided above, or they are
140        // not elided using iterators. Therefore, for clarity, the indexing method is used.
141        //
142        // Additionally, note that vectors 0 and 3 are reversed copies (+ negation for vector 0) of
143        // vectors 1 and 2, respectively. Pulling these copies out into a separate loop using
144        // iterators yielded no measureable difference either.
145    }
146}