signal_processing 0.3.0

A signal processing library.
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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

use core::ops::AddAssign;

use num::{Float, NumCast};
use option_trait::Maybe;

use crate::{quantities::{IntoList, List, ListOrSingle, OwnedListOrSingle}, transforms::bspline::BSplineEval};

pub trait QSpline1dEval<T, X, N, C>: IntoList<T, X, N>
where
    T: Float,
    X: ListOrSingle<T>,
    N: Maybe<usize>,
    C: List<T>
{
    fn qspline_1d_eval(self, numtaps: N, c: C) -> (X::Mapped<T>, X);
}

impl<T, C, X, XX, N> QSpline1dEval<T, X, N, C> for XX
where
    T: Float + AddAssign,
    XX: IntoList<T, X, N>,
    X: ListOrSingle<T>,
    X::Mapped<T>: OwnedListOrSingle<T>,
    N: Maybe<usize>,
    C: List<T>
{
    fn qspline_1d_eval(self, numtaps: N, c: C) -> (X::Mapped<T>, X)
    {
        let x = self.into_list(numtaps);
        let mut res = x.map_to_owned(|_| T::zero());
        if res.length() == 0
        {
            return (res, x)
        }
        let rr = res.as_mut_slice();

        let xx = x.as_view_slice();

        let c = c.into_vec();
        let n = c.len();
        let nf = T::from(n).unwrap();

        let zero = T::zero();
        let one = T::one();
        let two = one + one;
        let nm1 = nf - one;

        let cond1: Vec<_> = x.as_view_slice()
            .iter()
            .enumerate()
            .filter(|(_, x)| *x < &zero)
            .map(|(i, _)| i)
            .collect();
        let cond2: Vec<_> = x.as_view_slice()
            .iter()
            .enumerate()
            .filter(|(_, x)| *x > &nm1)
            .map(|(i, _)| i)
            .collect();
        let cond3: Vec<_> = x.as_view_slice()
            .iter()
            .enumerate()
            .filter(|(_, x)| !(*x < &zero || *x > &nm1))
            .map(|(i, _)| i)
            .collect();

        if cond1.len() > 0
        {
            let x: Vec<_> = cond1.iter()
                .map(|&i| -xx[i])
                .collect();
            let (y, _) = x.qspline_1d_eval((), c.as_slice());
            for (i, y) in cond1.into_iter()
                .zip(y)
            {
                rr[i] = y;
            }
        }
        if cond2.len() > 0
        {
            let x: Vec<_> = cond2.iter()
                .map(|&i| two*nm1 - xx[i])
                .collect();
            let (y, _) = x.qspline_1d_eval((), c.as_slice());
            for (i, y) in cond2.into_iter()
                .zip(y)
            {
                rr[i] = y;
            }
        }
        if cond3.len() > 0
        {
            let quadratic = |x: T| {
                if x.abs() > T::from(1.5).unwrap()
                {
                    return zero
                }
                const TT: [f64; 4] = [-1.5, -0.5, 0.5, 1.5];
                const CC: [f64; TT.len()*3 - 4] = const {
                    let mut c = [0.0; TT.len()*3 - 4];
                    c[TT.len() - 2] = 1.0;
                    c
                };
                x.bspline_eval((), TT.map(|t| T::from(t).unwrap()), CC.map(|c| T::from(c).unwrap()), TT.len() - 2).0
            };

            let x: Vec<_> = cond3.iter()
                .map(|&i| xx[i])
                .collect();
            let jlower: Vec<_> = x.iter()
                .map(|&x| <isize as NumCast>::from((x - T::from(1.5).unwrap()).floor() + one).unwrap())
                .collect();
            for i in 0..3
            {
                for ((r, &x), &j) in cond3.iter()
                    .copied()
                    .zip(x.iter())
                    .zip(jlower.iter())
                {
                    let thisj: isize = j + i;
                    let indj = (thisj.max(0) as usize).min(n - 1);
                    rr[r] += c[indj]*quadratic(x - T::from(thisj).unwrap())
                }
            }
        }

        (res, x)
    }
}