V_BESSELRATIO calculate the Bessel function ratio besseli(v+1,x)./besseli(v,x)
Inputs: x Bessel function argument (scalar or matrix)
v denominator Bessel function order [0]
p digits precision <=14 [5]
Outputs: y value of the Bessel function ratio besseli(v+1,x)./besseli(v,x)
Uses an algorithm from [1] of which a pseudocode version is given in [2] but which
has an misprint in the first line: "min" instead of "max".
The iteration count and minimum initial order are made dependent on the
required precision to improve efficiency.
[1] D. E. Amos. Computation of modified bessel functions and their ratios.
Mathematics of Computation, 28 (125): 239�251, jan 1974. doi: 10.1090/S0025-5718-1974-0333287-7.
[2] G. Kurz, I. Gilitschenski, and U. D. Hanebeck.
Recursive nonlinear filtering for angular data based on circular distributions.
In Proc American Control Conf, Washington, June 2013.0001 function y=v_besselratio(x,v,p) 0002 %V_BESSELRATIO calculate the Bessel function ratio besseli(v+1,x)./besseli(v,x) 0003 % 0004 % Inputs: x Bessel function argument (scalar or matrix) 0005 % v denominator Bessel function order [0] 0006 % p digits precision <=14 [5] 0007 % 0008 % Outputs: y value of the Bessel function ratio besseli(v+1,x)./besseli(v,x) 0009 % 0010 % Uses an algorithm from [1] of which a pseudocode version is given in [2] but which 0011 % has an misprint in the first line: "min" instead of "max". 0012 % The iteration count and minimum initial order are made dependent on the 0013 % required precision to improve efficiency. 0014 % 0015 % [1] D. E. Amos. Computation of modified bessel functions and their ratios. 0016 % Mathematics of Computation, 28 (125): 239�251, jan 1974. doi: 10.1090/S0025-5718-1974-0333287-7. 0017 % [2] G. Kurz, I. Gilitschenski, and U. D. Hanebeck. 0018 % Recursive nonlinear filtering for angular data based on circular distributions. 0019 % In Proc American Control Conf, Washington, June 2013. 0020 0021 % Copyright (C) Mike Brookes 2016-2017 0022 % Version: $Id: v_besselratio.m 10865 2018-09-21 17:22:45Z dmb $ 0023 % 0024 % VOICEBOX is a MATLAB toolbox for speech processing. 0025 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0026 % 0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0028 % This program is free software; you can redistribute it and/or modify 0029 % it under the terms of the GNU General Public License as published by 0030 % the Free Software Foundation; either version 2 of the License, or 0031 % (at your option) any later version. 0032 % 0033 % This program is distributed in the hope that it will be useful, 0034 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0035 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0036 % GNU General Public License for more details. 0037 % 0038 % You can obtain a copy of the GNU General Public License from 0039 % http://www.gnu.org/copyleft/gpl.html or by writing to 0040 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0041 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0042 wm=[1 1 2 4 10 20 14 21 16 22 18 23 20 25]; % lower bound for initial v as a function of p 0043 nm=[1 1 1 1 1 1 2 2 3 3 4 4 5 5]; % number of iterations as a function of p 0044 if nargin<3 || isempty(p) 0045 p=5; % default precision is 5 digits 0046 else 0047 p=min(max(ceil(p),1),14); % digits precision required 0048 end 0049 n=nm(p); % number of iterations 0050 if nargin<2 || isempty(v) 0051 v=0; % default to I1/I0 0052 end 0053 u=max(v,wm(p)); % minimum value of v for first stage 0054 s=size(x); 0055 x=x(:); 0056 r=zeros(numel(x),n+1); % intermediate estimates: one row per value of x 0057 for i=1:n+1 0058 r(:,i)=x./(u+i-0.5+sqrt((u+i+0.5)^2+x.^2)); % lower bound for ratio: Eqn (20a) from [1] 0059 end 0060 for i=1:n 0061 for k=1:n-i+1 % k is actually k+1 in (20b) of [1] 0062 r(:,k)=x./(u+k+sqrt((u+k)^2+x.^2.*r(:,k+1)./r(:,k))); % Eqn (20b) from [1] 0063 end 0064 end 0065 y=r(:,1); 0066 for i=u:-1:v+1 0067 y=1./(y+2*i./x); % backward recursion on v using Eqn (2) from [1] 0068 end 0069 y(x==0)=0; % fix up special case of x=0 0070 y=reshape(y,s); 0071 if ~nargout 0072 plot(x,y); 0073 xlabel('x'); 0074 ylabel(sprintf('I_{%d}(x)/I_{%d}(x)',v+1,v)); 0075 end