genhyper.m

我认为很不错的语音处理的ｍａｔｌａｂ源代码

function [pfq]=genHyper(a,b,z,lnpfq,ix,nsigfig);
% function [pfq]=genHyper(a,b,z,lnpfq,ix,nsigfig)
% Description : A numerical evaluator for the generalized hypergeometric
%               function for complex arguments with large magnitudes
%               using a direct summation of the Gauss series.
%               pFq isdefined by (borrowed from Maple):
%   pFq = sum(z^k / k! * product(pochhammer(n[i], k), i=1..p) /
%         product(pochhammer(d[j], k), j=1..q), k=0..infinity )
%
% INPUTS:       a => array containing numerator parameters
%               b => array containing denominator parameters
%               z => complex argument (scalar)
%           lnpfq => (optional) set to 1 if desired result is the natural
%                    log of pfq (default is 0)
%              ix => (optional) maximum number of terms in a,b (see below)
%         nsigfig => number of desired significant figures (default=10)
%
% OUPUT:      pfq => result
%
% EXAMPLES:     a=[1+i,1]; b=[2-i,3,3]; z=1.5;
%               >> genHyper(a,b,z)
%               ans =
%                          1.02992154295955 +     0.106416425916656i
%               or with more precision,
%               >> genHyper(a,b,z,0,0,15)
%               ans =
%                          1.02992154295896 +     0.106416425915575i
%               using the log option,
%               >> genHyper(a,b,z,1,0,15)
%               ans =
%                        0.0347923403326305 +     0.102959427435454i
%               >> exp(ans)
%               ans =
%                          1.02992154295896 +     0.106416425915575i
%
%
% Translated from the original fortran using f2matlab.m
%  by Ben E. Barrowes - [e-mail], 7/04.
%  

%% Original fortran documentation
%     ACPAPFQ.  A NUMERICAL EVALUATOR FOR THE GENERALIZED HYPERGEOMETRIC
%
%     1  SERIES.  W.F. PERGER, A. BHALLA, M. NARDIN.
%
%     REF. IN COMP. PHYS. COMMUN. 77 (1993) 249
%
%     ****************************************************************
%     *                                                              *
%     *    SOLUTION TO THE GENERALIZED HYPERGEOMETRIC FUNCTION       *
%     *                                                              *
%     *                           by                                 *
%     *                                                              *
%     *                      W. F. PERGER,                           *
%     *                                                              *
%     *              MARK NARDIN  and ATUL BHALLA                    *
%     *                                                              *
%     *                                                              *
%     *            Electrical Engineering Department                 *
%     *            Michigan Technological University                 *
%     *                  1400 Townsend Drive                         *
%     *                Houghton, MI  49931-1295   USA                *
%     *                     Copyright 1993                           *
%     *                                                              *
%     *               e-mail address: [e-mail]                    *
%     *                                                              *
%     *  Description : A numerical evaluator for the generalized     *
%     *    hypergeometric function for complex arguments with large  *
%     *    magnitudes using a direct summation of the Gauss series.  *
%     *    The method used allows an accuracy of up to thirteen      *
%     *    decimal places through the use of large integer arrays    *
%     *    and a single final division.                              *
%     *    (original subroutines for the confluent hypergeometric    *
%     *    written by Mark Nardin, 1989; modifications made to cal-  *
%     *    culate the generalized hypergeometric function were       *
%     *    written by W.F. Perger and A. Bhalla, June, 1990)         *
%     *                                                              *
%     *  The evaluation of the pFq series is accomplished by a func- *
%     *  ion call to PFQ, which is a double precision complex func-  *
%     *  tion.  The required input is:                               *
%     *  1. Double precision complex arrays A and B.  These are the  *
%     *     arrays containing the parameters in the numerator and de-*
%     *     nominator, respectively.                                 *
%     *  2. Integers IP and IQ.  These integers indicate the number  *
%     *     of numerator and denominator terms, respectively (these  *
%     *     are p and q in the pFq function).                        *
%     *  3. Double precision complex argument Z.                     *
%     *  4. Integer LNPFQ.  This integer should be set to '1' if the *
%     *     result from PFQ is to be returned as the natural logaritm*
%     *     of the series, or '0' if not.  The user can generally set*
%     *     LNPFQ = '0' and change it if required.                   *
%     *  5. Integer IX.  This integer should be set to '0' if the    *
%     *     user desires the program PFQ to estimate the number of   *
%     *     array terms (in A and B) to be used, or an integer       *
%     *     greater than zero specifying the number of integer pos-  *
%     *     itions to be used.  This input parameter is escpecially  *
%     *     useful as a means to check the results of a given run.   *
%     *     Specificially, if the user obtains a result for a given  *
%     *     set of parameters, then changes IX and re-runs the eval- *
%     *     uator, and if the number of array positions was insuffi- *
%     *     cient, then the two results will likely differ.  The rec-*
%     *     commended would be to generally set IX = '0' and then set*
%     *     it to 100 or so for a second run.  Note that the LENGTH  *
%     *     parameter currently sets the upper limit on IX to 777,   *
%     *     but that can easily be changed (it is a single PARAMETER *
%     *     statement) and the program recompiled.                   *
%     *  6. Integer NSIGFIG.  This integer specifies the requested   *
%     *     number of significant figures in the final result.  If   *
%     *     the user attempts to request more than the number of bits*
%     *     in the mantissa allows, the program will abort with an   *
%     *     appropriate error message.  The recommended value is 10. *
%     *                                                              *
%     *     Note: The variable NOUT is the file to which error mess- *
%     *           ages are written (default is 6).  This can be      *
%     *           changed in the FUNCTION PFQ to accomodate re-      *
%     *           of output to another file                          *
%     *                                                              *
%     *  Subprograms called: HYPER.                                  *
%     *                                                              *
%     ****************************************************************
%
%
%
%

