signal_processing 0.3.0

A signal processing library.
Documentation
use core::{iter::Sum, ops::Mul};

use ndarray::Array2;
use num::{complex::ComplexFloat, traits::FloatConst, Complex, NumCast, Zero};
use option_trait::Maybe;

use crate::quantities::{List, Matrix};

pub trait ToKw<T, K, N>: List<T>
where
    T: ComplexFloat,
    K: Matrix<Complex<T::Real>>,
    N: Maybe<usize>
{
    fn to_kw(self, n: N) -> K;
}

impl<T, W, const N: usize> ToKw<T, [[Complex<T::Real>; N]; N], ()> for W
where
    T: ComplexFloat,
    W: List<T>,
    Complex<T::Real>: Mul<T, Output = Complex<T::Real>> + Sum
{
    fn to_kw(self, (): ()) -> [[Complex<T::Real>; N]; N]
    {
        let mf = <T::Real as NumCast>::from(self.length()).unwrap();
        let nf = <T::Real as NumCast>::from(N).unwrap();
        let mut kw = [[Complex::zero(); N]; N];
        for i in 1..2*N
        {
            let dk = <T::Real as NumCast>::from(i).unwrap() - nf;
            let w = self.as_view_slice()
                .iter()
                .enumerate()
                .map(|(i, w)| {
                    let i = <T::Real as NumCast>::from(i).unwrap();
                    Complex::cis(T::Real::TAU()*i/mf*dk)**w
                }).sum::<Complex<T::Real>>()/mf;
            for k_ in i.saturating_sub(N)..i.min(N)
            {
                let k = k_ + N - i;
                kw[k][k_] = w;
            }
        }
        kw
    }
}

impl<T, W> ToKw<T, Array2<Complex<T::Real>>, usize> for W
where
    T: ComplexFloat,
    W: List<T>,
    Complex<T::Real>: Mul<T, Output = Complex<T::Real>> + Sum
{
    fn to_kw(self, n: usize) -> Array2<Complex<T::Real>>
    {
        let mf = <T::Real as NumCast>::from(self.length()).unwrap();
        let nf = <T::Real as NumCast>::from(n).unwrap();
        let mut kw = Array2::from_elem((n, n), Complex::zero());
        for i in 1..2*n
        {
            let dk = <T::Real as NumCast>::from(i).unwrap() - nf;
            let w = self.as_view_slice()
                .iter()
                .enumerate()
                .map(|(i, w)| {
                    let i = <T::Real as NumCast>::from(i).unwrap();
                    Complex::cis(T::Real::TAU()*i/mf*dk)**w
                }).sum::<Complex<T::Real>>()/mf;
            for k_ in i.saturating_sub(n)..i.min(n)
            {
                let k = k_ + n - i;
                kw[(k, k_)] = w;
            }
        }
        kw
    }
}

#[cfg(test)]
mod test
{
    use array_math::Array2dOps;

    use crate::{plot, windows::Hann, gen::window::{WindowGen, WindowRange}, transforms::window::ToKw};

    #[test]
    fn test()
    {
        const W: usize = 1024;
        const N: usize = 10;
        let w: [f64; W] = Hann.window_gen((), WindowRange::Symmetric);
        let kw: [[_; N]; _] = w.to_kw(());

        println!("{:?}", kw.diagonal_ref());

        plot::plot_parametric_curve_2d("K[i, j]", "plots/k_ij_to_kw.svg",
            core::array::from_fn::<_, N, _>(|i| i as f64),
            core::array::from_fn::<_, N, _>(|i| i as f64),
            |i, j| [i, j, kw[i as usize][j as usize].norm().log10()*20.0]

        ).unwrap()
    }
}