Mirror/aicodix___dsp
1
0
Fork 0
mirror of https://github.com/aicodix/dsp.git synced 2026-04-27 14:30:36 +00:00
Code Issues Projects Releases Packages Wiki Activity

simplified i0() and added some LaTeX code

Browse source
This commit is contained in:
Ahmet Inan 2018-03-03 19:49:56 +01:00
commit 2d4e1886c3
1 changed files with 10 additions and 2 deletions
Download patch file Download diff file Expand all files Collapse all files

12
window.hh
View file

@ -89,6 +89,14 @@ template <int TAPS, typename TYPE>
class Kaiser
{
TYPE w[TAPS];
/*
i0() implements the zero-th order modified Bessel function of the first kind:
https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1
$I_\alpha(x) = i^{-\alpha} J_\alpha(ix) = \sum_{m=0}^\infty \frac{1}{m!\, \Gamma(m+\alpha+1)}\left(\frac{x}{2}\right)^{2m+\alpha}$
$I_0(x) = J_0(ix) = \sum_{m=0}^\infty \frac{1}{m!\, \Gamma(m+1)}\left(\frac{x}{2}\right)^{2m} = \sum_{m=0}^\infty \left(\frac{x^m}{2^m\,m!}\right)^{2}$
We obviously can't use the factorial here, so let's get rid of it:
$= 1 + \left(\frac{x}{2 \cdot 1}\right)^2 + \left(\frac{x}{2 \cdot 1}\cdot \frac{x}{2 \cdot 2}\right)^2 + \left(\frac{x}{2 \cdot 1}\cdot \frac{x}{2 \cdot 2}\cdot \frac{x}{2 \cdot 3}\right)^2 + .. = 1 + \sum_{m=1}^\infty \left(\prod_{n=1}^m \frac{x}{2n}\right)^2$
*/
TYPE i0(TYPE x)
{
Kahan<TYPE> sum(1.0);
@ -97,8 +105,8 @@ class Kaiser
// float: 25 iterations
// double: 35 iterations
for (int n = 1; n < 35; ++n) {
TYPE tmp = x / TYPE(2 * n);
if (sum.same(val *= tmp * tmp))
val *= x / TYPE(2 * n);
if (sum.same(val * val))
return sum();
}
return sum();