Modified Bessel Function of the First Kind Matlab & C++ Code

http://mathworld.wolfram.com/ModifiedBesselFunctionoftheFirstKind.html

Modified Bessel Function of the First Kind는 다음과 같은 수식으로 정의 된다.

$*$ { I }_{ n }(z)={ i }^{ -n }{ J }_{ n }(ix) $*$

여기에서 $*${n}$*$이 실수 $*${v}$*$인 경우엔 아래와 같이 계산할 수 있다.

$*${ I }_{ v }(z)={ \left( \cfrac { 1 }{ 2 } z \right)  }^{ v }\sum _{ k=0 }^{ \infty  }{ \frac { { \left( \cfrac { 1 }{ 2 } z \right)  }^{ 2k } }{ k!\Gamma \left( v+k+1 \right)  }  } $*$

$*$\Gamma$*$는 gamma function이다.

여기에서 다시 $*${v}$*$가 정수라 하면 다음과 같이 정리할 수 있다.

$*${ I }_{ v }(z)={ \left( \cfrac { 1 }{ 2 } z \right)  }^{ v }\sum _{ k=0 }^{ \infty  }{ \frac { { \left( \cfrac { 1 }{ 2 } z \right)  }^{ 2k } }{ k!(k+v)! }  } $*$

간단히 $*${v}=0$*$라 하면 아래와 같으며, 이것은 Kaiser Window를 생성할때 사용하는 Zeroth Order Modified Bessel Function이다.

$*${ I }_{ 0 }(z)=\sum _{ k=0 }^{ \infty  }{ { \left( \frac { { \left( \cfrac { 1 }{ 2 } z \right)  }^{ k } }{ k! }  \right)  }^{ 2 } }  $*$

Matlab 코드로 구현하려면 생각해 보아야 할 것이. 일단 저 함수는 시그마에 무한대 기호가 붙어있다.

따라서 실제로 코드로 구현하려면, 적당한 조건을 주어서 유한한 연산이 되도록 해 주어야 한다.

먼저 최대 몇번째 항 까지 쓸지 정해야 한다. 코드에 Maximum Iteration이라고 되어있고, 이것에 따라서 for-loop연산을 하게 되어있다.

이것을 정해줌으로 해서, 적어도 이 함수가 200번 연산을 하고는 끝난다는것을 보장할 수 있다.

둘째로는 각 항을 더해감에 따라 $*$k!$*$에 의해서 값이 줄어드는데, 최종 결과에 영향을 미치지 않을 만큼 충분히 작은 값이 되면 더이상 더하지 않는것이다.

이미 더해진 값을 sum이라 하고, 현재 계산된 항의 값을 inc라 했을 때, inc가 sum의 0.1% 보다 작아지면 멈춘다.

	
Matlab Code
	
function sum = besseli0(x)

M = 200; % maximum interation

sum = 0;
for m = 0 : (M-1)
   factm = factorial(m);    % m!
   inc = (1/factm*(x*0.5)^m)^2; %  (1/(m!) * (x/2)^m)^2
   
   frac = inc/sum; % increment ratio
   sum = sum + inc;
   
   if( frac <0.001) % if the increment ratio is less than 0.1% then break
       break;
   end

end



function val = factorial(m)

if( m == 0 )
    val = 1;
    return;
end

val = 1;

for idx = 2 : m
    val = val*idx;
end 

• C++ Code

mbesseli0.cpp