sci_rs/signal/
resample.rs

1use nalgebra::Complex;
2use num_traits::{Float, Zero};
3use rustfft::FftNum;
4
5///
6/// Resample the data to the desired number of samples using the Fourier transform.
7///
8/// This method is similar but not exactly equivalent to the SciPy method of resampling:
9/// <https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.resample.html>
10///
11/// It skips some complexity of the SciPy method, such as windowing and handling odd vs. even-length signals.
12///
13/// Procedure:
14/// 1. Convert to the frequency-domain.
15///   a. If upsampling, pad higher frequency bins with 0
16///   b. If downsampling, truncate higher frequency bins
17/// 2. Convert back to the time-domain.
18///
19pub fn resample<F: Float + FftNum>(x: &[F], n: usize) -> Vec<F> {
20    // SciPy style 'Fourier' resampling
21    // 1. Compute FFT of x
22    // 2. Fill vec of zeros with the desired length, y.
23    // 3. Set the from beginning of y to the first half of x
24    // 4. Set the from end of y to the second half of x
25    // 5. Compute IFFT of y
26    // 6. Multiply y by (n / x.len())
27    // 7. Take the real part of y
28
29    // Compute FFT of x
30    let mut fft_planner = rustfft::FftPlanner::<F>::new();
31    let fft = fft_planner.plan_fft_forward(x.len());
32    let ifft = fft_planner.plan_fft_inverse(n);
33
34    let scratch_length = std::cmp::max(
35        fft.get_inplace_scratch_len(),
36        ifft.get_inplace_scratch_len(),
37    );
38    let mut scratch = vec![Complex::zero(); scratch_length];
39    let mut x = x
40        .into_iter()
41        .map(|x| Complex::new(*x, F::zero()))
42        .collect::<Vec<_>>();
43    fft.process_with_scratch(&mut x, &mut scratch);
44
45    // Fill y with halfs of x
46    let mut y = vec![Complex::zero(); n];
47    let bins = std::cmp::min(x.len(), n);
48    let half_spectrum = bins / 2;
49    y[..half_spectrum].copy_from_slice(&x[..half_spectrum]);
50    y[n - half_spectrum..].copy_from_slice(&x[x.len() - half_spectrum..]);
51
52    // Compute iFFT of y
53    ifft.process_with_scratch(&mut y, &mut scratch);
54
55    // Take the scaled real domain as the resampled result
56    let scale_factor = F::from(1.0 / x.len() as f64).unwrap();
57    let y = y.iter().map(|x| x.re * scale_factor).collect::<Vec<_>>();
58
59    y
60}
61
62#[cfg(test)]
63mod tests {
64    use approx::assert_relative_eq;
65    use rand::Rng;
66
67    use super::*;
68
69    #[test]
70    #[should_panic]
71    fn can_resample_like_scipy() {
72        let x = vec![1., 2., 3., 4., 5., 6., 7., 8., 9.];
73        let y = resample(&x, 5);
74        let expected = vec![3., 2.18649851, 5.01849831, 5.98150169, 8.81350149];
75        assert_eq!(y.len(), expected.len());
76
77        // Fails due to slight algorithmic differences
78        println!("y: {:?}, scipy: {:?}", y, expected);
79        for (y, expected) in y.iter().zip(expected.iter()) {
80            assert_relative_eq!(y, expected, epsilon = 1e-10);
81        }
82    }
83
84    #[test]
85    fn can_resample_to_exact_number() {
86        // Randomly generate 1000 vectors of length (10,1000)
87        // Resample each to length 100
88        // Check that each resampled vector is length 100
89
90        let mut rng = rand::thread_rng();
91        for i in (0..100) {
92            let len = rng.gen_range((10..50));
93            let x = (0..len)
94                .map(|_| rng.gen_range((-100.0..100.)))
95                .collect::<Vec<_>>();
96            let y = resample(&x, 100);
97            assert_eq!(y.len(), 100);
98        }
99
100        for i in (0..50) {
101            let len = rng.gen_range((200..10000));
102            let target_len = rng.gen_range((50..50000));
103            let x = (0..len)
104                .map(|_| rng.gen_range((-100.0..100.)))
105                .collect::<Vec<_>>();
106            let y = resample(&x, target_len);
107            assert_eq!(y.len(), target_len);
108        }
109    }
110}
sci_rs/signal/resample.rs

sci_rs/signal/
resample.rs
