Mathematical Operations

The math module implements several useful functions for digital signal processing including transcendental function not necessarily in the standard C library, windowing functions, and polynomial manipulation methods.

Transcendental Functions∞

This section describes the implementation and interface to transcendental functions not in the C standard library including a full arrangement of Gamma and Bessel functions.[ref:tab-math-transcendentals] summarizes the interfaces provided in liquid .

Table [tab-math-transcendentals]. Summary of Transcendental Math Interfaces

Function	Interface
\(\ln\Gamma(z)\)	liquid_lngammaf(z)
\( \Gamma(z)\)	liquid_gammaf(z)
\(\ln\gamma(z,\alpha)\)	liquid_lnlowergammaf(z,alpha)
\( \gamma(z,\alpha)\)	liquid_lowergammaf(z,alpha)
\(\ln\Gamma(z,\alpha)\)	liquid_lnuppergammaf(z,alpha)
\( \Gamma(z,\alpha)\)	liquid_uppergammaf(z,alpha)
\(n!\)	liquid_factorialf(n)
\(\ln I_\nu(z)\)	liquid_lnbesselif(nu,z)
\( I_\nu(z)\)	liquid_besselif(nu,z)
\( I_0 (z)\)	liquid_besseli0f(z)
\( J_\nu(z)\)	liquid_besseljf(nu,z)
\( J_0 (z)\)	liquid_besselj0f(z)
\(Q(z)\)	liquid_Qf(z)
\(Q_M(\alpha,\beta)\)	liquid_MarcumQf(M,alpha,beta)
\(Q_1(\alpha,\beta)\)	liquid_MarcumQ1f(alpha,beta)
\(\sinc(z)\)	liquid_sincf(z)
\(\lceil\log_2(n)\rceil\)	liquid_nextpow2(n)
\({n \choose k}\)	liquid_nchoosek(n,k)

liquid_gammaf(z), liquid_lngammaf(z)∞

liquid computes \(\Gamma(z)\) from \(\ln\Gamma(z)\) (see below) due to its steep, exponential response to \(z\) . The complete Gamma function is defined as

$$ \Gamma(z) \triangleq \int_0^\infty{t^{z-1}e^{-t}dt} $$

The upper an lower incomplete Gamma functions are described in Sections[ref:section-math-transcendentals-uppergamma] and[ref:section-math-transcendentals-lowergamma] , respectively. The natural log of the complete Gamma function is computed by splitting into discrete piecewise sections:

$$ \ln\left[ \Gamma(z) \right] \approx \begin{cases} % undefined for z \lt 0 \text{undefined} & z \lt 0 \\ % low signal approximation (recursion) \ln\Gamma(z+1) - \ln(z) & 0 \le z \lt 10 \\ % high signal approximation \frac{z}{2} \ln\left( \frac{2\pi}{z} \right) \left( \ln\left(z + \frac{1}{12 z - 0.1/z} \right) - 1 \right) & z \ge 0.6 \end{cases} $$

liquid_lowergammaf(z,a), liquid_lnlowergammaf(z,a) (lower incomplete Gamma)∞

Like \(\Gamma(z)\) , liquid computes the lower incomplete gamma function\(\gamma(z,\alpha)\) from its logarithm \(\ln\gamma(z,\alpha)\) due to its steep, exponential response to \(z\) . The lower incomplete Gamma function is defined as

$$ \gamma(z,\alpha) \triangleq \int_0^\alpha{ t^{z-1}e^{-t}dt } $$

liquid computes the log of lower incomplete Gamma function as

$$ \ln\gamma(z,\alpha) = z \ln(\alpha) + \ln\Gamma(z) - \alpha + \ln\left[ \sum_{k=0}^\infty{ \frac{\alpha^k}{\Gamma(z + k + 1)} } \right] $$

liquid_uppergammaf(z,a), liquid_lnuppergammaf(z,a) (upper incomplete Gamma)∞

Like \(\Gamma(z)\) , liquid computes the upper incomplete gamma function\(\Gamma(z,\alpha)\) from \(\ln\Gamma(z,\alpha)\) due to its steep, exponential response to \(z\) . The complete Gamma function is defined as

$$ \Gamma(z,\alpha) \triangleq \int_\alpha^\infty{t^{z-1}e^{-t}dt} $$

By definition the sum of the lower and upper incomplete gamma functions is the complete Gamma function:\(\Gamma(z) = \gamma(z,\alpha) + \Gamma(z,\alpha)\) . As such, liquid computes the upper incomplete Gamma function as

$$ \Gamma(z,\alpha) = \Gamma(z) - \gamma(z,\alpha) $$

liquid_factorialf(n)∞

liquid computes \(n!=n\cdot(n-1)\cdot(n-2)\cdots3\cdot2\cdot1\) iteratively for small values of \(n\) , and with the Gamma function for larger values. Specifically, \(n! = \Gamma(n+1)\) .

