nmr-schedule 0.2.1

Algorithms for NMR Non-Uniform Sampling
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

//! Implements the Point Spread Function
//!
//! See [`crate::Schedule::point_spread`]

use core::ops::Deref;

use alloc::{borrow::ToOwned, boxed::Box, vec::Vec};
use ndarray::Ix1;
use once_cell::race::OnceBox;

use rustfft::{num_complex::Complex, FftPlanner};

use crate::DisplayMode;

use super::Schedule;

#[derive(Debug)]
/// Represents a Point Spread Function; see [`crate::schedule::point_spread`]
pub struct PointSpread {
    unrotated: Box<[Complex<f64>]>,
    rotated: OnceBox<Box<[Complex<f64>]>>,
}

impl Deref for PointSpread {
    type Target = [Complex<f64>];

    fn deref(&self) -> &Self::Target {
        &self.unrotated
    }
}

impl PointSpread {
    /// 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.
    pub fn into_inner(self) -> Box<[Complex<f64>]> {
        self.unrotated
    }

    /// 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.
    ///
    /// This is not appropriate for computation because it is not the direct FFT of the sampling schedule.
    pub fn rotated(&self) -> &[Complex<f64>] {
        self.rotated.get_or_init(|| {
            let mut rotated = self.unrotated.to_owned();

            rotated.rotate_right(rotated.len() / 2);

            rotated.into()
        })
    }

    pub(crate) fn inner_mut(&mut self) -> &mut [Complex<f64>] {
        self.rotated = OnceBox::new();
        &mut self.unrotated
    }

    pub(super) fn calculate(sched: &Schedule<Ix1>) -> PointSpread {
        let mut planner = FftPlanner::new();
        let fft = planner.plan_fft_forward(sched.len());

        let mut data = sched
            .iter()
            .map(|v| Complex {
                re: if *v { 1. } else { 0. },
                im: 0.,
            })
            .collect::<Vec<_>>();
        fft.process(&mut data);

        for v in data.iter_mut() {
            *v /= (sched.len() as f64).sqrt()
        }

        PointSpread {
            unrotated: data.into(),
            rotated: OnceBox::new(),
        }
    }

    /// 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.
    ///
    /// This is useful to ignore line-broadening of the the central peak caused by the PDF.
    pub fn central_peak_radius(&self, mode: DisplayMode) -> usize {
        self.windows(2)
            .take_while(|v| mode.magnitude(v[0]) > mode.magnitude(v[1]))
            .count()
            + 1
    }

    /// Calculate the ratio between the magnitude of the central peak and the magnitude of the largest feature other than the central peak.
    ///
    /// The central peak line-broadening caused by the non-uniformity of the schedule is removed according to [`PointSpread::central_peak_radius`].
    pub fn peak_to_sidelobe_ratio(&self, mode: DisplayMode) -> f64 {
        let peak = mode.magnitude(self[0]);

        self.iter()
            .take(self.len() / 2 + 1)
            .skip(self.central_peak_radius(mode) + 1)
            .map(|v| mode.magnitude(*v))
            .max_by(|a, b| a.total_cmp(b))
            .unwrap_or(0.)
            / peak
    }
}

#[cfg(feature = "terminal-viz")]
impl crate::terminal_viz::Plot for PointSpread {
    type BarChart<'a> = <[Complex<f64>] as crate::terminal_viz::Plot>::BarChart<'a>
    where
        Self: 'a;

    type RowChart<'a> = <[Complex<f64>] as crate::terminal_viz::Plot>::RowChart<'a>
    where
        Self: 'a;

    fn bar_chart(&self, options: crate::terminal_viz::BarChartOptions) -> Self::BarChart<'_> {
        self.rotated().bar_chart(options)
    }

    fn row_chart(&self, options: crate::terminal_viz::RowChartOptions) -> Self::RowChart<'_> {
        self.rotated().row_chart(options)
    }
}