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

use core::iter::Sum;

use num::{complex::ComplexFloat, Complex, Float, Zero};

use crate::{
    analysis::{FreqZ, ImpZ},
    quantities::{List, ListOrSingle, Lists, MaybeList, MaybeLists},
    System,
    systems::Tf
};

pub trait FilterNorm<'a>: System
{
    type Output: ListOrSingle<<Self::Set as ComplexFloat>::Real>;

    fn filternorm(&'a self, p: <Self::Set as ComplexFloat>::Real) -> Self::Output;
}

const FILTER_INF_NORM_RES: usize = 1024;

impl<'a, T, B, A> FilterNorm<'a> for Tf<T, B, A>
where
    T: ComplexFloat<Real: Sum>,
    B: MaybeLists<T>,
    A: MaybeList<T>,
    B::RowsMapped<Vec<T>>: for<'b> Lists<T, RowView<'b>: List<T>, RowsMapped<T::Real> = B::RowsMapped<T::Real>>,
    B::RowsMapped<[Complex<T::Real>; FILTER_INF_NORM_RES]>: for<'b> Lists<Complex<T::Real>, RowView<'b>: List<Complex<T::Real>>, RowsMapped<T::Real> = B::RowsMapped<T::Real>>, 
    Self: ImpZ<'a, B::RowsMapped<Vec<T>>, Vec<T::Real>, ()> + FreqZ<'a, B::RowsMapped<[Complex<T::Real>; FILTER_INF_NORM_RES]>, [T::Real; FILTER_INF_NORM_RES], ()> + System<Set = T>
{
    type Output = B::RowsMapped<T::Real>;

    fn filternorm(&'a self, p: T::Real) -> Self::Output
    {
        if Float::abs(p) > Float::sqrt(<T::Real as Float>::max_value())
        {
            let (h, _): (_, [_; FILTER_INF_NORM_RES]) = self.freqz((), false);
            h.map_rows_to_owned(|h| h.as_view_slice()
                .iter()
                .map(|&h| h.abs())
                .reduce(if p.is_sign_positive() {Float::max} else {Float::min})
                .unwrap_or_else(Zero::zero)
            )
        }
        else
        {
            let (h, _) = self.impz((), None);
            let mut inf = false;
            let n = h.map_rows_to_owned(|h| {
                let n = Float::powf(h.as_view_slice()
                    .iter()
                    .map(|&h| Float::powf(h.abs(), p))
                    .sum::<T::Real>(),
                    Float::recip(p)
                );
                if Float::is_infinite(n)
                {
                    inf = true;
                }
                n
            });
            if inf
            {
                return self.filternorm(Float::copysign(Float::infinity(), p))
            }
            n
        }
    }
}

#[cfg(test)]
mod test
{
    use crate::{
        gen::filter::{Butter, FilterGenPlane, FilterGenType},
        analysis::FilterNorm,
        systems::Tf
    };

    #[test]
    fn test()
    {
        let fs = 1000.0;
        
        let h: Tf::<f64, _, _> = Tf::butter(15, [220.0], FilterGenType::LowPass, FilterGenPlane::Z { sampling_frequency: Some(fs) })
            .unwrap();
        
        let n2 = h.filternorm(2.0);

        println!("n2 = {}", n2);
        
        let n_inf = h.filternorm(f64::INFINITY);

        println!("n_inf = {}", n_inf);
    }
}