if nargin<6
 nsigfig=10;
elseif isempty(nsigfig)
 nsigfig=10;
end
if nargin<5
 ix=0;
elseif isempty(ix)
 ix=0;
end
if nargin<4
 lnpfq=0;
elseif isempty(lnpfq)
 lnpfq=0;
end
ip=length(a);
iq=length(b);

global  zero   half   one   two   ten   eps;
[zero , half , one , two , ten , eps]=deal(0.0d0,0.5d0,1.0d0,2.0d0,10.0d0,1.0d-10);
global  nout;
%
%
%
%
a1=zeros(2,1);b1=zeros(1,1);gam1=0;gam2=0;gam3=0;gam4=0;gam5=0;gam6=0;gam7=0;hyper1=0;hyper2=0;z1=0;
argi=0;argr=0;diff=0;dnum=0;precis=0;
%
%
i=0;
%
%
%
nout=6;
if ((lnpfq~=0) & (lnpfq~=1)) ;
 ' error in input arguments: lnpfq ~= 0 or 1',
 error('stop encountered in original fortran code');
end;
if ((ip>iq) & (abs(z)>one)) ;
 ip , iq , abs(z),
 %format [,1x,'ip=',1i2,3x,'iq=',1i2,3x,'and abs(z)=',1e12.5,2x,./,' which is greater than one--series does',' not converge');
 error('stop encountered in original fortran code');
