welch_sde/
lib.rs

1//! # Spectral density and power spectrum estimation with Welch method
2//!
3//! The Welch's method for the estimation of the periodogram of a *stationnary and
4//! zero-mean* signal  involves splitting the signal in overlapping segments, multiplying
5//! each segment with a predeterming window and computing the discrete Fourier transform
6//! of each windowed segment.
7//! The periodogram is then given by the average of the squared magnitude
8//! of each Fourier transform.
9//! Only the halve of the power spectrum that corresponds to positive frequencies is returned,
10//! both halves beeing symmetric with respect to the zero frequency.
11//!
12//!
13//! From the periodogram, one can derive either the **spectral density** or the **power spectrum**.
14//! Both differs with respect to the scaling of the periodogram.
15//! For the **spectral density**, the periodogram is divided by the product of the sampling frequency with the sum of the squared window samples.
16//! For the **power spectrum**, the periodogram is divided by the square of the sum of the window samples.
17//!
18//! The Welch algorithm is implemented in the [Welch] structure.
19//! For convenience, 2 new types (that encapsulates [Welch]) are provided, [SpectralDensity] and [PowerSpectrum], to compute  the **spectral density** and the **power spectrum**, respectively.
20//! [SpectralDensity] uses a [Hann] window whereas [PowerSpectrum] does not use any window.
21//!
22//! Custom windows can be used with [Welch] if they implement the [Window] trait.
23//! The signal is either a [single](f32) or [double](f64) floating point array.
24//!
25//! ## Examples
26//! ### Power spectrum
27//!```
28//!use rand::prelude::*;
29//!use rand_distr::StandardNormal;
30//!use std::time::Instant;
31//!use welch_sde::{Build, PowerSpectrum};
32//!
33//!fn main() {
34//!    let n = 1e6 as usize;
35//!    let signal: Vec<f64> = (0..n)
36//!        .map(|_| 1.234_f64.sqrt() * thread_rng().sample::<f64, StandardNormal>(StandardNormal))
37//!        .collect();
38//!    {
39//!        let mean = signal.iter().cloned().sum::<f64>() / n as f64;
40//!        let variance = signal.iter().map(|s| *s - mean).map(|x| x * x).sum::<f64>() / n as f64;
41//!        println!("Signal variance: {:.3}", variance);
42//!    }
43//!
44//!    let welch: PowerSpectrum<f64> = PowerSpectrum::builder(&signal).build();
45//!    println!("{}", welch);
46//!
47//!    let now = Instant::now();
48//!    let ps = welch.periodogram();
49//!    println!(
50//!        "Power spectrum estimated in {}ms",
51//!        now.elapsed().as_millis()
52//!    );
53//!    {
54//!        let variance = 2. * ps.iter().sum::<f64>();
55//!        println!("Signal variance from power spectrum: {:.3}", variance);
56//!    }
57//!}
58//!```
59//! ### Spectral density
60//! ```
61//!use rand::prelude::*;
62//!use rand_distr::StandardNormal;
63//!use std::time::Instant;
64//!use welch_sde::{SpectralDensity, Build};
65//!
66//!fn main() {
67//!    let n = 1e5 as usize;
68//!    let fs = 10e3_f64;
69//!    let amp = 2. * 2f64.sqrt();
70//!    let freq = 1550f64;
71//!    let noise_power = 0.001 * fs / 2.;
72//!    let sigma = noise_power.sqrt();
73//!    let signal: Vec<f64> = (0..n)
74//!        .map(|i| i as f64 / fs)
75//!        .map(|t| {
76//!            amp * (2. * std::f64::consts::PI * freq * t).sin()
77//!                + thread_rng().sample::<f64, StandardNormal>(StandardNormal) * sigma
78//!        })
79//!        .collect();
80//!
81//!    let welch: SpectralDensity<f64> =
82//!        SpectralDensity::<f64>::builder(&signal, fs).build();
83//!    println!("{}", welch);
84//!    let now = Instant::now();
85//!    let sd = welch.periodogram();
86//!    println!(
87//!        "Spectral density estimated in {}ms",
88//!        now.elapsed().as_millis()
89//!    );
90//!    let noise_floor = sd.iter().cloned().sum::<f64>() / sd.len() as f64;
91//!    println!("Noise floor: {:.3}", noise_floor);
92//!
93//!    let _: complot::LinLog = (
94//!        sd.frequency()
95//!            .into_iter()
96//!            .zip(&(*sd))
97//!            .map(|(x, &y)| (x, vec![y])),
98//!        complot::complot!(
99//!            "spectral_density.png",
100//!            xlabel = "Frequency [Hz]",
101//!            ylabel = "Spectral density [s^2/Hz]"
102//!        ),
103//!    )
104//!        .into();
105//!}
106//!```
107
108mod builder;
109mod periodogram;
110mod power_spectrum;
111mod spectral_density;
112mod welch;
113mod window;
114pub use builder::Builder;
115use num_traits::Float;
116pub use periodogram::{Periodogram, PowerSpectrumPeriodogram, SpectralDensityPeriodogram};
117pub use power_spectrum::PowerSpectrum;
118use rustfft::FftNum;
119pub use spectral_density::SpectralDensity;
120pub use welch::Welch;
121pub use window::{Hann, One, Window};
122
123/// The trait the signal type `T` must implement
124pub trait Signal:
125    Float + FftNum + std::iter::Sum + std::ops::SubAssign + std::ops::AddAssign
126{
127}
128impl Signal for f64 {}
129impl Signal for f32 {}
130
131/// [Builder] trait
132pub trait Build<T: Signal, W: Window<T>, E> {
133    /// Returns a struct `E` initialized according to the [Builder] settings
134    fn build(&self) -> E;
135}
welch_sde/lib.rs

welch_sde/
lib.rs
