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.
The inverse of r=v_besselratio(k,0) is available as the function k=v_besratinv0(r).
[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.
Revision History:
2018/09/21 Original Version
2022/04/21 Corrected error for x=Inf. Revised comments.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 has an misprint 0011 % in the first line: "min" instead of "max". The iteration count and minimum initial order are made 0012 % dependent on the required precision to improve efficiency. 0013 % The inverse of r=v_besselratio(k,0) is available as the function k=v_besratinv0(r). 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 % Revision History: 0022 % 2018/09/21 Original Version 0023 % 2022/04/21 Corrected error for x=Inf. Revised comments. 0024 0025 % Copyright (C) Mike Brookes 2016-2017 0026 % Version: $Id: v_besselratio.m 10865 2018-09-21 17:22:45Z dmb $ 0027 % 0028 % VOICEBOX is a MATLAB toolbox for speech processing. 0029 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0030 % 0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0032 % This program is free software; you can redistribute it and/or modify 0033 % it under the terms of the GNU General Public License as published by 0034 % the Free Software Foundation; either version 2 of the License, or 0035 % (at your option) any later version. 0036 % 0037 % This program is distributed in the hope that it will be useful, 0038 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0039 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0040 % GNU General Public License for more details. 0041 % 0042 % You can obtain a copy of the GNU General Public License from 0043 % http://www.gnu.org/copyleft/gpl.html or by writing to 0044 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0046 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 0047 nm=[1 1 1 1 1 1 2 2 3 3 4 4 5 5]; % number of iterations as a function of p 0048 if nargin<3 || isempty(p) 0049 p=5; % default precision is 5 digits 0050 else 0051 p=min(max(ceil(p),1),14); % digits precision required 0052 end 0053 n=nm(p); % number of iterations 0054 if nargin<2 || isempty(v) 0055 v=0; % default to I1/I0 0056 end 0057 u=max(v,wm(p)); % minimum value of v for first stage 0058 s=size(x); 0059 x=x(:); 0060 r=zeros(numel(x),n+1); % intermediate estimates: one row per value of x 0061 for i=1:n+1 0062 r(:,i)=x./(u+i-0.5+sqrt((u+i+0.5)^2+x.^2)); % lower bound for ratio: Eqn (20a) from [1] 0063 end 0064 for i=1:n 0065 for k=1:n-i+1 % k is actually k+1 in (20b) of [1] 0066 r(:,k)=x./(u+k+sqrt((u+k)^2+x.^2.*r(:,k+1)./r(:,k))); % Eqn (20b) from [1] 0067 end 0068 end 0069 y=r(:,1); 0070 for i=u:-1:v+1 0071 y=1./(y+2*i./x); % backward recursion on v using Eqn (2) from [1] 0072 end 0073 y(x==0)=0; % fix up special case of x=0 0074 y(x==Inf)=1; % fix up special case of x=Inf 0075 y=reshape(y,s); 0076 if ~nargout 0077 plot(x,y); 0078 xlabel('x'); 0079 ylabel(sprintf('I_{%d}(x)/I_{%d}(x)',v+1,v)); 0080 end