fixed_dsp/transform/
mfcc.rs

1use crate::basic::{
2    div_i16, div_i32, dot_i16, dot_i32, offset_i32, scale_i16, scale_i32, shift_i32,
3    vlog_i32_in_place,
4};
5use crate::complex::{cmplx_mag_i16, cmplx_mag_i32};
6use crate::matrix::{Matrix, mat_vec_mul_i16, mat_vec_mul_i32};
7use crate::statistics::{absmax_i16, absmax_i32};
8use crate::transform::{RfftI16, RfftI32};
9use crate::{sat_i16, sat_i32};
10
11const LOG2TOLOG_Q15: i32 = 0x02C5_C860;
12const MICRO_Q15: i64 = 0x0000_0219;
13const SHIFT_MELFILTER_SATURATION_Q15: i32 = 10;
14
15const LOG2TOLOG_Q31: i32 = 0x02C5_C860;
16const MICRO_Q31: i64 = 0x0863_7BD0;
17const SHIFT_MELFILTER_SATURATION_Q31: i32 = 10;
18
19#[inline]
20fn view_u16_as_i16(data: &[u16]) -> &[i16] {
21    unsafe { core::slice::from_raw_parts(data.as_ptr() as *const i16, data.len()) }
22}
23
24#[inline]
25fn view_u32_as_i32(data: &[u32]) -> &[i32] {
26    unsafe { core::slice::from_raw_parts(data.as_ptr() as *const i32, data.len()) }
27}
28
29#[inline]
30fn view_i16_as_i32_mut(data: &mut [i16], min_i32_len: usize) -> &mut [i32] {
31    let (_, body, _) = unsafe { data.align_to_mut::<i32>() };
32    assert!(
33        body.len() >= min_i32_len,
34        "tmp alignment/length must allow at least min_i32_len i32 values"
35    );
36    body
37}
38
39pub struct MfccI16 {
40    pub n_fft: usize,
41    pub n_mels: usize,
42    pub n_mfcc: usize,
43    pub dct: &'static [u16],
44    pub filter: &'static [u16],
45    pub filter_pos: &'static [u16],
46    pub filter_len: &'static [u16],
47    pub window: &'static [u16],
48    pub rfft: RfftI16,
49}
50
51impl MfccI16 {
52    pub const fn new(
53        n_fft: usize,               // FFT size
54        n_mels: usize,              // Number of Mel filters
55        n_mfcc: usize,              // Number of MFCC coefficients
56        dct: &'static [u16],        // DCT matrix for MFCC computation
57        filter: &'static [u16],     // Mel filter bank coefficients
58        filter_pos: &'static [u16], // Starting positions of each filter in the FFT output
59        filter_len: &'static [u16], // Lengths of each filter
60        window: &'static [u16],     // Window function coefficients
61    ) -> Self {
62        assert!(
63            dct.len() == n_mfcc * n_mels,
64            "DCT matrix size must be n_mfcc * n_mels"
65        );
66        assert!(
67            filter_pos.len() == n_mels,
68            "Filter position array size must be n_mels"
69        );
70        assert!(
71            filter_len.len() == n_mels,
72            "Filter length array size must be n_mels"
73        );
74        assert!(window.len() == n_fft, "Window size must be equal to n_fft");
75
76        Self {
77            n_fft,
78            n_mels,
79            n_mfcc,
80            dct,
81            filter,
82            filter_pos,
83            filter_len,
84            window,
85            rfft: RfftI16::new(n_fft, false, true),
86        }
87    }
88
89    pub fn run(&self, input: &mut [i16], output: &mut [i16], tmp: &mut [i16]) {
90        assert!(
91            tmp.len() >= self.n_fft * 2,
92            "tmp length must be at least 2 * n_fft for q15 MFCC"
93        );
94
95        let dct = view_u16_as_i16(self.dct);
96        let filter = view_u16_as_i16(self.filter);
97
98        // 1) Front-end normalization to reduce overflow risk in fixed-point stages.
99        let (max_abs, _) = absmax_i16(input, self.n_fft);
100
101        if max_abs != 0 && max_abs != i16::MAX {
102            let (quotient, shift) = div_i16(i16::MAX, max_abs).expect("division must succeed");
103            scale_i16(input, quotient, shift as i8);
104        }
105
106        for (sample, &window) in input.iter_mut().zip(self.window) {
107            *sample = sat_i16(((*sample as i32) * (window as i32)) >> 15);
108        }
109
110        // 2) Time-domain to frequency-domain: RFFT then magnitude spectrum.
111        {
112            let tmp_fft = &mut tmp[..self.n_fft * 2];
113            self.rfft.run(input, tmp_fft);
114
115            let filter_limit = 1 + (self.n_fft >> 1);
116            cmplx_mag_i16(&tmp_fft[..filter_limit * 2], &mut input[..filter_limit]);
117        }
118
119        // 3) Apply packed Mel filterbank (dot products over magnitude bins).
120        let tmp_q31 = view_i16_as_i32_mut(tmp, self.n_mels);
121
122        let mut packed_pos = 0usize;
123        for ((dst, &filter_pos), &filter_len) in tmp_q31
124            .iter_mut()
125            .zip(self.filter_pos)
126            .zip(self.filter_len)
127            .take(self.n_mels)
128        {
129            let filter_pos = usize::from(filter_pos);
130            let filter_len = usize::from(filter_len);
131            let filter_end = packed_pos + filter_len;
132
133            let acc = dot_i16(
134                &input[filter_pos..filter_pos + filter_len],
135                &filter[packed_pos..filter_end],
136            );
137            packed_pos = filter_end;
138
139            let acc = (acc + MICRO_Q15) >> SHIFT_MELFILTER_SATURATION_Q15;
140            *dst = sat_i32(acc);
141        }
142
143        // 4) Log-energy domain conversion with FFT/saturation compensation.
144        if max_abs != 0 && max_abs != i16::MAX {
145            scale_i32(&mut tmp_q31[..self.n_mels], (max_abs as i32) << 16, 0);
146        }
147
148        vlog_i32_in_place(&mut tmp_q31[..self.n_mels]);
149
150        let fft_shift = self.n_fft.trailing_zeros() as i32;
151        let log_exponent =
152            (fft_shift + 2 + SHIFT_MELFILTER_SATURATION_Q15).wrapping_mul(LOG2TOLOG_Q15);
153        offset_i32(&mut tmp_q31[..self.n_mels], log_exponent);
154        shift_i32(&mut tmp_q31[..self.n_mels], -19);
155
156        for (dst, &src) in input.iter_mut().zip(tmp_q31.iter()).take(self.n_mels) {
157            *dst = sat_i16(src);
158        }
159
160        // 5) DCT projection from Mel log energies to MFCC coefficients.
161        let dct = Matrix {
162            rows: self.n_mfcc,
163            cols: self.n_mels,
164            data: dct.as_ptr() as *mut i16,
165        };
166        mat_vec_mul_i16(dct, &input[..self.n_mels], output);
167    }
168}
169
170pub struct MfccI32 {
171    pub n_fft: usize,
172    pub n_mels: usize,
173    pub n_mfcc: usize,
174    pub dct: &'static [u32],
175    pub filter: &'static [u32],
176    pub filter_pos: &'static [u16],
177    pub filter_len: &'static [u16],
178    pub window: &'static [u32],
179    pub rfft: RfftI32,
180}
181
182impl MfccI32 {
183    pub const fn new(
184        n_fft: usize,
185        n_mels: usize,
186        n_mfcc: usize,
187        dct: &'static [u32],
188        filter: &'static [u32],
189        filter_pos: &'static [u16],
190        filter_len: &'static [u16],
191        window: &'static [u32],
192    ) -> Self {
193        assert!(
194            dct.len() == n_mfcc * n_mels,
195            "DCT matrix size must be n_mfcc * n_mels"
196        );
197        assert!(
198            filter_pos.len() == n_mels,
199            "Filter position array size must be n_mels"
200        );
201        assert!(
202            filter_len.len() == n_mels,
203            "Filter length array size must be n_mels"
204        );
205        assert!(window.len() == n_fft, "Window size must be equal to n_fft");
206
207        Self {
208            n_fft,
209            n_mels,
210            n_mfcc,
211            dct,
212            filter,
213            filter_pos,
214            filter_len,
215            window,
216            rfft: RfftI32::new(n_fft, false, true),
217        }
218    }
219
220    pub fn run(&self, input: &mut [i32], output: &mut [i32], tmp: &mut [i32]) {
221        let dct = view_u32_as_i32(self.dct);
222        let filter = view_u32_as_i32(self.filter);
223
224        assert!(
225            tmp.len() >= self.n_fft * 2,
226            "tmp length must be at least 2 * n_fft for q31 MFCC"
227        );
228
229        // 1) Front-end normalization to reduce overflow risk in fixed-point stages.
230        let (max_abs, _) = absmax_i32(input, self.n_fft);
231
232        if max_abs != 0 && max_abs != i32::MAX {
233            let (quotient, shift) = div_i32(i32::MAX, max_abs).expect("division must succeed");
234            scale_i32(input, quotient, shift as i8);
235        }
236
237        for (sample, &window) in input.iter_mut().zip(self.window) {
238            *sample = sat_i32((((*sample as i64) * (window as i64)) >> 32) << 1);
239        }
240
241        // 2) Time-domain to frequency-domain: RFFT then magnitude spectrum.
242        self.rfft.run(input, &mut tmp[..self.n_fft * 2]);
243
244        let filter_limit = 1 + (self.n_fft >> 1);
245        cmplx_mag_i32(&tmp[..filter_limit * 2], &mut input[..filter_limit]);
246
247        // 3) Apply packed Mel filterbank (dot products over magnitude bins).
248        let mut packed_pos = 0usize;
249        for ((dst, &filter_pos), &filter_len) in tmp
250            .iter_mut()
251            .zip(self.filter_pos)
252            .zip(self.filter_len)
253            .take(self.n_mels)
254        {
255            let filter_pos = usize::from(filter_pos);
256            let filter_len = usize::from(filter_len);
257            let filter_end = packed_pos + filter_len;
258
259            let mut acc = dot_i32(
260                &input[filter_pos..filter_pos + filter_len],
261                &filter[packed_pos..filter_end],
262            );
263            packed_pos = filter_end;
264
265            acc = (acc + MICRO_Q31) >> (SHIFT_MELFILTER_SATURATION_Q31 + 18);
266            *dst = sat_i32(acc);
267        }
268
269        // 4) Log-energy domain conversion with FFT/saturation compensation.
270        if max_abs != 0 && max_abs != i32::MAX {
271            scale_i32(&mut tmp[..self.n_mels], max_abs, 0);
272        }
273
274        vlog_i32_in_place(&mut tmp[..self.n_mels]);
275
276        let fft_shift = 31 - (self.n_fft as u32).leading_zeros() as i32;
277        let log_exponent =
278            (fft_shift + 2 + SHIFT_MELFILTER_SATURATION_Q31).wrapping_mul(LOG2TOLOG_Q31);
279        offset_i32(&mut tmp[..self.n_mels], log_exponent);
280        shift_i32(&mut tmp[..self.n_mels], -3);
281
282        // 5) DCT projection from Mel log energies to MFCC coefficients.
283        let dct = Matrix {
284            rows: self.n_mfcc,
285            cols: self.n_mels,
286            data: dct.as_ptr() as *mut i32,
287        };
288        mat_vec_mul_i32(dct, &tmp[..self.n_mels], output);
289    }
290}
fixed_dsp/transform/mfcc.rs

fixed_dsp/transform/
mfcc.rs
