comnum_waveform/cpm/
pulse.rs

1//! The frequency pulses of CPM, the closed forms of the Gaussian pulse of GMSK and the
2//! power spectral density of MSK.
3
4use comnum_core::{Error, Result};
5use comnum_math::{normal_pdf, qfunc};
6use std::f64::consts::{PI, TAU};
7
8/// The frequency pulse $g(t)$ of a CPM signal (see the module documentation).
9#[non_exhaustive]
10#[derive(Clone, Copy, Debug, PartialEq)]
11pub enum FrequencyPulse {
12    /// $L$REC: constant over $L$ symbols (1REC is CPFSK; with $h = 1/2$, MSK).
13    Rectangular,
14    /// $L$RC: a raised cosine over $L$ symbols.
15    RaisedCosine,
16    /// GMSK: the Gaussian pulse of Murota & Hirade with the normalised bandwidth `bt`,
17    /// truncated to $L$ symbols.
18    Gaussian {
19        /// The 3 dB bandwidth of the Gaussian filter times the symbol duration, $BT$.
20        bt: f64,
21    },
22}
23
24impl FrequencyPulse {
25    /// The short name of the pulse, `rec`, `rc` or `gmsk` (read back by
26    /// [`FrequencyPulse::from_name`]).
27    pub fn name(&self) -> &'static str {
28        match self {
29            FrequencyPulse::Rectangular => "rec",
30            FrequencyPulse::RaisedCosine => "rc",
31            FrequencyPulse::Gaussian { .. } => "gmsk",
32        }
33    }
34
35    /// The phase response $q(t)$ of the pulse of `length` symbols at the time `t` in
36    /// symbol durations: 0 before 0, $1/2$ after `length`.
37    ///
38    /// # Examples
39    /// ```
40    /// use comnum_waveform::cpm::FrequencyPulse;
41    /// assert_eq!(FrequencyPulse::Rectangular.phase_response(0.5, 1), 0.25);
42    /// assert_eq!(FrequencyPulse::RaisedCosine.phase_response(3.0, 3), 0.5);
43    /// let g = FrequencyPulse::Gaussian { bt: 0.3 };
44    /// assert!((g.phase_response(1.5, 3) - 0.25).abs() < 1e-12);     // symmetric pulse
45    /// ```
46    pub fn phase_response(&self, t: f64, length: usize) -> f64 {
47        let l = length as f64;
48        if t <= 0.0 {
49            return 0.0;
50        }
51        if t >= l {
52            return 0.5;
53        }
54        match *self {
55            FrequencyPulse::Rectangular => t / (2.0 * l),
56            FrequencyPulse::RaisedCosine => t / (2.0 * l) - (TAU * t / l).sin() / (2.0 * TAU),
57            FrequencyPulse::Gaussian { bt } => {
58                // the pulse centred at L/2, truncated to [0, L] and scaled to the area 1/2
59                let lo = gmsk_phase(-l / 2.0, bt);
60                let hi = gmsk_phase(l / 2.0, bt);
61                0.5 * (gmsk_phase(t - l / 2.0, bt) - lo) / (hi - lo)
62            }
63        }
64    }
65
66    /// Parse `rec` (`rectangular`), `rc` (`raised-cosine`) or `gmsk` (`gaussian`) with
67    /// the bandwidth `bt` for the Gaussian pulse.
68    ///
69    /// # Errors
70    /// [`Error::Config`] for another name.
71    pub fn from_name(name: &str, bt: f64) -> Result<FrequencyPulse> {
72        match name.to_ascii_lowercase().as_str() {
73            "rec" | "rectangular" => Ok(FrequencyPulse::Rectangular),
74            "rc" | "raised-cosine" => Ok(FrequencyPulse::RaisedCosine),
75            "gmsk" | "gaussian" => Ok(FrequencyPulse::Gaussian { bt }),
76            other => Err(Error::Config(format!("pulse must be rec, rc or gmsk, got {other:?}"))),
77        }
78    }
79}
80
81/// $\int Q(x) \\, dx = x Q(x) - \phi(x)$.
82fn q_integral(x: f64) -> f64 {
83    x * qfunc(x) - normal_pdf(x)
84}
85
86/// The Gaussian frequency pulse of GMSK, centred at 0, with $T = 1$ (Murota & Hirade,
87/// 1981): the response of the Gaussian low-pass filter of 3 dB bandwidth $B$ to a
88/// rectangular pulse of one symbol,
89///
90/// $$
91/// g(t) = \frac{1}{2}\left[Q\left(\frac{2 \pi B}{\sqrt{\ln 2}}\left(t - \frac{1}{2}\right)\right) -
92/// Q\left(\frac{2 \pi B}{\sqrt{\ln 2}}\left(t + \frac{1}{2}\right)\right)\right], \qquad
93/// \int_{-\infty}^{\infty} g(t) \\, dt = \frac{1}{2}.
94/// $$
95///
96/// # Examples
97/// ```
98/// use comnum_waveform::cpm::gmsk_pulse;
99/// // the pulse is symmetric and its area is 1/2 (Simpson's rule over ±4 symbols)
100/// let area = comnum_math::simpson(|t| gmsk_pulse(t, 0.3), -4.0, 4.0, 800);
101/// assert!((area - 0.5).abs() < 1e-9);
102/// assert!((gmsk_pulse(0.7, 0.3) - gmsk_pulse(-0.7, 0.3)).abs() < 1e-15);
103/// ```
104pub fn gmsk_pulse(t: f64, bt: f64) -> f64 {
105    let kappa = TAU * bt / 2f64.ln().sqrt();
106    0.5 * (qfunc(kappa * (t - 0.5)) - qfunc(kappa * (t + 0.5)))
107}
108
109/// The phase response of the Gaussian pulse, $q(t) = \int_{-\infty}^{t} g(\tau) \\, d\tau$
110/// in closed form (0 at $-\infty$, $1/4$ at 0, $1/2$ at $+\infty$).
111///
112/// # Examples
113/// ```
114/// use comnum_waveform::cpm::gmsk_phase;
115/// assert!((gmsk_phase(0.0, 0.3) - 0.25).abs() < 1e-12);
116/// assert!((gmsk_phase(10.0, 0.3) - 0.5).abs() < 1e-12);
117/// ```
118pub fn gmsk_phase(t: f64, bt: f64) -> f64 {
119    let kappa = TAU * bt / 2f64.ln().sqrt();
120    0.5 + (q_integral(kappa * (t - 0.5)) - q_integral(kappa * (t + 0.5))) / (2.0 * kappa)
121}
122
123/// The power spectral density of MSK at the normalised frequency `ft` $= f T$, with
124/// $T = 1$ and a unit power (Proakis & Salehi, §3.4-3):
125///
126/// $$
127/// S(f) = \frac{16 T}{\pi^{2}} \left[\frac{\cos 2 \pi f T}{1 - 16 f^{2} T^{2}}\right]^{2},
128/// $$
129///
130/// finite at $f T = \pm 1/4$ (where it equals $T$), zero at $f T = \pm 3/4, \pm 5/4, \dots$
131///
132/// # Examples
133/// ```
134/// use comnum_waveform::cpm::msk_psd;
135/// assert!((msk_psd(0.0) - 16.0 / std::f64::consts::PI.powi(2)).abs() < 1e-12);
136/// assert!((msk_psd(0.25) - 1.0).abs() < 1e-6);
137/// assert!(msk_psd(0.75).abs() < 1e-12);
138/// ```
139pub fn msk_psd(ft: f64) -> f64 {
140    let d = 1.0 - 16.0 * ft * ft;
141    let ratio = if d.abs() < 1e-7 { PI / 4.0 } else { (TAU * ft).cos() / d };
142    16.0 / (PI * PI) * ratio * ratio
143}
comnum_waveform/cpm/pulse.rs

comnum_waveform/cpm/
pulse.rs