end;
if (ip==2 & iq==1 & abs(z)>0.9) ;
 if (lnpfq~=1) ;
  %
  %      Check to see if the Gamma function arguments are o.k.; if not,
  %
  %      then the series will have to be used.
  %
  %
  %
  %      PRECIS - MACHINE PRECISION
  %
  %
  precis=one;
  precis=precis./two;
  dnum=precis+one;
  while (dnum>one);
   precis=precis./two;
   dnum=precis+one;
  end;
  precis=two.*precis;
  for i=1 : 6;
   if (i==1) ;
    argi=imag(b(1));
    argr=real(b(1));
   elseif (i==2);
    argi=imag(b(1)-a(1)-a(2));
    argr=real(b(1)-a(1)-a(2));
   elseif (i==3);
    argi=imag(b(1)-a(1));
    argr=real(b(1)-a(1));
   elseif (i==4);
    argi=imag(a(1)+a(2)-b(1));
    argr=real(a(1)+a(2)-b(1));
   elseif (i==5);
    argi=imag(a(1));
    argr=real(a(1));
   elseif (i==6);
    argi=imag(a(2));
    argr=real(a(2));
   end;
   %
   %       CASES WHERE THE ARGUMENT IS REAL
   %
   %
   if (argi==0.0) ;
    %
    %        CASES WHERE THE ARGUMENT IS REAL AND NEGATIVE
    %
    %
    if (argr<=0.0) ;
     %
     %         USE THE SERIES EXPANSION IF THE ARGUMENT IS TOO NEAR A POLE
     %
     %
     diff=abs(real(round(argr))-argr);
     if (diff<=two.*precis) ;
      pfq=hyper(a,b,ip,iq,z,lnpfq,ix,nsigfig);
      return;
     end;
    end;
   end;
  end;
  gam1=cgamma(b(1),lnpfq);
  gam2=cgamma(b(1)-a(1)-a(2),lnpfq);
  gam3=cgamma(b(1)-a(1),lnpfq);
  gam4=cgamma(b(1)-a(2),lnpfq);
  gam5=cgamma(a(1)+a(2)-b(1),lnpfq);
  gam6=cgamma(a(1),lnpfq);
  gam7=cgamma(a(2),lnpfq);
  a1(1)=a(1);
  a1(2)=a(2);
  b1(1)=a(1)+a(2)-b(1)+one;
  z1=one-z;
  hyper1=hyper(a1,b1,ip,iq,z1,lnpfq,ix,nsigfig);
  a1(1)=b(1)-a(1);
  a1(2)=b(1)-a(2);
  b1(1)=b(1)-a(1)-a(2)+one;
  hyper2=hyper(a1,b1,ip,iq,z1,lnpfq,ix,nsigfig);
  pfq=gam1.*gam2.*hyper1./(gam3.*gam4)+(one-z).^(b(1)-a(1)-a(2)).*gam1.*gam5.*hyper2./(gam6.*gam7);
  return;
 end;
end;
pfq=hyper(a,b,ip,iq,z,lnpfq,ix,nsigfig);
return;


%     ****************************************************************
%     *                                                              *
%     *                   FUNCTION BITS                              *
%     *                                                              *
%     *                                                              *
%     *  Description : Determines the number of significant figures  *
%     *    of machine precision to arrive at the size of the array   *
%     *    the numbers must be stored in to get the accuracy of the  *
%     *    solution.                                                 *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [bits]=bits;
%
bit2=0;
%
%
%
bit=1.0;
nnz=0;
nnz=nnz+1;
bit2=bit.*2.0;
bit=bit2+1.0;
while ((bit-bit2)~=0.0);
 nnz=nnz+1;
 bit2=bit.*2.0;
 bit=bit2+1.0;
end;
bits=nnz-3;


