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    /// TODO: Citation
90    ///
91    /// Calculate the ratio between the magnitude of the central peak and the magnitude of the largest feature other than the central peak.
92    ///
93    /// The central peak line-broadening caused by the non-uniformity of the schedule is removed according to [`PointSpread::central_peak_radius`].
94    pub fn peak_to_sidelobe_ratio(&self, mode: DisplayMode) -> f64 {
95        let peak = mode.magnitude(self[0]);
96
97        self.iter()
98            .take(self.len() / 2 + 1)
99            .skip(self.central_peak_radius(mode) + 1)
100            .map(|v| mode.magnitude(*v))
101            .max_by(|a, b| a.total_cmp(b))
102            .unwrap_or(0.)
103            / peak
104    }
105}
106
107#[cfg(feature = "terminal-viz")]
108impl crate::terminal_viz::Plot for PointSpread {
109    type BarChart<'a> = <[Complex<f64>] as crate::terminal_viz::Plot>::BarChart<'a>
110    where
111        Self: 'a;
112
113    type RowChart<'a> = <[Complex<f64>] as crate::terminal_viz::Plot>::RowChart<'a>
114    where
115        Self: 'a;
116
117    fn bar_chart(&self, options: crate::terminal_viz::BarChartOptions) -> Self::BarChart<'_> {
118        self.rotated().bar_chart(options)
119    }
120
121    fn row_chart(&self, options: crate::terminal_viz::RowChartOptions) -> Self::RowChart<'_> {
122        self.rotated().row_chart(options)
123    }
124}
nmr_schedule/schedule/point_spread.rs

nmr_schedule/schedule/
point_spread.rs
