hilbert_transform 0.1.2

An implementation of Hilbert Transformation like Matlab/Octave hilbert function
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

use rustfft::{FftNum, FftPlanner, num_complex::Complex};

/// Hilbert_transform is a library written in Rust to perform the hilbert transformation like
/// Matlab/Octave or scipy.signals.hilbert.
///
/// Hilbert_transform is implemented based on scipy implementation of same function.
///
/// # Usage
///
/// ```
/// use hilbert_transform::{hilbert};
///
/// let input = vec![1.0, 2.0, 3.0, 4.0];     
/// let hilbert_output = hilbert(&input);
/// ```
///
/// 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 }]
///

pub fn hilbert<T: FftNum>(input: &[T]) -> Vec<Complex<T>> {

    let len = input.len();
    let mut planner = FftPlanner::<T>::new();
    let fft = planner.plan_fft_forward(len);

    let mut fft_complex = input
        .iter()
        .map(|&val| Complex::new(val, T::zero()))
        .collect::<Vec<_>>();

    fft.process(&mut fft_complex);

    let mut h_spectrum = vec![Complex::new(T::zero(), T::zero()); len];

    let two = T::from_usize(2).unwrap();

    if len % 2 == 0 {
        h_spectrum[0] = Complex::new(T::one(), T::zero());
        h_spectrum[len / 2] = Complex::new(T::one(), T::zero());
        for i in 1..(len / 2) {
            h_spectrum[i] = Complex::new(two, T::zero());
        }
    } else {
        h_spectrum[0] = Complex::new(T::one(), T::zero());
        for i in 1..((len + 1) / 2) {
        h_spectrum[i] = Complex::new(two, T::zero());
        }
    }

    for i in 0..len {
        fft_complex[i] = fft_complex[i] * h_spectrum[i];
    }

    let mut ifft_complex = fft_complex.clone();
    let ifft = planner.plan_fft_inverse(len);
    ifft.process(&mut ifft_complex);

    let len = T::from_usize(len).expect("failed to convert len to T");
    for val in ifft_complex.iter_mut() {
        *val = *val / len;
    }

    ifft_complex
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_hilbert_f64() {
        let input = vec![1.0, 2.0, 3.0, 4.0];
        let expected_output = vec![
            Complex::new(1.0, 1.0), 
            Complex::new(2.0, -1.0), 
            Complex::new(3.0, -1.0), 
            Complex::new(4.0, 1.0)
        ];
        assert_eq!(hilbert(&input), expected_output);
    }


    #[test]
    fn test_hilbert_f32() {
        let input: Vec<f32> = vec![1.0, 2.0, 3.0, 4.0];
        let expected_output = vec![
            Complex::new(1.0, 1.0), 
            Complex::new(2.0, -1.0), 
            Complex::new(3.0, -1.0), 
            Complex::new(4.0, 1.0)
        ];
        assert_eq!(hilbert(&input), expected_output);
    }
}