Modified Bessel Functions

Definition
Modified Bessel Functions are a special form of Bessel Functions, which were defined by Friedrich Bessel and Daniel Bernoulli as solutions, y(x), to Bessel's Differential Equation:


 * $$x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - v^2)y = 0$$

In this equation v is an unrestricted constant which can be real or complex, and defines the "order" of the function. A special case of the Bessel Function occurs when v is an integer. Like regular Bessel Functions, Modified Bessel Functions appear in various applications in physics and engineering, namely in dealing with stochastic processes such as noise in communication systems.

The nth Order Modified Bessel Function of the First Kind is defined as:


 * $$ I_n(x) = \frac{1}{\pi}\int_0^{\pi} e^{x\cos{\theta}}\cos{(n\theta)}  d\theta$$

Modified Bessel Functions, unlike regular Bessel Functions which behave like decaying sines or cosines, are monotonically increasing or decreasing functions that are defined as 1 for x=0.

Derivation
The Modified Bessel Equation is derived by allowing the argument, x, to be purely imaginary. Thus, by letting x = ix, the equation becomes:  : 


 * $$ (ix)^2 \frac{d^2 y}{d(ix)^2} + ix \frac{dy}{d(ix)} + ((ix)^2 - v^2)y = 0 $$


 * $$ -x^2 \frac{d^2 y}{-dx^2} + ix \frac{dy}{idx} + (-x^2 - v^2)y = 0 $$


 * $$ x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} - (x^2 + v^2)y = 0 $$



This differential equation is then solved by:


 * $$ y(x) = AJ_v(ix) + BJ_{-v}(ix) $$


 * $$ y(x) = AJ_v(ix) + BY_v(ix) $$

where $$J_v(x)$$ is a Bessel Function of the First Kind, and $$Y_v(x)$$ is a Bessel Function of the Second Kind.

If v is not an integer, and by substituting $$ ix = xe^{i\frac{\pi}{2}}$$, the Modified Bessel Function of the First Kind is defined as:


 * $$ I_v(x) = J_v(xe^{i\frac{\pi}{2}})e^{-i\frac{\pi}{2}v}$$


 * $$ y(x) = AI_v(x) + BI_{-v}(x) $$

In the case where v is a positive integer, i.e. v=n, then the Modified Bessel Function of the First Kind of nth order can be defined as:


 * $$ I_n(x) = (i)^{-n}J_n(ix)$$

where, as with regular Bessel Functions, it can be shown that


 * $$ I_n(x) = I_{-n}(x)$$


 * $$ y(x) = AI_n(x) + BK_n(x) $$

and $$K_n(x)$$ is defined as the Modified Bessel Function of the Second Kind, of nth order.

When evaluated for real positive values of x, the Modified Bessel Function of the First Kind can be defined as:


 * $$ I_n(x) = \frac{1}{\pi}\int_0^{\pi} e^{x\cos{\theta}}\cos{(n\theta)}  d\theta$$

Approximations for Modified Bessel Functions of the First Kind
In many applications, such as dealing with white nosie, it is important to be able to accurately evaluate Modified Bessel Function for various positive arguments, x. Two methods for efficient approximation of different ranges of x will be discussed in this article.

Power Series Expansion
The Power Series Expansion for Modified Bessel Functions of the First Kind for an integer order, n is given by


 * $$ I_n(x) = \sum_{m=0}^\infty \frac{x^{2m+n}}{ 2^{2m+n} \cdot m! (m+n)!} $$

The first terms of this sum can be used to accurately evaluate Modified Bessel Functions for an argument x. The more terms that are used, the more accurate the approximation.

This is an appropriate approximation for small values of x, but as the argument gets larger the numerator becomes significantly larger

Zero Order

 * $$ I_0(x) = \sum_{m=0}^\infty \frac{x^{2m}}{ 2^{2m}\cdot (m!)^2} $$


 * $$ I_0(x) = 1 + \frac{x^2}{ 2^2(1!)^2} + \frac{x^4}{2^4(2!)^2 } + \frac{x^6}{2^6(3!)^2} + ...... $$

Using the first five terms for the approximation gives:


 * $$ I_0(x) \approx 1 + \frac{x^2}{4} + \frac{x^4}{64} + \frac{x^6}{2304} + \frac{x^8}{147456}  $$

Example
Evaluating the five term approximation for $$ I_0(x)$$ for a small value, x=2, yields 2.279513889. The actual value of the function 2.279585302, which means there is a 0.007% error.

For x=5, this approximation gives 26.44640435, where the actual value is 27.23987182, an error of 2.9%. Adding one more term to the sum yields 27.1086782, and an error of 0.48%.

First Order

 * $$ I_1(x) = \sum_{m=0}^\infty \frac{x^{2m+1}}{ 2^{2m+1} m! (m+1)!} $$


 * $$ I_1(x) = \frac{x}{2} + \frac{x^3}{ 2^3\cdot 1!\cdot2!} + \frac{x^5}{2^5\cdot2!\cdot3! } + \frac{x^7}{2^7\cdot3!\cdot4!} + ...... $$

Once again using the first five terms for the approximation, this becomes:


 * $$ I_1(x) \approx \frac{x}{2} + \frac{x^3}{16} + \frac{x^5}{384} + \frac{x^7}{18432} + \frac{x^9}{1474560}  $$

Higher Orders
The Power Series expansion is an efficient way for approximating small values of x with a degree of error that can be minimized by adding more terms from the sum. This becomes impractical for large values of x as more terms are needed to keep the error within a tolerable range. Thus, this method should only be used for values of x less than or equal to 5.

