# v_besselratio

## PURPOSE V_BESSELRATIO calculate the Bessel function ratio besseli(v+1,x)./besseli(v,x)

## SYNOPSIS function y=v_besselratio(x,v,p)

## DESCRIPTION ```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 
p digits precision <=14 

Outputs: y value of the Bessel function ratio besseli(v+1,x)./besseli(v,x)

Uses an algorithm from  of which a pseudocode version is given in  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.

    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.
    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.```

## CROSS-REFERENCE INFORMATION This function calls:
This function is called by:
• v_besselratioi V_BESSELRATIOI calculate the inverse Bessel function ratio

## SOURCE CODE ```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 
0006 %          p digits precision <=14 
0007 %
0008 % Outputs: y value of the Bessel function ratio besseli(v+1,x)./besseli(v,x)
0009 %
0010 % Uses an algorithm from  of which a pseudocode version is given in  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 %     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 %     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.
0026 %
0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0028 %   This program is free software; you can redistribute it and/or modify
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 
0059 end
0060 for i=1:n
0061     for k=1:n-i+1                   % k is actually k+1 in (20b) of 
0062         r(:,k)=x./(u+k+sqrt((u+k)^2+x.^2.*r(:,k+1)./r(:,k)));  % Eqn (20b) from 
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 
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```