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}