Home > voicebox > besselratioi.m

besselratioi

PURPOSE ^

BESSELRATIOI calculate the inverse Bessel function ratio

SYNOPSIS ^

function s=besselratioi(r,v,p)

DESCRIPTION ^

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/

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function s=besselratioi(r,v,p)
0002 %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: besselratioi.m 10395 2018-02-22 17:00: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+(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)+(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

Generated on Wed 28-Mar-2018 15:33:24 by m2html © 2003