%     ****************************************************************
%     *                                                              *
%     *                   FUNCTION HYPER                             *
%     *                                                              *
%     *                                                              *
%     *  Description : Function that sums the Gauss series.          *
%     *                                                              *
%     *  Subprograms called: ARMULT, ARYDIV, BITS, CMPADD, CMPMUL,   *
%     *                      IPREMAX.                                *
%     *                                                              *
%     ****************************************************************
%
function [hyper]=hyper(a,b,ip,iq,z,lnpfq,ix,nsigfig);
%
%
% PARAMETER definitions
%
sumr=[];sumi=[];denomr=[];denomi=[];final=[];l=[];rmax=[];ibit=[];temp=[];cr=[];i1=[];ci=[];qr1=[];qi1=[];wk1=[];wk2=[];wk3=[];wk4=[];wk5=[];wk6=[];cr2=[];ci2=[];qr2=[];qi2=[];foo1=[];cnt=[];foo2=[];sigfig=[];numr=[];numi=[];ar=[];ai=[];ar2=[];ai2=[];xr=[];xi=[];xr2=[];xi2=[];bar1=[];bar2=[];
length=0;
length=777;
%
%
global  zero   half   one   two   ten   eps;
global  nout;
%
%
%
%
accy=0;ai=zeros(10,1);ai2=zeros(10,1);ar=zeros(10,1);ar2=zeros(10,1);ci=zeros(10,1);ci2=zeros(10,1);cnt=0;cr=zeros(10,1);cr2=zeros(10,1);creal=0;denomi=zeros(length+2,1);denomr=zeros(length+2,1);dum1=0;dum2=0;expon=0;log2=0;mx1=0;mx2=0;numi=zeros(length+2,1);numr=zeros(length+2,1);qi1=zeros(length+2,1);qi2=zeros(length+2,1);qr1=zeros(length+2,1);qr2=zeros(length+2,1);ri10=0;rmax=0;rr10=0;sigfig=0;sumi=zeros(length+2,1);sumr=zeros(length+2,1);wk1=zeros(length+2,1);wk2=zeros(length+2,1);wk3=zeros(length+2,1);wk4=zeros(length+2,1);wk5=zeros(length+2,1);wk6=zeros(length+2,1);x=0;xi=0;xi2=0;xl=0;xr=0;xr2=0;
%
cdum1=0;cdum2=0;final=0;oldtemp=0;temp=0;temp1=0;
%
%
%
i=0;i1=0;ibit=0;icount=0;ii10=0;ir10=0;ixcnt=0;l=0;lmax=0;nmach=0;rexp=0;
%
%
%
goon1=0;
foo1=zeros(length+2,1);foo2=zeros(length+2,1);bar1=zeros(length+2,1);bar2=zeros(length+2,1);
%
%
zero=0.0d0;
log2=log10(two);
ibit=fix(bits);
rmax=two.^(fix(ibit./2));
sigfig=two.^(fix(ibit./4));
%
for i1=1 : ip;
 ar2(i1)=real(a(i1)).*sigfig;
 ar(i1)=fix(ar2(i1));
 ar2(i1)=round((ar2(i1)-ar(i1)).*rmax);
 ai2(i1)=imag(a(i1)).*sigfig;
 ai(i1)=fix(ai2(i1));
 ai2(i1)=round((ai2(i1)-ai(i1)).*rmax);
end;
for i1=1 : iq;
 cr2(i1)=real(b(i1)).*sigfig;
 cr(i1)=fix(cr2(i1));
 cr2(i1)=round((cr2(i1)-cr(i1)).*rmax);
 ci2(i1)=imag(b(i1)).*sigfig;
 ci(i1)=fix(ci2(i1));
 ci2(i1)=round((ci2(i1)-ci(i1)).*rmax);
end;
xr2=real(z).*sigfig;
xr=fix(xr2);
xr2=round((xr2-xr).*rmax);
xi2=imag(z).*sigfig;
xi=fix(xi2);
xi2=round((xi2-xi).*rmax);
%
%     WARN THE USER THAT THE INPUT VALUE WAS SO CLOSE TO ZERO THAT IT
%     WAS SET EQUAL TO ZERO.
%
for i1=1 : ip;
 if ((real(a(i1))~=0.0) & (ar(i1)==0.0) & (ar2(i1)==0.0));
  i1,
 end;
 %format (1x,'warning - real part of a(',1i2,') was set to zero');
 if ((imag(a(i1))~=0.0) & (ai(i1)==0.0) & (ai2(i1)==0.0));
  i1,
 end;
 %format (1x,'warning - imag part of a(',1i2,') was set to zero');
end;
for i1=1 : iq;
 if ((real(b(i1))~=0.0) & (cr(i1)==0.0) & (cr2(i1)==0.0));
  i1,
 end;
 %format (1x,'warning - real part of b(',1i2,') was set to zero');
 if ((imag(b(i1))~=0.0) & (ci(i1)==0.0) & (ci2(i1)==0.0));
  i1,
 end;
 %format (1x,'warning - imag part of b(',1i2,') was set to zero');
end;
if ((real(z)~=0.0) & (xr==0.0) & (xr2==0.0)) ;
 ' warning - real part of z was set to zero',
 z=complex(0.0,imag(z));
end;
if ((imag(z)~=0.0) & (xi==0.0) & (xi2==0.0)) ;
 ' warning - imag part of z was set to zero',
 z=complex(real(z),0.0);