Asymptotic Expansion
In order to compute values of x larger than 5, the Modified Bessel Function of the First Kind is defined by its Asymptotic Expansion, as follows for zero and first order.

Zero Order

 * $$ I_0(x) \approx \frac{e^x}{\sqrt{2\pi x}} (1 + \frac{1^2}{1!\cdot8x} + \frac{1^2 3^2}{2!(8x)^2} + \frac{1^2 3^2 5^2}{3!(8x)^3} + \frac{1^2 3^2 5^2 7^2}{4!(8x)^4} + ... )  $$

Using only the first five terms of the inner sum gives:


 * $$ I_0(x) \approx \frac{e^x}{\sqrt{2\pi x}} (1 + \frac{1}{8x} + \frac{9}{128x^2} + \frac{225}{3072x^3} + \frac{11025}{98304x^4}) $$

Example
Evaluating the five term approximation for $$ I_0(x)$$ for x=20 gives 43,558,279. The value given in tables is 43,558,283, which gives an error of 0.0009%.

For x=100 this yields 1.07374171 x 10^42, which is the approximate value given in tables.

First Order

 * $$ I_1(x) \approx \frac{e^x}{\sqrt{2\pi x}} (1 - \frac{1 \cdot 3}{1!\cdot 8x} - \frac{1^2\cdot3\cdot 5}{2!(8x)^2} - \frac{1^2\cdot 3^2 \cdot 5\cdot 7}{3!(8x)^3} - ... )$$

Using only the first five terms of the inner sum gives:


 * $$ I_1(x) \approx \frac{e^x}{\sqrt{2\pi x}} (1 - \frac{3}{8x} - \frac{15}{128x^2} - \frac{315}{3072x^3} - \frac{14175}{98304x^4}) $$

When x becomes very large (larger than 150), the approximation is simplified to:


 * $$ I_1(x) \approx \frac{e^x}{\sqrt{2\pi x}} $$

Higher Orders
For any higher order Modified Bessel Functions, the general approximation given can be used so long as x is much, much greater than n.


 * $$ I_n(x) \approx \frac{e^x}{\sqrt{2\pi x}} $$

More accurate calculation for any order of Modified Bessel Functions can be performed at High Precision Calculation of Modified Bessel Function

Some Applications
There are many applications which require accurate computation of the Modified Bessel Function. Two of these applications are in Radar Detection Theory and Windowing in Digital Signal Processing.

Radar Detection Theory
In Radar Detection Theory, the returned target signal is defined as v(t) = s(t) + n(t), where s(t) is the return signal from a target and n(t) is additive zero mean white Gaussian noise with variance $$\sigma^2$$. The function v(t) is defined as the output of an IF Band Pass Filter after being received by the antenna and passing through a low noise amplifier. This function is then passed through an envelope detector and a low pass filter which produces r(t), the envelope of v(t).  : 

The probability density function of the random variables r(t) and $$\phi$$(t), the random noise phase of v(t), is given by:


 * $$ f(r,\phi<>) = \frac{r}{2\pi \sigma^2} e^{ -\frac{r^2 + A^2 }{2 \sigma^2} } e^{ \frac{rA\cos \phi}{\sigma^2} }$$

where A is the amplitude of the sine wave from the target return signal.

Then, dropping the phase dependence by integrating from 0 to 2pi over d$$\phi$$, the function becomes:


 * $$ f(r) = \int_0^{2\pi} f(r,\phi) d\phi = \frac{r}{\sigma^2} e^{-\frac{r^2+A^2}{2 \sigma^2} } \frac{1}{2 \pi} \int_0^{2\pi} e^{ \frac{rA \cos \phi}{\sigma^2} } d\phi$$

where the integral on the right hand side is a zero-order Modified Bessel Function. Thus, f(r) can be rewritten as


 * $$ f(r) = \frac{r}{\sigma^2} e^{-\frac{r^2+A^2}{2 \sigma^2} } I_0(\frac{rA}{\sigma^2} $$

which follows a Ricean probability density function for typical values of $$\frac{rA}{\sigma^2}$$

If there is no target return signal, then the argument of the Modified Bessel Function is 0, which evaluates to 1, and the function simplifies to


 * $$ f(r) = \frac{r}{\sigma^2} e^{-\frac{r^2+A^2}{2 \sigma^2} } $$

Thus, in the absence of a target return signal, the probability density function follows a Rayleigh distribution.

When there is a strong target signal return, then $$\frac{rA}{\sigma^2}$$ is very large, and the Modified Bessel Function can be approximated by $$ \frac{ e^{ \frac{rA}{\sigma^2} } }{ \sqrt{ \frac{2 \pi rA}{\sigma^2} } }$$. After simplification, it is seen that when there is a strong target return the function follows a Gaussian probability density function.


 * $$ f(r) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(r-A)^2}{2 \sigma^2} } $$

Kaiser Window
The Kaiser window is a one-parameter family of window functions used for Digital Signal Processing, and is defined by the formula

 : 


 * $$ w_n = \frac{I_0( \alpha \sqrt{ 1 - (\frac{n}{m})^2 )} } {I_0(\alpha)} $$


 * $$ for: -M \leq n \leq M $$

where:
 * I0 is the zeroth order Modified Bessel function of the first kind.
 * α is an arbitrary real number that determines the shape of the window;
 * in the frequency domain, it determines the trade-off between main-lobe width and side lobe level, which is a central decision in window design.
 * M is an integer, and the length of the sequence is N = M+1.

When N is an odd number, the peak value of the window is wM/2 = 1. And when N is even, the peak values are wN/2-1 = wN/2 &lt; 1.

This function requires accurate calculation of Modified Bessel Functions