# 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 [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.

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- v_besratinv0 V_BESRATINV0 Inverse function of the Modified Bessel Ratio I1(k)/I0(k)
- 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 [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 Lesser General Public License as published by
    0034 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0041 %
    0042 %   You can obtain a copy of the GNU Lesser General Public License from
    0043 %   https://www.gnu.org/licenses/ .
    0044 %
    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