yagi 0.1.3

Batteries-included DSP library
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

use crate::error::{Error, Result};
use crate::math::gamma::{lowergammaf, lngammaf};
use crate::random::randgammaf;

pub fn randnakmf(m: f32, omega: f32) -> Result<f32> {
    // validate input
    if m < 0.5 {
        return Err(Error::Config("m cannot be less than 0.5".to_string()));
    } else if omega <= 0.0 {
        return Err(Error::Config("omega must be greater than zero".to_string()));
    }

    // generate Gamma random variable
    let alpha = m;
    let beta = omega / m;
    let x = randgammaf(alpha, beta)?;

    // sqrt(x) ~ Nakagami(m,omega)
    Ok(x.sqrt())
}

/// Nakagami-m distribution probability distribution function
/// Nakagami-m
///  f(x) = (2/Gamma(m)) (m/omega)^m x^(2m-1) exp{-(m/omega)x^2}
/// where
///      m       : shape parameter, m >= 0.5
///      omega   : spread parameter, omega > 0
///      Gamma(z): regular complete gamma function
///      x >= 0
pub fn randnakmf_pdf(x: f32, m: f32, omega: f32) -> Result<f32> {
    // validate input
    if m < 0.5 {
        return Err(Error::Config("m cannot be less than 0.5".to_string()));
    } else if omega <= 0.0 {
        return Err(Error::Config("omega must be greater than zero".to_string()));
    }

    if x <= 0.0 {
        return Ok(0.0);
    }

    let t0 = lngammaf(m);
    let t1 = m * (m / omega).ln();
    let t2 = (2.0 * m - 1.0) * x.ln();
    let t3 = -(m / omega) * x * x;

    Ok(2.0 * (-t0 + t1 + t2 + t3).exp())
}

// Nakagami-m distribution cumulative distribution function
//  F(x) = gamma(m, x^2 m / omega) / Gamma(m)
//  where
//      gamma(z,a) = lower incomplete gamma function
//      Gamma(z)   = regular gamma function
//
pub fn randnakmf_cdf(x: f32, m: f32, omega: f32) -> Result<f32> {
    // validate input
    if m < 0.5 {
        return Err(Error::Config("m cannot be less than 0.5".to_string()));
    } else if omega <= 0.0 {
        return Err(Error::Config("omega must be greater than zero".to_string()));
    }

    if x <= 0.0 {
        return Ok(0.0);
    }

    let t0 = lowergammaf(m, x * x * m / omega);
    let t1 = lngammaf(m);
    Ok((t0 - t1).exp())
}