mfcc/
freqs.rs

1use std::f64;
2use std::sync::Arc;
3use num_complex::Complex64;
4#[cfg(feature = "fftrust")]
5use rustfft::{FFT, FFTplanner};
6#[cfg(feature = "fftrust")]
7use rustfft::num_complex::Complex;
8#[cfg(feature = "fftextern")]
9use fftw::plan::*;
10#[cfg(feature = "fftextern")]
11use fftw::types::*;
12
13#[cfg(feature = "fftrust")]
14pub struct InverseCosineTransform {
15    ifft: Arc<dyn FFT<f64>>,
16    buf: Vec<Complex64>,
17    buf2: Vec<Complex64>
18}
19
20#[cfg(feature = "fftrust")]
21impl InverseCosineTransform {
22    pub fn new(size: usize) -> InverseCosineTransform {
23        let mut planner = FFTplanner::new(true);
24        InverseCosineTransform {
25            ifft: planner.plan_fft(size),
26            buf: vec![Complex::i(); size],
27            buf2: vec![Complex::i(); size]
28        }
29    }
30
31    pub fn transform(&mut self, input: &[f64], output: &mut [f64]) {
32        let length = self.ifft.len();
33        let norm = 2.0 * length as f64;
34
35        for i in 0..length {
36            let theta = i as f64 / norm * f64::consts::PI;
37
38            self.buf[i] = Complex::from_polar(&(input[i] * norm.sqrt()), &theta);
39        }
40
41        self.buf[0] = self.buf[0].unscale(2.0_f64.sqrt());
42
43        //dbg!(&self.buf);
44        self.ifft.process(&mut self.buf, &mut self.buf2);
45
46        for i in 0..length/2 {
47            output[i*2] = self.buf2[i].re / (length as f64);
48            output[i*2+1] = self.buf2[length - i - 1].re  / (length as f64);
49        }
50    }
51}
52
53#[cfg(feature = "fftrust")]
54pub struct ForwardRealFourier {
55    fft: Arc<dyn FFT<f64>>,
56    buf: Vec<Complex64>,
57    buf2: Vec<Complex64>
58}
59
60#[cfg(feature = "fftrust")]
61impl ForwardRealFourier {
62    pub fn new(size: usize) -> ForwardRealFourier {
63        let mut planner = FFTplanner::new(false);
64        ForwardRealFourier {
65            fft: planner.plan_fft(size/2),
66            buf: vec![Complex::i(); size/2],
67            buf2: vec![Complex::i(); size/2+1]
68        }
69    }
70
71    pub fn transform(&mut self, input: &[f64], output: &mut [Complex64]) {
72        let length = self.fft.len();
73
74        for i in 0..length {
75            self.buf[i] = Complex::new(input[i*2], input[i*2 + 1]);
76        }
77
78        self.fft.process(&mut self.buf, &mut self.buf2[0..length]);
79
80        let first = self.buf2[0].clone();
81        self.buf2[length] = self.buf2[0];
82        for i in 0..length {
83            let cplx = Complex64::from_polar(&1.0, &(-f64::consts::PI / (length as f64) * (i as f64) + f64::consts::PI / 2.0));
84            output[i] = 0.5 * (self.buf2[i] + self.buf2[length - i].conj()) + 0.5 * cplx * (-self.buf2[i] + self.buf2[length - i].conj());
85        }
86
87        output[length] = Complex64::new(first.re - first.im, 0.0);
88    }
89}
90
91#[cfg(feature = "fftextern")]
92pub struct InverseCosineTransform {
93    dct_state: R2RPlan64
94}
95
96#[cfg(feature = "fftextern")]
97impl InverseCosineTransform {
98    pub fn new(size: usize) -> InverseCosineTransform {
99        InverseCosineTransform {
100            dct_state: R2RPlan::aligned(&[size], R2RKind::FFTW_REDFT01, Flag::Estimate).unwrap()
101        }
102    }
103
104    pub fn transform(&mut self, input: &mut [f64], output: &mut [f64]) {
105        input[0] *= 2.0f64.sqrt();
106
107        self.dct_state.r2r(input, output).unwrap();
108
109        for x in output {
110            *x = *x / (2.0 * input.len() as f64).sqrt();
111        }
112    }
113}
114
115#[cfg(feature = "fftextern")]
116pub struct ForwardRealFourier {
117    fft_state: R2CPlan64
118}
119
120#[cfg(feature = "fftextern")]
121impl ForwardRealFourier {
122    pub fn new(size: usize) -> ForwardRealFourier {
123        ForwardRealFourier {
124            fft_state: R2CPlan::aligned(&[size], Flag::Estimate).unwrap(),
125        }
126    }
127
128    pub fn transform(&mut self, input: &mut [f64], output: &mut [Complex64]) {
129        self.fft_state.r2c(input, output).unwrap();
130    }
131}
132
133