end;
%
%
%     SCREENING OF NUMERATOR ARGUMENTS FOR NEGATIVE INTEGERS OR ZERO.
%     ICOUNT WILL FORCE THE SERIES TO TERMINATE CORRECTLY.
%
nmach=fix(log10(two.^fix(bits)));
icount=-1;
for i1=1 : ip;
 if ((ar2(i1)==0.0) & (ar(i1)==0.0) & (ai2(i1)==0.0) &(ai(i1)==0.0)) ;
  hyper=complex(one,0.0);
  return;
 end;
 if ((ai(i1)==0.0) & (ai2(i1)==0.0) & (real(a(i1))<0.0));
  if (abs(real(a(i1))-real(round(real(a(i1)))))<ten.^(-nmach)) ;
   if (icount~=-1) ;
    icount=min([icount,-round(real(a(i1)))]);
   else;
    icount=-round(real(a(i1)));
   end;
  end;
 end;
end;
%
%     SCREENING OF DENOMINATOR ARGUMENTS FOR ZEROES OR NEGATIVE INTEGERS
%     .
%
for i1=1 : iq;
 if ((cr(i1)==0.0) & (cr2(i1)==0.0) & (ci(i1)==0.0) &(ci2(i1)==0.0)) ;
  i1,
  %format (1x,'error - argument b(',1i2,') was equal to zero');
  error('stop encountered in original fortran code');
 end;
 if ((ci(i1)==0.0) & (ci2(i1)==0.0) & (real(b(i1))<0.0));
  if ((abs(real(b(i1))-real(round(real(b(i1)))))<ten.^(-nmach)) &(icount>=-round(real(b(i1))) | icount==-1)) ;
   i1,
   %format (1x,'error - argument b(',1i2,') was a negative',' integer');
   error('stop encountered in original fortran code');
  end;
 end;
end;
%
nmach=fix(log10(two.^ibit));
nsigfig=min([nsigfig,fix(log10(two.^ibit))]);
accy=ten.^(-nsigfig);
l=ipremax(a,b,ip,iq,z);
if (l~=1) ;
 %
 %     First, estimate the exponent of the maximum term in the pFq series
 %     .
 %
 expon=0.0;
 xl=real(l);
 for i=1 : ip;
  expon=expon+real(factor(a(i)+xl-one))-real(factor(a(i)-one));
 end;
 for i=1 : iq;
  expon=expon-real(factor(b(i)+xl-one))+real(factor(b(i)-one));
 end;
 expon=expon+xl.*real(log(z))-real(factor(complex(xl,0.0)));
 lmax=fix(log10(exp(one)).*expon);
 l=lmax;
 %
 %     Now, estimate the exponent of where the pFq series will terminate.
 %
 temp1=complex(one,0.0);
 creal=one;
 for i1=1 : ip;
  temp1=temp1.*complex(ar(i1),ai(i1))./sigfig;
 end;
 for i1=1 : iq;
  temp1=temp1./(complex(cr(i1),ci(i1))./sigfig);
  creal=creal.*cr(i1);
 end;
 temp1=temp1.*complex(xr,xi);
 %
 %     Triple it to make sure.
 %
 l=3.*l;
 %
 %     Divide the number of significant figures necessary by the number
 %     of
 %     digits available per array position.
 %
 %
 l=fix((2.*l+nsigfig)./nmach)+2;
end;
%
%     Make sure there are at least 5 array positions used.
%
l=max([l,5]);
l=max([l,ix]);
%     write (6,*) ' Estimated value of L=',L
if ((l<0) | (l>length)) ;
 length,
 %format (1x,['error in fn hyper: l must be < '],1i4);
 error('stop encountered in original fortran code');
end;
if (nsigfig>nmach) ;
 nmach,
 %format (1x,' warning--the number of significant figures requ','ested',./,'is greater than the machine precision--','final answer',./,'will be accurate to only',i3,' digits');
end;
%
sumr(-1+2)=one;
sumi(-1+2)=one;
numr(-1+2)=one;
numi(-1+2)=one;
denomr(-1+2)=one;
denomi(-1+2)=one;
for i=0 : l+1;
 sumr(i+2)=0.0;
 sumi(i+2)=0.0;
 numr(i+2)=0.0;
 numi(i+2)=0.0;
 denomr(i+2)=0.0;
 denomi(i+2)=0.0;
end;