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

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

## 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 [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
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 % [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 %      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 [1]
0059 end
0060 for i=1:n
0061     for k=1:n-i+1                   % k is actually k+1 in (20b) of [1]
0062         r(:,k)=x./(u+k+sqrt((u+k)^2+x.^2.*r(:,k+1)./r(:,k)));  % Eqn (20b) from [1]
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 [1]
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```