liquid_nchoosek()∞

liquid computes binomial coefficients using the liquid_nchoosek() method:

$$ {n \choose k} = \frac{n!}{(n-k)!k!} $$

Because the arguments can explode for relatively large values of \(n\) and\(k\) , liquid uses the following approximation under certain conditions:

$$ {n \choose k} \approx \exp\Bigl\{ \ln\Gamma(n+1) - \ln\Gamma(n-k+1) - \ln\Gamma(k+1) \Bigr\} $$

liquid_nextpow2()∞

computes \(\lceil \log_2(x) \rceil\)

liquid_sinc(z)∞

The \(\sinc\) function is defined as

$$ \sinc(z) = \frac{\sin(\pi z)}{\pi z} $$

Simply evaluating the above equation with finite precision for \(z\) results in a discontinuity for small \(z\) , and is approximated by expanding the first few terms of the series

$$ \sinc(z) = \prod_{k=1}^{\infty}{ \cos\left( 2^{-k} \pi z \right) } $$

liquid_lnbesselif(), liquid_besselif(), liquid_besseli0f()∞

\(I_\nu(z)\) is the modified Bessel function of the first kind and is particularly useful for filter design. An iterative method for computing \(I_\nu\) comes from Gross(1995),

$$ I_\nu(z) = \left(\frac{z}{2}\right)^\nu \sum_{k=0}^{\infty}{\frac{\left(\frac{1}{4}z^2\right)^k}{k!\Gamma(k+\nu+1)}} $$

Due to its steep response to \(z\) it is often useful to compute\(I_\nu(z)\) by first computing \(\ln I_\nu(z)\) as

..eqnarray [eqn-math-lnbesseli]

\ln I_\nu(z) & = &
    \nu\ln(z/2) + \ln\left[
        \sum_{k=0}^\infty{
            \frac{
                \left(\frac{1}{4}z^2\right)^k
            }{
                k! \Gamma(\nu + k + 1)
            }
        }
    \right] \\
    & = &
    \nu\ln(z/2) + \ln\left[
        \sum_{k=0}^\infty{ \exp\Bigl\{
            2 k \ln(z/2) - \ln\Gamma(k+1) - \ln\Gamma(\nu+k+1)
        \Bigr\} }
    \right] \\

For \(\nu=0\) a good approximation can be derived by using piecewise polynomials,

$$ \ln\Bigl[\ln\left(I_0(z)\right)\Bigr] \approx c_0 + c_1 t + c_2 t ^2 + c_3 t^3 $$

where \(t=\ln(z)\) and

$$ \left\{c_0,c_1,c_2,c_3\right\} = \begin{cases} \left\{\text{-1.52624, 1.9597, -9.4287e-03, -7.0471e-04}\right\} & t \lt 0.5 \\ \left\{\text{-1.5531, 1.8936, -0.07972, -0.01333}\right\} & 0.5 \le t \lt 2.3 \\ \left\{\text{-1.2958, 1.7693, -0.1175, 0.006341}\right\} & \text{else}. \end{cases} $$

This is a particularly useful approximation for the Kaiser window in fixed-point math where \(w[n]\) is computed as the ratio of two large numbers.

liquid_lnbesseljf(), liquid_besselj0f()∞

\(J_\nu(z)\) is the Bessel function of the first kind and is found in Doppler filter design. liquid computes \(J_\nu(z)\) using the series expansion

$$ J_\nu(z) = \sum_{k=0}^\infty{ \frac{ (-1)^k }{ 2^{2k+|v|} k! \left(|v|+k\right)! } z^{2k+|v|} } $$

liquid_Qf(), liquid_MarcumQf(), liquid_MarcumQ1f()∞

The \(Q\) -function is commonly used in signal processing and is defined as

.. eqnarray
Q(z) &=& \frac{1}{2}\left(1 - \textup{erf}(z/\sqrt{2})\right) \\
     &=& \frac{1}{\sqrt{2\pi}} \int_{z}^{\infty} { \exp\left\{-u^2/2 \right\} du }

Similarly Marcum's \(Q\) -function is defined as the following, with an appropriate expansion:

.. eqnarray

Q_M(\alpha,\beta) & = &
    \int_{\beta}^{\infty}{
        u\left(\frac{u}{\alpha}\right)^{M-1}
        \exp\left\{ -\frac{u^2 + \alpha^2}{2}\right\}
        I_{M-1}(\alpha u)
        du
    }\\
    & = &
    \exp\left\{-\frac{\alpha^2 + \beta^2}{2}\right\}
    \sum_{k=1-M}^{\infty}{
        \left(\frac{\alpha}{\beta}\right)^k
        I_k(\alpha\beta)
    }

where \(I_\nu\) is the modified Bessel function of the first kind (see [ref:section-math-transcendentals-besseli] ). liquid implements \(Q_M(\alpha,\beta)\) with the function liquid_MarcumQf(M,a,b) using the approximation {cite:Helstrom:1992((25))}

