Add stft

Author: sigurd4

src/analysis/cceps.rs

@@ -10,16 +10,26 @@ use crate::{quantities::{List, ListOrSingle, Lists}, util::TruncateIm};
 #[derive(Debug, Clone, Copy, PartialEq, Error)]
 pub enum CepsError
 {
-    #[error("Sequence has one or more zero-valued fourier coefficient.")]
+    #[error("Sequence has one or more zero-valued fourier coefficients.")]
     ZeroInFourier
 }
 
+/// A trait for calculating the complex cepstrum of a sequence or several sequences.
 pub trait CCeps<'a, T, C, N>: Lists<T>
 where
     T: ComplexFloat,
     C: List<T>,
     N: Maybe<usize>,
 {
+    /// Computes the complex cepstrum of a sequence or several sequences.
+    /// 
+    /// # Arguments
+    /// 
+    /// * `numtaps`: Number of taps if cepstrum is dynamically sized, otherwise inferred from return type.
+    /// 
+    /// # Returns
+    /// 
+    /// * `cepstrum`: Complex cepstrum in the form of a sequence.
     fn cceps(&'a self, numtaps: N) -> Result<Self::RowsMapped<C>, CepsError>;
 }
 

src/analysis/cpsd.rs

@@ -3,6 +3,7 @@ use option_trait::{Maybe, NotVoid, StaticMaybe};
 
 use crate::{quantities::List, util::MaybeLenEq, analysis::{PWelch, PWelchDetrend}};
 
+/// A trait for computing the cross power spectral density of two sequences.
 pub trait CPsd<T, Y, YY, W, WW, WWW, WL, N, S>: List<T> + MaybeLenEq<YY, true>
 where
     T: ComplexFloat,
@@ -15,6 +16,25 @@ where
     N: Maybe<usize>,
     S: Maybe<bool>
 {
+    /// Computes the cross power spectral density of two sequences using Welch's averaged periodogram method.
+    /// 
+    /// # Arguments
+    /// 
+    /// * `y` - The other sequence.
+    /// * `window` - A window sequence. If none, a [Hamming](crate::windows::Hamming) window with lenght `window_length` will be used.
+    /// * `window_length` - A window length if no `window` is given.
+    /// * `overlap` - Overlap length in samples. If none, defaults to `window_length/2`.
+    /// * `nfft` - Length used for FFT if no `window` is given.
+    /// * `sampling_frequency` - Sampling frequency. If none, defaults to `1.0`.
+    /// * `confidence` - Probability in the range `[0.0, 1.0]` for the confidence intervals for the PSD estimate. If none, defaults to `0.95`.
+    /// * `detrend` - Method of detrending the signals before or after analysis. If none, [LongMean](PWelchDetrend::LongMean) is used.
+    /// * `sloppy` - If true, `nfft` will be set to the next power of two. Only applicable if no `window` is given.
+    /// * `shift` - If true, data will be shifted to center-DC.
+    /// 
+    /// # Returns
+    /// 
+    /// * `pxy` - Cross power spectral density.
+    /// * `frequencies` - Frequencies of the cross power spectral density.
     #[doc(alias = "csd")]
     fn cpsd<O, FS, CONF, DT, F>(
         self,
@@ -114,9 +134,9 @@ mod test
 
         const M: usize = 500;
         let w: [_; M] = Triangular.window_gen((), WindowRange::Symmetric);
-        let nov = 250;
+        let noverlap = 250;
 
-        let (pxy, f): (_, [_; M/2 + 1]) = x.real_cpsd(y, w, (), nov, (), (), (), (), ());
+        let (pxy, f): (_, [_; M/2 + 1]) = x.real_cpsd(y, w, (), noverlap, (), (), (), (), ());
 
         plot::plot_curves("P_xy(e^jw)", "plots/pxy_z_cpsd.png", [
             &f.zip(pxy.map(|p| 10.0*p.norm().log10()))

src/analysis/filternorm.rs

@@ -9,11 +9,21 @@ use crate::{
     systems::Tf
 };
 
+/// A trait for computing the `Lp`-norm or infinity-norm of a digital filter.
 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;
+    /// Computes the `Lp`-norm or infinity-norm of a digital filter.
+    /// 
+    /// # Arguments
+    /// 
+    /// * `pnorm` - The norm power. If infinite, the infinity-norm of the filter is returned.
+    /// 
+    /// # Returns
+    /// 
+    /// * `norm` - The `Lp`-norm or infinity-norm of the digital filter.
+    fn filternorm(&'a self, pnorm: <Self::Set as ComplexFloat>::Real) -> Self::Output;
 }
 
 const FILTER_INF_NORM_RES: usize = 1024;
@@ -29,15 +39,15 @@ where
 {
     type Output = B::RowsMapped<T::Real>;
 
-    fn filternorm(&'a self, p: T::Real) -> Self::Output
+    fn filternorm(&'a self, pnorm: T::Real) -> Self::Output
     {
-        if Float::abs(p) > Float::sqrt(<T::Real as Float>::max_value())
+        if Float::abs(pnorm) > 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})
+                .reduce(if pnorm.is_sign_positive() {Float::max} else {Float::min})
                 .unwrap_or_else(Zero::zero)
             )
         }
@@ -48,9 +58,9 @@ where
             let n = h.map_rows_to_owned(|h| {
                 let n = Float::powf(h.as_view_slice()
                     .iter()
-                    .map(|&h| Float::powf(h.abs(), p))
+                    .map(|&h| Float::powf(h.abs(), pnorm))
                     .sum::<T::Real>(),
-                    Float::recip(p)
+                    Float::recip(pnorm)
                 );
                 if Float::is_infinite(n)
                 {
@@ -60,7 +70,7 @@ where
             });
             if inf
             {
-                return self.filternorm(Float::copysign(Float::infinity(), p))
+                return self.filternorm(Float::copysign(Float::infinity(), pnorm))
             }
             n
         }

src/analysis/specgram.rs

@@ -49,9 +49,9 @@ where
         let x: &[T] = self.as_view_slice();
         let window: &[W] = window.as_view_slice();
         let n = window.len();
-        let step = n - overlap;
+        let step = n - overlap % n;
 
-        let s: Vec<_> = (0..x.len() - n).step_by(step)
+        let s: Vec<_> = (0..x.len() + 1 - n).step_by(step)
             .map(|i| {
                 let mut y: Vec<_> = x[i..].iter()
                     .zip(window.iter())
@@ -65,7 +65,7 @@ where
         let fs = sampling_frequency.into_option()
             .unwrap_or_else(FloatConst::TAU);
 
-        let t = (0..x.len() - n).step_by(step)
+        let t = (0..x.len() + 1 - n).step_by(step)
             .map(|i| <T::Real as NumCast>::from(i).unwrap()/fs)
             .collect();
         let f = (0..n).map(|i| <T::Real as NumCast>::from(i).unwrap()/NumCast::from(n).unwrap()*fs)
@@ -82,7 +82,7 @@ where
     WW: List<W, Length = usize>,
     X: List<T, Length = usize>,
     Complex<T::Real>: AddAssign + MulAssign,
-    //[(); 0 - (X::LENGTH - WW::LENGTH) % L]:
+    [(); WW::LENGTH - (X::LENGTH - WW::LENGTH)/L - 1]:
 {
     fn specgram<FS>(
         &self,
@@ -157,7 +157,8 @@ where
     X: List<T, Length = usize>,
     [T::Real; W]: List<T::Real, Length = usize>,
     Complex<T::Real>: AddAssign + MulAssign,
-    [(); <[T::Real; W]>::LENGTH]:
+    [(); <[T::Real; W]>::LENGTH]:,
+    [(); <[T::Real; W]>::LENGTH - (X::LENGTH - <[T::Real; W]>::LENGTH)/L - 1]:
 {
     fn specgram<FS>(
         &self,

src/transforms/fourier/mod.rs

@@ -19,6 +19,7 @@ moddef::moddef!(
         idwht,
         sdft_kw,
         sdft,
+        stft,
         swift
     }
 );