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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
#![feature(new_uninit)]
#![feature(is_sorted)]
#![feature(generic_arg_infer)]
#![feature(let_chains)]
#![feature(array_methods)]
#![feature(trait_alias)]
#![feature(unsize)]
#![feature(associated_type_bounds)]
#![feature(unboxed_closures)]
#![feature(array_chunks)]
#![feature(associated_type_defaults)]
#![feature(auto_traits)]
#![feature(negative_impls)]
#![feature(extract_if)]
#![feature(associated_const_equality)]
#![feature(split_array)]
#![feature(iterator_try_collect)]
#![feature(slice_flatten)]
#![feature(float_gamma)]
#![feature(fn_traits)]
#![feature(lazy_cell)]
#![feature(coerce_unsized)]
#![feature(decl_macro)]
#![feature(array_try_map)]

#![allow(internal_features)]
#![allow(incomplete_features)]

#![feature(adt_const_params)]
#![feature(core_intrinsics)]
#![feature(inherent_associated_types)]
#![feature(generic_const_exprs)]
#![feature(specialization)]

moddef::moddef!(
    flat(pub) mod {
        analysis,
        decompositions,
        error,
        gen,
        identification,
        operations,
        quantities,
        systems,
        transforms,
        util,

        rtf_or_system,
        system,
        validate_filter_bands
    },
    pub mod {
        window,
    },
    mod {
        plot for cfg(test)
    }
);

#[cfg(test)]
mod tests
{
    use array_math::ArrayOps;
    use linspace::LinspaceArray;
    use num::Complex;

    use crate::{plot, BesselAP, BesselF, ButtAP, Butter, Cheb1AP, Cheb2AP, EllipAP, FilterGenPlane, FilterGenType, FreqS, FreqZ, RealFreqZ, Tf, Zpk};

    #[test]
    fn testt()
    {
        let fs = 1e3f64;
        let f_p = 50.0;
        let f_s = 10.0;
        let dp = 3.0;
        let ds = 80.0;

        let (n, wn, _, t) = crate::buttord([f_p], [f_s], dp, ds, FilterGenPlane::Z { sampling_frequency: Some(fs) })
            .unwrap();
        let ba = Tf::butter(n, wn, t, FilterGenPlane::Z { sampling_frequency: None })
            .unwrap();

        const N: usize = 1024;
        let (h, w): ([_; N], _) = ba.real_freqz(());

        plot::plot_curves("H(e^jw)", "temp/butter_hf.png", [&w.zip(h.map(|h| 10.0*h.norm_sqr().log10()))])
            .unwrap()
    }

    #[test]
    fn test_ap()
    {
        let fs = 2.0f64;
        let o = 6;
        let h_bessel = Zpk::besselap(o);
        let h_butter = Zpk::buttap(o);
        let h_cheb1 = Zpk::cheb1ap(o, 5.0);
        let h_cheb2 = Zpk::cheb2ap(o, 5.0);
        let h_ellip = Zpk::ellipap(o, 5.0, 5.0);
        
        const N: usize = 1024;
        let w: [_; N] = (0.0..fs).linspace_array();
        let h_f_bessel = h_bessel.freqs(w.map(|w| Complex::new(0.0, w)));
        let h_f_butter = h_butter.freqs(w.map(|w| Complex::new(0.0, w)));
        let h_f_cheb1 = h_cheb1.freqs(w.map(|w| Complex::new(0.0, w)));
        let h_f_cheb2 = h_cheb2.freqs(w.map(|w| Complex::new(0.0, w)));
        let h_f_ellip = h_ellip.freqs(w.map(|w| Complex::new(0.0, w)));

        plot::plot_curves("H(jw)", "plots/h_s_ap.png",
            [
                &w.zip(h_f_bessel.map(|h| 10.0*h.norm_sqr().log10())),
                &w.zip(h_f_butter.map(|h| 10.0*h.norm_sqr().log10())),
                &w.zip(h_f_cheb1.map(|h| 10.0*h.norm_sqr().log10())),
                &w.zip(h_f_cheb2.map(|h| 10.0*h.norm_sqr().log10())),
                &w.zip(h_f_ellip.map(|h| 10.0*h.norm_sqr().log10())),
            ]
        ).unwrap();
    }

    #[test]
    fn test()
    {
        let fs: f64 = 44100.0;
        let h1 = Tf::<_, [_; 2], [_; 2]>::besself((), [800.0], FilterGenType::HighPass, crate::FilterGenPlane::Z { sampling_frequency: Some(fs) }).unwrap();
        let h2 = Tf::<_, [_; 5], [_; 5]>::butter((), [8000.0], FilterGenType::LowPass, crate::FilterGenPlane::Z { sampling_frequency: Some(fs) }).unwrap();
        let h3 = Tf::<_, [_; 3], [_; 3]>::butter((), [20000.0], FilterGenType::HighPass, crate::FilterGenPlane::Z { sampling_frequency: Some(fs) }).unwrap();

        let h = h1*&h2 + h3;

        const N: usize = 1024;

        //let h: Sos<f64, _> = h.to_sos((), ());
        let (h_f, w): ([_; N], [_; N]) = h.freqz((), false);

        plot::plot_curves("H(e^jw)", "plots/h_z.png", [&w.zip(h_f.map(|h| h.norm())), &w.zip(h_f.map(|h| h.arg()))]).unwrap();

        const M: usize = 256;

        let z_r = Complex::new(1.0, 1.0);
        let z = (0..M).map(|i| (0..M).map(|j| Complex::new((2.0*i as f64/(M - 1) as f64 - 1.0)*z_r.re, (2.0*j as f64/(M - 1) as f64 - 1.0)*z_r.im)).collect()).collect();
        let h_z: Vec<Vec<Complex<f64>>> = FreqS::<Complex<f64>, Vec<Vec<Complex<f64>>>>::freqs(&h, z);

        let o = <[_; M]>::fill(|j| (2.0*j as f64/(M - 1) as f64 - 1.0)*z_r.re);
        let w = <[_; M]>::fill(|i| (2.0*i as f64/(M - 1) as f64 - 1.0)*z_r.im);

        let i: [f64; M] = (0..M).linspace_array().map(|i| i as f64);

        plot::plot_parametric_curve_2d("H(z)", "plots/h_z.svg",
            i,
            i,
            |i, j| [w[i as usize], o[j as usize], (h_z[i as usize][j as usize].norm().log10()*20.0).min(10.0).max(-10.0)]
        ).unwrap();
    }
}