$$ Q_M(\alpha,\beta) \approx \textup{erfc}(u), \,\,\, u=\frac{\beta-\alpha-M}{\sigma^2}, \,\,\, \sigma = M + 2\alpha $$

which works over a reasonable range of \(M\) , \(\alpha\) , and \(\beta\) . The special case for \(M=1\) is implemented in liquid using the function liquid_MarcumQ1f(M,a,b) using the expansion {cite:Helstrom:1960},

$$ Q_1(\alpha,\beta) = \exp\left\{-\frac{\alpha^2 + \beta^2}{2}\right\} \sum_{k=0}^{\infty}{ \left(\frac{\alpha}{\beta}\right)^k I_k(\alpha\beta) } $$

which converges quickly with a few iterations.

Complex Trigonometry∞

This section describes the implementation and interface to complex trigonometric functions not in the C standard library.[ref:tab-math-transcendentals] summarizes the interfaces provided in liquid .

Table [tab-math-complex]. Summary of Complex Trigonometric Math Interfaces

Function	Interface
\(\sqrt{z}\)	liquid_csqrtf(z)
\(e^{z}\)	liquid_cexpf(z)
\(\ln(z)\)	liquid_clogf(z)
\(\sin^{-1}(z)\)	liquid_casinf(z)
\(\cos^{-1}(z)\)	liquid_cacosf(z)
\(\tan^{-1}(z)\)	liquid_catanf(z)

liquid_csqrtf()∞

The function liquid_csqrtf(z) computes the complex square root of a number

$$ \sqrt{z} = \sqrt{\frac{r+a}{2}} + j\text{sgn}\bigl(\Im\{z\}\bigr) \sqrt{\frac{r-a}{2}} $$

where \(r=|z|\) , \(a=\Re\{z\}\) , and \(\text{sgn}(t)=t/|t|\) .

liquid_cexpf()∞

The function liquid_cexpf(z) computes the complex exponential of a number

$$ e^{z} = \exp\bigl\{a\bigr\} \bigl( \cos(b) + j\sin(b) \bigr) $$

where \(a=\Re\{z\}\) and \(b=\Im\{z\}\) .

liquid_clogf()∞

The function liquid_clogf(z) computes the complex natural logarithm of a number.

$$ \log(z) = \log(|z|) + j\arg(z) $$

liquid_cacosf()∞

The function liquid_cacosf(z) computes the complex \(\arccos\) of a number

$$ \arccos(z) = \begin{cases} -j \log\bigl( z + \sqrt{z^2 - 1} \bigr) & \text{sgn}\bigl(\Re\{z\}\bigr) = \text{sgn}\bigl(\Im\{z\}\bigr) \\ -j \log\bigl( z - \sqrt{z^2 - 1} \bigr) & \text{otherwise} \end{cases} $$

liquid_casinf()∞

The function liquid_casinf(z) computes the complex \(\arcsin\) of a number

$$ \arcsin(z) = \frac{\pi}{2} - \arccos(z) $$

liquid_catanf()∞

The function liquid_catanf(z) computes the complex \(\arctan\) of a number

$$ \arctan(z) = \frac{j}{2} \log\left( \frac{1-jz}{1+jz} \right) $$

Modular Arithmetic∞

liquid_is_prime(n)∞

Returns 1 if \(n\) is prime, 0 otherwise.

liquid_factor(n,*factors,*num_factors)∞

Computes all the prime factors of a number \(n\) in ascending order. Example:

unsigned int factors[LIQUID_MAX_FACTORS];
unsigned int num_factors;
liquid_factor(280, factors, &num_factors);

>> num_factors = 5>> factors = { 2, 2, 2, 5, 7 }

liquid_unique_factor(n,*factors,*num_factors)∞

Computes all the unique prime factors of a number \(n\) in ascending order. Example:

unsigned int factors[LIQUID_MAX_FACTORS];
unsigned int num_factors;
liquid_unique_factor(280, factors, &num_factors);

>> num_factors = 3>> factors = { 2, 5, 7 }

liquid_modpow(base,exp,n)∞

Computes \(c = b^x \pmod n\)

liquid_primitive_root(n)∞

Finds and returns smallest primitive root of a number \(n\) .

liquid_primitive_root_prime(n)∞

Finds and returns smallest primitive root of a prime number \(n\) .

liquid_totient(n)∞

Euler's totient function \(\varphi(n)\) computes the number of positive integers less than or equal to \(n\) that are relatively prime to \(n\) . The totient function is computed as

$$ \varphi(n) = n \prod_{p|n} {\left(1 - \frac{1}{p}\right)} $$

where the product is computed over the unique factors of \(n\) . For example, if \(n=24 = 2^3 \cdot 3\) , then\(\varphi(24) = 24\left(1-\frac{1}{2}\right)\left(1-\frac{1}{3}\right)=8\) .