nmr_schedule/schedule/
point_spread.rs

1//! Implements the Point Spread Function
2//!
3//! See [`crate::Schedule::point_spread`]
4
5use core::ops::Deref;
6
7use alloc::{borrow::ToOwned, boxed::Box, vec::Vec};
8use ndarray::Ix1;
9use once_cell::race::OnceBox;
10
11use rustfft::{num_complex::Complex, FftPlanner};
12
13use crate::DisplayMode;
14
15use super::Schedule;
16
17#[derive(Debug)]
18/// Represents a Point Spread Function; see [`crate::schedule::point_spread`]
19pub struct PointSpread {
20    unrotated: Box<[Complex<f64>]>,
21    rotated: OnceBox<Box<[Complex<f64>]>>,
22}
23
24impl Deref for PointSpread {
25    type Target = [Complex<f64>];
26
27    fn deref(&self) -> &Self::Target {
28        &self.unrotated
29    }
30}
31
32impl PointSpread {
33    /// Unwrap the PSF as a list of complex numbers. This is the direct FFT of the sampling schedule so the central peak appears at the beginning and end when plotted.
34    pub fn into_inner(self) -> Box<[Complex<f64>]> {
35        self.unrotated
36    }
37
38    /// Get the PSF in rotated form such that the "central" peak is actually in the center. This is what people normally expect a PSF to look like when visualized.
39    ///
40    /// This is not appropriate for computation because it is not the direct FFT of the sampling schedule.
41    pub fn rotated(&self) -> &[Complex<f64>] {
42        self.rotated.get_or_init(|| {
43            let mut rotated = self.unrotated.to_owned();
44
45            rotated.rotate_right(rotated.len() / 2);
46
47            rotated.into()
48        })
49    }
50
51    pub(crate) fn inner_mut(&mut self) -> &mut [Complex<f64>] {
52        self.rotated = OnceBox::new();
53        &mut self.unrotated
54    }
55
56    pub(super) fn calculate(sched: &Schedule<Ix1>) -> PointSpread {
57        let mut planner = FftPlanner::new();
58        let fft = planner.plan_fft_forward(sched.len());
59
60        let mut data = sched
61            .iter()
62            .map(|v| Complex {
63                re: if *v { 1. } else { 0. },
64                im: 0.,
65            })
66            .collect::<Vec<_>>();
67        fft.process(&mut data);
68
69        for v in data.iter_mut() {
70            *v /= (sched.len() as f64).sqrt()
71        }
72
73        PointSpread {
74            unrotated: data.into(),
75            rotated: OnceBox::new(),
76        }
77    }
78
79    /// Calculate the radius of the central peak. This method counts the number of terms that are monotonically decreasing as you walk away from the constant term, including the constant term itself.
80    ///
81    /// This is useful to ignore line-broadening of the the central peak caused by the PDF.
82    pub fn central_peak_radius(&self, mode: DisplayMode) -> usize {
83        self.windows(2)
84            .take_while(|v| mode.magnitude(v[0]) > mode.magnitude(v[1]))
85            .count()
86            + 1
87    }
88
89    /// Calculate the ratio between the magnitude of the central peak and the magnitude of the largest feature other than the central peak.
90    ///
91    /// The central peak line-broadening caused by the non-uniformity of the schedule is removed according to [`PointSpread::central_peak_radius`].
92    pub fn peak_to_sidelobe_ratio(&self, mode: DisplayMode) -> f64 {
93        let peak = mode.magnitude(self[0]);
94
95        self.iter()
96            .take(self.len() / 2 + 1)
97            .skip(self.central_peak_radius(mode) + 1)
98            .map(|v| mode.magnitude(*v))
99            .max_by(|a, b| a.total_cmp(b))
100            .unwrap_or(0.)
101            / peak
102    }
103}
104
105#[cfg(feature = "terminal-viz")]
106impl crate::terminal_viz::Plot for PointSpread {
107    type BarChart<'a> = <[Complex<f64>] as crate::terminal_viz::Plot>::BarChart<'a>
108    where
109        Self: 'a;
110
111    type RowChart<'a> = <[Complex<f64>] as crate::terminal_viz::Plot>::RowChart<'a>
112    where
113        Self: 'a;
114
115    fn bar_chart(&self, options: crate::terminal_viz::BarChartOptions) -> Self::BarChart<'_> {
116        self.rotated().bar_chart(options)
117    }
118
119    fn row_chart(&self, options: crate::terminal_viz::RowChartOptions) -> Self::RowChart<'_> {
120        self.rotated().row_chart(options)
121    }
122}
nmr_schedule/schedule/point_spread.rs

nmr_schedule/schedule/
point_spread.rs
