V_BESSELRATIOI calculate the inverse Bessel function ratio
Inputs:
r value of the Bessel function ratio besseli(v+1,s)./besseli(v,s) (scalar or matrix)
v denominator Bessel function order [0] [for now v=0 always]
p digits precision <=14 [5]
Outputs:
s value of s such that r=besseli(v+1,s)./besseli(v,s)
This function implements function VKAPPA from [1] for which a FORTRAN
implementation is available in [2].
Refs:
[1] G. W. Hill, Evaluation and Inversion of the Ratios of Modified Bessel
Functions, I 1 (x)/I 0 (x) and I 1.5 (x)/I 0.5 (x),
ACM Transactions on Mathematical Software (TOMS), Vol 7, pp199-208, 1981
[2] ACM Collected Algorithms (CALGO), http://calgo.acm.org/0001 function s=v_besselratioi(r,v,p) 0002 %V_BESSELRATIOI calculate the inverse Bessel function ratio 0003 % Inputs: 0004 % r value of the Bessel function ratio besseli(v+1,s)./besseli(v,s) (scalar or matrix) 0005 % v denominator Bessel function order [0] [for now v=0 always] 0006 % p digits precision <=14 [5] 0007 % 0008 % Outputs: 0009 % s value of s such that r=besseli(v+1,s)./besseli(v,s) 0010 % 0011 % This function implements function VKAPPA from [1] for which a FORTRAN 0012 % implementation is available in [2]. 0013 % 0014 % Refs: 0015 % [1] G. W. Hill, Evaluation and Inversion of the Ratios of Modified Bessel 0016 % Functions, I 1 (x)/I 0 (x) and I 1.5 (x)/I 0.5 (x), 0017 % ACM Transactions on Mathematical Software (TOMS), Vol 7, pp199-208, 1981 0018 % [2] ACM Collected Algorithms (CALGO), http://calgo.acm.org/ 0019 0020 % Copyright (C) Mike Brookes 2017-2018 0021 % Version: $Id: v_besselratioi.m 10865 2018-09-21 17:22:45Z dmb $ 0022 % 0023 % VOICEBOX is a MATLAB toolbox for speech processing. 0024 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0025 % 0026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0027 % This program is free software; you can redistribute it and/or modify 0028 % it under the terms of the GNU General Public License as published by 0029 % the Free Software Foundation; either version 2 of the License, or 0030 % (at your option) any later version. 0031 % 0032 % This program is distributed in the hope that it will be useful, 0033 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0034 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0035 % GNU General Public License for more details. 0036 % 0037 % You can obtain a copy of the GNU General Public License from 0038 % http://www.gnu.org/copyleft/gpl.html or by writing to 0039 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0040 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0041 if nargin<3 || isempty(p) 0042 p=5; 0043 else 0044 p=min(max(ceil(p),1),13); % digits precision required 0045 end 0046 if nargin<2 || isempty(v) 0047 v=0; % default to I1/I0 0048 end 0049 sr=size(r); 0050 r=r(:); % convert to a column vector 0051 nr=prod(sr); 0052 if v~=0 0053 error('v must be zero (for now)'); 0054 end 0055 s=repmat(-1,nr,1); % default is -1 for r<0 or r>=1 0056 y=2./(1-r); 0057 m1=r>=0 & r<1; % m1 is mask for valid range 0058 mn=m1; % mn is mask for using Newton iteration 0059 m=m1 & r<=0.85; % use inverse taylor series for 0<=r<=0.85 0060 if any(m) 0061 rm=r(m); 0062 mn(m)= rm>=0.642; % we will need a Newton iteration for 0.642<=r<=0.85 0063 xm=rm.^2; 0064 sm=(((rm-5.6076).*rm+5.0797).*rm-4.6494).*y(m).*xm-1; 0065 s(m)=((((sm.*xm+15.0).*xm+60.0).*xm/360.0+1.0).*xm-2.0).*rm./(xm-1); 0066 end 0067 m=m1 & ~m; % use continued fraction for 0.85<r<1 0068 if any(m) 0069 rm=r(m); 0070 mn(m)=rm<0.95; % we will need a Newton iteration for 0.85<r<0.95 0071 ym=y(m); 0072 mc=rm<0.95; 0073 if any(mc) 0074 rm(mc)=(-2326.0.*rm(mc)+4317.5526).*rm(mc) - 2001.035224; 0075 end 0076 if any(~mc) 0077 rm(~mc)=32.0./(120.0*rm(~mc)-131.5+ym(~mc)); 0078 end 0079 s(m)=(ym+1.0+3.0./(ym-5.0-12.0./(ym-10.0-rm)))*0.25; 0080 end 0081 if any(mn) % we need to do one or two Newton iterations 0082 rmn=r(mn); 0083 smn=s(mn); 0084 ymn=y(mn); 0085 ymn= ((0.00048*ymn-0.1589).*ymn+0.744).*ymn - 4.2932; % Gradient approximation ds/dr 0086 smn=smn+(v_besselratio(smn,0,p+1)-rmn).*ymn; % do a Newton iteration 0087 mr=(rmn>=0.75) & (rmn<=0.875); % need a second Newton iteration for 0.75<=r<=0.875 0088 if any(mr) 0089 smn(mr)=smn(mr)+(v_besselratio(smn(mr),0,p+1)-rmn(mr)).*ymn(mr); % do a second Newton iteration 0090 end 0091 s(mn)=smn; % update the main s() array 0092 end 0093 s=reshape(s,sr); % put back into the same shape is the input r 0094 if ~nargout 0095 plot(r,s); 0096 ylabel('x'); 0097 xlabel(sprintf('I_{%d}(x)/I_{%d}(x)',v+1,v)); 0098 end