diff --git a/src/bin/psd.rs b/src/bin/psd.rs index 8050a1a..8d7f8a4 100644 --- a/src/bin/psd.rs +++ b/src/bin/psd.rs @@ -77,9 +77,9 @@ impl AcqOpts { fn merge_opts(&self) -> MergeOpts { MergeOpts { - remove_overlap: !self.keep_overlap, + keep_overlap: self.keep_overlap, min_count: self.avg_min, - remove_transition_band: !self.keep_transition_band, + keep_transition_band: self.keep_transition_band, } } } @@ -92,7 +92,7 @@ enum Cmd { } /// Trapezoidal integrator for irregular sampling -#[derive(Default, Copy, Clone, Debug, PartialEq, PartialOrd)] +#[derive(Default, Clone, Debug)] struct Trapezoidal { x: f32, y: f32, @@ -121,15 +121,15 @@ struct Trace { } impl Trace { - fn into_plot(self, acq: &AcqOpts) -> (String, f32, Vec<[f64; 2]>, Vec) { + fn plot(&self, acq: &AcqOpts) -> (f32, Vec<[f64; 2]>) { let logfs = acq.fs.log10(); let mut p0 = Trapezoidal::default(); let mut pi = 0.0; let plot = self .psd - .into_iter() - .zip(self.frequencies) - .filter_map(|(p, f)| { + .iter() + .zip(self.frequencies.iter()) + .filter_map(|(&p, &f)| { // TODO: check at stage breaks let dp = p0.push(f, p); if (acq.integral_start..=acq.integral_end).contains(&(acq.fs * f)) { @@ -150,7 +150,7 @@ impl Trace { } }) .collect(); - (self.name, pi.sqrt(), plot, self.breaks) + (pi.sqrt(), plot) } } @@ -277,8 +277,8 @@ impl eframe::App for App { self.current = trace .into_iter() .map(|t| { - let (name, pi, xy, b) = t.into_plot(&self.acq); - (format!("{name}: {pi:.2e}"), xy, b) + let (pi, xy) = t.plot(&self.acq); + (format!("{}: {pi:.2e}", t.name), xy, t.breaks) }) .collect(); ctx.request_repaint_after(Duration::from_secs_f32(self.repaint)); @@ -302,16 +302,13 @@ impl App { fn plot(&mut self, ui: &mut Ui) { Plot::new("stages") .link_axis("plots", true, false) - .show_axes([false, true]) - .y_axis_min_width(4.0) - .y_axis_label("X") - .y_axis_formatter(|_, _| "".to_string()) + .show_axes([false; 2]) .show_grid(false) + .show_background(false) .show_x(false) .show_y(false) .include_y(0.0) .include_y(1.0) - .show_background(false) .allow_boxed_zoom(false) .allow_double_click_reset(false) .allow_drag(false) @@ -320,7 +317,7 @@ impl App { .height(20.0) .show(ui, |plot_ui| { let ldfs = self.acq.fs.log10() as f64; - // TODO: single-trace data + // TODO: single-trace data/common modulation axis for (_name, _line, breaks) in self.current.iter().take(1) { let mut end = f32::NEG_INFINITY; let mut texts = Vec::with_capacity(breaks.len()); @@ -328,7 +325,7 @@ impl App { breaks .iter() .map(|b| { - let bins = b.bins(); + let bins = b.bins.clone(); let rbw = b.rbw(); let start = (bins.start.max(1) as f32 * rbw).log10().max(end); end = (bins.end as f32 * rbw).log10(); @@ -375,8 +372,8 @@ impl App { .x_axis_formatter(log10_axis_formatter) .link_axis("plots", false, false) .auto_bounds([true; 2].into()) - .y_axis_min_width(4.0) - .y_axis_label("Power spectral density (dB/Hz) or integrated RMS") + .y_axis_min_width(30.0) + .y_axis_label("Power spectral density (dB/Hz) or integrated RMS (1)") .legend(Legend::default()) .label_formatter(log10x_formatter) .show(ui, |plot_ui| { @@ -468,7 +465,7 @@ impl App { .on_hover_text("Do not remove bins in the filter transition band"); ui.separator(); ui.checkbox(&mut self.acq.integrate, "Integrate") - .on_hover_text("Show integrated PSD as linear cumulative sum"); + .on_hover_text("Show integrated PSD as linear cumulative sum (i.e. RMS)"); ui.separator(); ui.add( Slider::new(&mut self.acq.integral_start, 0.0..=self.acq.integral_end) diff --git a/src/psd.rs b/src/psd.rs index 0aa1506..7f89284 100644 --- a/src/psd.rs +++ b/src/psd.rs @@ -70,17 +70,12 @@ pub enum Detrend { Linear, } -impl core::fmt::Display for Detrend { - fn fmt(&self, f: &mut core::fmt::Formatter) -> core::fmt::Result { - core::fmt::Debug::fmt(self, f) - } -} - impl Detrend { pub fn apply(&self, x: &[f32; N], win: &Window) -> [Complex; N] { // apply detrending, window, make complex let mut c = [Complex::default(); N]; + // Unfortunately it doesn't optimize well with the match inside the loop match self { Detrend::None => { for ((c, x), w) in c.iter_mut().zip(x.iter()).zip(win.win.iter()) { @@ -291,7 +286,7 @@ impl PsdStage for Psd { } /// Stage break information -#[derive(Copy, Clone, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)] +#[derive(Clone, Debug, PartialEq, Eq, Hash)] pub struct Break { /// Start index in PSD and frequencies pub start: usize, @@ -302,10 +297,10 @@ pub struct Break { /// Averaging limit pub avg: u32, /// FFT bin - pub bins: (usize, usize), + pub bins: Range, /// FFT size pub fft_size: usize, - /// The decimation power of two + /// The decimation pub decimation: usize, /// Unprocessed number of input samples (includes overlap) pub pending: usize, @@ -321,7 +316,7 @@ impl Break { .filter(|bi| bi.include) .flat_map(|bi| { let rbw = bi.rbw(); - bi.bins().map(move |f| f as f32 * rbw) + bi.bins.clone().map(move |f| f as f32 * rbw) }) .collect(); debug_assert_eq!(f.first(), Some(&0.0)); @@ -330,36 +325,32 @@ impl Break { } pub fn effective_fft_size(&self) -> usize { - self.fft_size << self.decimation + self.fft_size * self.decimation } /// Absolute resolution bandwidth pub fn rbw(&self) -> f32 { 1.0 / self.effective_fft_size() as f32 } - - pub fn bins(&self) -> Range { - self.bins.0..self.bins.1 - } } /// PSD segment merge options #[derive(Copy, Clone, Debug, PartialEq, Eq)] pub struct MergeOpts { - /// Remove low resolution bins - pub remove_overlap: bool, + /// Keep low resolution bins + pub keep_overlap: bool, /// Minimum averaging level pub min_count: u32, - /// Remove decimation filter transition bands - pub remove_transition_band: bool, + /// Keep decimation filter transition bands + pub keep_transition_band: bool, } impl Default for MergeOpts { fn default() -> Self { Self { - remove_overlap: true, + keep_overlap: false, min_count: 1, - remove_transition_band: true, + keep_transition_band: false, } } } @@ -488,26 +479,27 @@ impl PsdCascade { pub fn psd(&self, opts: &MergeOpts) -> (Vec, Vec) { let mut p = Vec::with_capacity(self.stages.len() * (N / 2 + 1)); let mut b = Vec::with_capacity(self.stages.len()); - let mut decimation = self - .stages - .iter() - .rev() - .map(|stage| stage.hbf.depth()) - .sum(); + let mut decimation = 1 + << self + .stages + .iter() + .rev() + .map(|stage| stage.hbf.depth()) + .sum::(); let mut end = 0; for stage in self.stages.iter().rev() { - decimation -= stage.hbf.depth(); + decimation >>= stage.hbf.depth(); // a stage yields frequency bins 0..N/2 from DC up to its nyquist // 0..floor(0.4*N) is its passband if it was preceeded by a decimator // 0..floor(0.4*N)/R is the passband of the next lower stage // hence take bins ceil(floor(0.4*N)/R)..floor(0.4*N) from a non-edge stage - let start = if opts.remove_overlap { + let start = if !opts.keep_overlap { // remove low f bins, ceil (end + (1 << stage.hbf.depth()) - 1) >> stage.hbf.depth() } else { 0 }; - end = if decimation > 0 && opts.remove_transition_band { + end = if decimation > 1 && !opts.keep_transition_band { // remove transition band of higher stage's decimator 2 * N / 5 // 0.4, floor } else { @@ -519,7 +511,7 @@ impl PsdCascade { count: stage.count(), include, avg: stage.avg, - bins: (start, end), + bins: start..end, fft_size: N, decimation, processed: N * stage.count() as usize @@ -527,7 +519,7 @@ impl PsdCascade { pending: stage.buf().len(), }); if include { - let g = (1 << decimation) as f32 / stage.gain(); + let g = decimation as f32 / stage.gain(); p.extend(stage.spectrum()[start..end].iter().map(|pi| pi * g)); } else { end = start; @@ -590,9 +582,9 @@ mod test { let mut s = PsdCascade::::new(1); s.process(&x); let merge_opts = MergeOpts { - remove_overlap: false, + keep_overlap: true, min_count: 0, - remove_transition_band: true, + ..Default::default() }; let (p, b) = s.psd(&merge_opts); let f = Break::frequencies(&b); diff --git a/src/source.rs b/src/source.rs index 059b8a0..2dd8410 100644 --- a/src/source.rs +++ b/src/source.rs @@ -107,6 +107,7 @@ impl Source { .map(|mut x| { x = (x - 0.5) * 12.0f32.sqrt(); // zero mean, RMS = 1 state.iter_mut().fold(x, |mut x, s| { + // TODO: branch optimization barrier (x, *s) = if *diff { (x - *s, x) } else { (*s, x + *s) }; x })