Skip to main content

rustradio/
hilbert.rs

1/*! Hilbert transform.
2
3[Wikipedia][wiki] has a bunch of math, but one use case for it is to
4convert floating point values (think audio waveform) into upper
5sideband.
6
7Then again I guess you can do the same with a `FloatToComplex` plus
8`FftFilter`.
9
10This implementation is a pretty inefficient.
11
12[wiki]: https://en.wikipedia.org/wiki/Hilbert_transform
13*/
14
15use crate::block::{Block, BlockRet};
16use crate::fir::Fir;
17use crate::stream::{ReadStream, WriteStream};
18use crate::window::WindowType;
19use crate::{Complex, Float, Result};
20
21/// Hilbert transformer block.
22#[derive(rustradio_macros::Block)]
23#[rustradio(crate)]
24pub struct Hilbert {
25    #[rustradio(in)]
26    src: ReadStream<Float>,
27    #[rustradio(out)]
28    dst: WriteStream<Complex>,
29    history: Vec<Float>,
30    filter: Fir<Float>,
31    ntaps: usize,
32}
33
34impl Hilbert {
35    /// Create new hilber transformer with this many taps.
36    #[must_use]
37    pub fn new(
38        src: ReadStream<Float>,
39        ntaps: usize,
40        window_type: &WindowType,
41    ) -> (Self, ReadStream<Complex>) {
42        // TODO: take window function.
43        assert!(ntaps & 1 == 1, "hilbert filter len must be odd");
44        let taps = crate::fir::hilbert(&window_type.make_window(ntaps));
45        let (dst, dr) = crate::stream::new_stream();
46        (
47            Self {
48                src,
49                ntaps,
50                dst,
51                history: vec![0.0; ntaps],
52                filter: Fir::new(&taps),
53            },
54            dr,
55        )
56    }
57}
58
59impl Block for Hilbert {
60    fn work(&mut self) -> Result<BlockRet<'_>> {
61        debug_assert_eq!(self.ntaps, self.history.len());
62        let (ii, tags) = self.src.read_buf()?;
63        let i = ii.slice();
64        if i.is_empty() {
65            return Ok(BlockRet::WaitForStream(&self.src, 1));
66        }
67        let mut oo = self.dst.write_buf()?;
68        let o = oo.slice();
69        if o.is_empty() {
70            return Ok(BlockRet::WaitForStream(&self.dst, 1));
71        }
72
73        let inout = std::cmp::min(i.len(), o.len());
74        let len = self.history.len() + inout;
75        let n = len - self.ntaps;
76
77        if n == 0 {
78            // TODO: use the right number instead of 1.
79            return Ok(BlockRet::WaitForStream(
80                if i.len() < o.len() {
81                    &self.src
82                } else {
83                    &self.dst
84                },
85                1,
86            ));
87        }
88
89        // TODO: Probably needless copy.
90        let mut iv = Vec::with_capacity(len);
91        iv.extend(&self.history);
92        iv.extend(i.iter().take(inout).copied());
93
94        {
95            use rayon::prelude::*;
96            o.par_iter_mut().take(n).enumerate().for_each(|(i, val)| {
97                let t = &iv[i..(i + self.ntaps)];
98                *val = Complex::new(iv[i + self.ntaps / 2], self.filter.filter_float(t));
99            });
100        }
101
102        oo.produce(n, &tags);
103
104        self.history[..self.ntaps].clone_from_slice(&iv[n..len]);
105        ii.consume(n);
106        Ok(BlockRet::Again)
107    }
108}