Skip to main content

hilbert_transform/
lib.rs

1use rustfft::{FftNum, FftPlanner, num_complex::Complex};
2
3/// Hilbert_transform is a library written in Rust to perform the hilbert transformation like
4/// Matlab/Octave or scipy.signals.hilbert.
5///
6/// Hilbert_transform is implemented based on scipy implementation of same function.
7///
8/// # Usage
9///
10/// ```
11/// use hilbert_transform::{hilbert};
12///
13/// let input = vec![1.0, 2.0, 3.0, 4.0];     
14/// let hilbert_output = hilbert(&input);
15/// ```
16///
17/// hilbert_output will be equal to: [Complex { re: 1.0, im: 1.0 }, Complex { re: 2.0, im: -1.0 }, Complex { re: 3.0, im: -1.0 }, Complex { re: 4.0, im: 1.0 }]
18///
19
20pub fn hilbert<T: FftNum>(input: &[T]) -> Vec<Complex<T>> {
21
22    let len = input.len();
23    let mut planner = FftPlanner::<T>::new();
24    let fft = planner.plan_fft_forward(len);
25
26    let mut fft_complex = input
27        .iter()
28        .map(|&val| Complex::new(val, T::zero()))
29        .collect::<Vec<_>>();
30
31    fft.process(&mut fft_complex);
32
33    let mut h_spectrum = vec![Complex::new(T::zero(), T::zero()); len];
34
35    let two = T::from_usize(2).unwrap();
36
37    if len % 2 == 0 {
38        h_spectrum[0] = Complex::new(T::one(), T::zero());
39        h_spectrum[len / 2] = Complex::new(T::one(), T::zero());
40        for i in 1..(len / 2) {
41            h_spectrum[i] = Complex::new(two, T::zero());
42        }
43    } else {
44        h_spectrum[0] = Complex::new(T::one(), T::zero());
45        for i in 1..((len + 1) / 2) {
46        h_spectrum[i] = Complex::new(two, T::zero());
47        }
48    }
49
50    for i in 0..len {
51        fft_complex[i] = fft_complex[i] * h_spectrum[i];
52    }
53
54    let mut ifft_complex = fft_complex.clone();
55    let ifft = planner.plan_fft_inverse(len);
56    ifft.process(&mut ifft_complex);
57
58    let len = T::from_usize(len).expect("failed to convert len to T");
59    for val in ifft_complex.iter_mut() {
60        *val = *val / len;
61    }
62
63    ifft_complex
64}
65
66#[cfg(test)]
67mod tests {
68    use super::*;
69
70    #[test]
71    fn test_hilbert_f64() {
72        let input = vec![1.0, 2.0, 3.0, 4.0];
73        let expected_output = vec![
74            Complex::new(1.0, 1.0), 
75            Complex::new(2.0, -1.0), 
76            Complex::new(3.0, -1.0), 
77            Complex::new(4.0, 1.0)
78        ];
79        assert_eq!(hilbert(&input), expected_output);
80    }
81
82
83    #[test]
84    fn test_hilbert_f32() {
85        let input: Vec<f32> = vec![1.0, 2.0, 3.0, 4.0];
86        let expected_output = vec![
87            Complex::new(1.0, 1.0), 
88            Complex::new(2.0, -1.0), 
89            Complex::new(3.0, -1.0), 
90            Complex::new(4.0, 1.0)
91        ];
92        assert_eq!(hilbert(&input), expected_output);
93    }
94}