Home > voicebox > besselratio.m

besselratio

PURPOSE ^

BESSELRATIO calculate the Bessel function ratio besseli(v+1,x)./besseli(v,x)

SYNOPSIS ^

function y=besselratio(x,v,p)

DESCRIPTION ^

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): 239251, 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:

SOURCE CODE ^

0001 function y=besselratio(x,v,p)
0002 %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): 239251, 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 1997-2011
0022 %      Version: $Id: besselratio.m 10207 2017-10-09 10:17:27Z dmb $
0023 %
0024 %   VOICEBOX is a MATLAB toolbox for speech processing.
0025 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0026 %
0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0028 %   This program is free software; you can redistribute it and/or modify
0029 %   it under the terms of the GNU General Public License as published by
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);

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003