genhyper.m

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

%
cae(1,1)=real(cn);
cae(1,2)=0.0;
while (1);
 if (abs(cae(1,1))<ten) ;
  while (1);
   if ((abs(cae(1,1))>=one) | (cae(1,1)==0.0)) ;
    cae(2,1)=imag(cn);
    cae(2,2)=0.0;
    while (1);
     if (abs(cae(2,1))<ten) ;
      while ((abs(cae(2,1))<one) & (cae(2,1)~=0.0));
       cae(2,1)=cae(2,1).*ten;
       cae(2,2)=cae(2,2)-one;
      end;
      break;
     else;
      cae(2,1)=cae(2,1)./ten;
      cae(2,2)=cae(2,2)+one;
     end;
    end;
    break;
   else;
    cae(1,1)=cae(1,1).*ten;
    cae(1,2)=cae(1,2)-one;
   end;
  end;
  break;
 else;
  cae(1,1)=cae(1,1)./ten;
  cae(1,2)=cae(1,2)+one;
 end;
end;


%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE CONV21                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Converts a number represented in a 2x2 real   *
%     *    array to the form of a complex number.                    *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [cae,cn]=conv21(cae,cn);
%
%
%
global  zero   half   one   two   ten   eps;
global  nout;
%
%
%
dnum=0;tenmax=0;
itnmax=0;
%
%
%     TENMAX - MAXIMUM SIZE OF EXPONENT OF 10
%
itnmax=1;
dnum=0.1d0;
itnmax=itnmax+1;
dnum=dnum.*0.1d0;
while  (dnum>0.0);
 itnmax=itnmax+1;
 dnum=dnum.*0.1d0;
end;
itnmax=itnmax-2;
tenmax=real(itnmax);
%
if (cae(1,2)>tenmax | cae(2,2)>tenmax) ;
 %      CN=CMPLX(TENMAX,TENMAX)
 %
 itnmax,
 %format (' error - value of exponent required for summation',' was larger',./,' than the maximum machine exponent ',1i3,./,[' suggestions:'],./,' 1) re-run using lnpfq=1.',./,' 2) if you are using a vax, try using the',' fortran./g_floating option');
 error('stop encountered in original fortran code');
elseif (cae(2,2)<-tenmax);
 cn=complex(cae(1,1).*(10.^cae(1,2)),0.0);
else;
 cn=complex(cae(1,1).*(10.^cae(1,2)),cae(2,1).*(10.^cae(2,2)));
end;
return;


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ECPMUL                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Multiplies two numbers which are each         *
%     *    represented in the form of a two by two array and returns *
%     *    the solution in the same form.                            *
%     *                                                              *
%     *  Subprograms called: EMULT, ESUB, EADD                       *
%     *                                                              *
%     ****************************************************************
%
function [a,b,c]=ecpmul(a,b,c);
%
%
n1=[];e1=[];n2=[];e2=[];c2=[];
c2=zeros(2,2);e1=0;e2=0;n1=0;n2=0;
%
%
[a(1,1),a(1,2),b(1,1),b(1,2),n1,e1]=emult(a(1,1),a(1,2),b(1,1),b(1,2),n1,e1);
[a(2,1),a(2,2),b(2,1),b(2,2),n2,e2]=emult(a(2,1),a(2,2),b(2,1),b(2,2),n2,e2);
[n1,e1,n2,e2,c2(1,1),c2(1,2)]=esub(n1,e1,n2,e2,c2(1,1),c2(1,2));
[a(1,1),a(1,2),b(2,1),b(2,2),n1,e1]=emult(a(1,1),a(1,2),b(2,1),b(2,2),n1,e1);
[a(2,1),a(2,2),b(1,1),b(1,2),n2,e2]=emult(a(2,1),a(2,2),b(1,1),b(1,2),n2,e2);
[n1,e1,n2,e2,c(2,1),c(2,2)]=eadd(n1,e1,n2,e2,c(2,1),c(2,2));
c(1,1)=c2(1,1);
c(1,2)=c2(1,2);


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ECPDIV                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Divides two numbers and returns the solution. *
%     *    All numbers are represented by a 2x2 array.               *
%     *                                                              *
%     *  Subprograms called: EADD, ECPMUL, EDIV, EMULT               *
%     *                                                              *
%     ****************************************************************
%
function [a,b,c]=ecpdiv(a,b,c);
%
%
b2=[];c2=[];n1=[];e1=[];n2=[];e2=[];n3=[];e3=[];
global  zero   half   one   two   ten   eps;
%
b2=zeros(2,2);c2=zeros(2,2);e1=0;e2=0;e3=0;n1=0;n2=0;n3=0;
%
%
b2(1,1)=b(1,1);
b2(1,2)=b(1,2);
b2(2,1)=-one.*b(2,1);
b2(2,2)=b(2,2);
[a,b2,c2]=ecpmul(a,b2,c2);
[b(1,1),b(1,2),b(1,1),b(1,2),n1,e1]=emult(b(1,1),b(1,2),b(1,1),b(1,2),n1,e1);
[b(2,1),b(2,2),b(2,1),b(2,2),n2,e2]=emult(b(2,1),b(2,2),b(2,1),b(2,2),n2,e2);
[n1,e1,n2,e2,n3,e3]=eadd(n1,e1,n2,e2,n3,e3);
[c2(1,1),c2(1,2),n3,e3,c(1,1),c(1,2)]=ediv(c2(1,1),c2(1,2),n3,e3,c(1,1),c(1,2));
[c2(2,1),c2(2,2),n3,e3,c(2,1),c(2,2)]=ediv(c2(2,1),c2(2,2),n3,e3,c(2,1),c(2,2));


%     ****************************************************************
%     *                                                              *
%     *                   FUNCTION IPREMAX                           *
%     *                                                              *
%     *                                                              *
%     *  Description : Predicts the maximum term in the pFq series   *
%     *    via a simple scanning of arguments.                       *
%     *                                                              *
%     *  Subprograms called: none.                                   *
%     *                                                              *
%     ****************************************************************
%
function [ipremax]=ipremax(a,b,ip,iq,z);
%
%
%
global  zero   half   one   two   ten   eps;
global  nout;
%
%
%
%
%
expon=0;xl=0;xmax=0;xterm=0;
%
i=0;j=0;
%
xterm=0;
for j=1 : 100000;
 %
 %      Estimate the exponent of the maximum term in the pFq series.
 %
 %
 expon=zero;
 xl=real(j);
 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,zero)));
 xmax=log10(exp(one)).*expon;
 if ((xmax<xterm) & (j>2)) ;
  ipremax=j;
  return;
 end;
 xterm=max([xmax,xterm]);
end;
' error in ipremax--did not find maximum exponent',
error('stop encountered in original fortran code');


%     ****************************************************************
%     *                                                              *
%     *                   FUNCTION FACTOR                            *
%     *                                                              *
%     *                                                              *
%     *  Description : This function is the log of the factorial.    *
%     *                                                              *
%     *  Subprograms called: none.                                   *
%     *                                                              *
%     ****************************************************************
%
function [factor]=factor(z);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
pi=0;
%
if (((real(z)==one) & (imag(z)==zero)) | (abs(z)==zero)) ;
 factor=complex(zero,zero);
 return;
end;
pi=two.*two.*atan(one);
factor=half.*log(two.*pi)+(z+half).*log(z)-z+(one./(12.0d0.*z)).*(one-(one./(30.d0.*z.*z)).*(one-(two./(7.0d0.*z.*z))));


%     ****************************************************************
%     *                                                              *
%     *                   FUNCTION CGAMMA                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Calculates the complex gamma function.  Based *
%     *     on a program written by F.A. Parpia published in Computer*
%     *     Physics Communications as the `GRASP2' program (public   *
%     *     domain).                                                 *
%     *                                                              *
%     *                                                              *
%     *  Subprograms called: none.                                   *
%     *                                                              *
%     ****************************************************************
function [cgamma]=cgamma(arg,lnpfq);
%
%
%
%
%
global  zero   half   one   two   ten   eps;
global  nout;
%
%
%
argi=0;argr=0;argui=0;argui2=0;argum=0;argur=0;argur2=0;clngi=0;clngr=0;diff=0;dnum=0;expmax=0;fac=0;facneg=0;fd=zeros(7,1);fn=zeros(7,1);hlntpi=0;obasq=0;obasqi=0;obasqr=0;ovlfac=0;ovlfi=0;ovlfr=0;pi=0;precis=0;resi=0;resr=0;tenmax=0;tenth=0;termi=0;termr=0;twoi=0;zfaci=0;zfacr=0;
%
first=0;negarg=0;cgamma=0;
i=0;itnmax=0;
%
%
%
%----------------------------------------------------------------------*
%     *
%     THESE ARE THE BERNOULLI NUMBERS B02, B04, ..., B14, EXPRESSED AS *
%     RATIONAL NUMBERS. FROM ABRAMOWITZ AND STEGUN, P. 810.            *
%     *
fn=[1.0d00,-1.0d00,1.0d00,-1.0d00,5.0d00,-691.0d00,7.0d00];
fd=[6.0d00,30.0d00,42.0d00,30.0d00,66.0d00,2730.0d00,6.0d00];
%
%----------------------------------------------------------------------*
%
hlntpi=[1.0d00];
%
first=[true];
%
tenth=[0.1d00];
%
argr=real(arg);
argi=imag(arg);
%
%     ON THE FIRST ENTRY TO THIS ROUTINE, SET UP THE CONSTANTS REQUIRED
%     FOR THE REFLECTION FORMULA (CF. ABRAMOWITZ AND STEGUN 6.1.17) AND
%     STIRLING'S APPROXIMATION (CF. ABRAMOWITZ AND STEGUN 6.1.40).
%
if (first) ;
 pi=4.0d0.*atan(one);
 %
 %      SET THE MACHINE-DEPENDENT PARAMETERS:
 %
 %
 %      TENMAX - MAXIMUM SIZE OF EXPONENT OF 10
 %
 %
 itnmax=1;
 dnum=tenth;
 itnmax=itnmax+1;
 dnum=dnum.*tenth;
 while (dnum>0.0);
  itnmax=itnmax+1;
  dnum=dnum.*tenth;
 end;
 itnmax=itnmax-1;
 tenmax=real(itnmax);
 %
 %      EXPMAX - MAXIMUM SIZE OF EXPONENT OF E
 %
 %
 dnum=tenth.^itnmax;
 expmax=-log(dnum);
 %
 %      PRECIS - MACHINE PRECISION
 %
 %
 precis=one;
 precis=precis./two;
 dnum=precis+one;
 while (dnum>one);
  precis=precis./two;
  dnum=precis+one;
 end;
 precis=two.*precis;
 %
 hlntpi=half.*log(two.*pi);
 %
 for i=1 : 7;
  fn(i)=fn(i)./fd(i);
  twoi=two.*real(i);
  fn(i)=fn(i)./(twoi.*(twoi-one));
 end;
 %
 first=false;
 %
end;
%
%     CASES WHERE THE ARGUMENT IS REAL
%
if (argi==0.0) ;
 %
 %      CASES WHERE THE ARGUMENT IS REAL AND NEGATIVE
 %
 %
 if (argr<=0.0) ;
  %
  %       STOP WITH AN ERROR MESSAGE IF THE ARGUMENT IS TOO NEAR A POLE
  %
  %
  diff=abs(real(round(argr))-argr);
  if (diff<=two.*precis) ;
   ,
   argr , argi,
   %format (' argument (',1p,1d14.7,',',1d14.7,') too close to a',' pole.');
   error('stop encountered in original fortran code');
  else;
   %
   %        OTHERWISE USE THE REFLECTION FORMULA (ABRAMOWITZ AND STEGUN 6.1
   %        .17)
   %        TO ENSURE THAT THE ARGUMENT IS SUITABLE FOR STIRLING'S
   %
   %        FORMULA
   %
   %
   argum=pi./(-argr.*sin(pi.*argr));
   if (argum<0.0) ;
    argum=-argum;
    clngi=pi;
   else;
    clngi=0.0;
   end;
   facneg=log(argum);
   argur=-argr;
   negarg=true;
   %
  end;
  %
  %       CASES WHERE THE ARGUMENT IS REAL AND POSITIVE
  %
  %
 else;
  %
  clngi=0.0;
  argur=argr;
  negarg=false;
  %
 end;
 %
 %      USE ABRAMOWITZ AND STEGUN FORMULA 6.1.15 TO ENSURE THAT
 %
 %      THE ARGUMENT IN STIRLING'S FORMULA IS GREATER THAN 10
 %
 %
 ovlfac=one;
 while (argur<ten);
  ovlfac=ovlfac.*argur;
  argur=argur+one;
 end;
 %
 %      NOW USE STIRLING'S FORMULA TO COMPUTE LOG (GAMMA (ARGUM))
 %
 %
 clngr=(argur-half).*log(argur)-argur+hlntpi;
 fac=argur;
 obasq=one./(argur.*argur);
 for i=1 : 7;
  fac=fac.*obasq;
  clngr=clngr+fn(i).*fac;
 end;
 %
 %      INCLUDE THE CONTRIBUTIONS FROM THE RECURRENCE AND REFLECTION
 %
 %      FORMULAE
 %
 %
 clngr=clngr-log(ovlfac);
 if (negarg) ;
  clngr=facneg-clngr;
 end;
 %
else;
 %
 %      CASES WHERE THE ARGUMENT IS COMPLEX
 %
 %
 argur=argr;
 argui=argi;
 argui2=argui.*argui;
 %
 %      USE THE RECURRENCE FORMULA (ABRAMOWITZ AND STEGUN 6.1.15)
 %
 %      TO ENSURE THAT THE MAGNITUDE OF THE ARGUMENT IN STIRLING'S
 %
 %      FORMULA IS GREATER THAN 10
 %
 %
 ovlfr=one;
 ovlfi=0.0;
 argum=sqrt(argur.*argur+argui2);
 while (argum<ten);
  termr=ovlfr.*argur-ovlfi.*argui;
  termi=ovlfr.*argui+ovlfi.*argur;
  ovlfr=termr;
  ovlfi=termi;
  argur=argur+one;
  argum=sqrt(argur.*argur+argui2);
 end;
 %
 %      NOW USE STIRLING'S FORMULA TO COMPUTE LOG (GAMMA (ARGUM))
 %
 %
 argur2=argur.*argur;
 termr=half.*log(argur2+argui2);
 termi=atan2(argui,argur);
 clngr=(argur-half).*termr-argui.*termi-argur+hlntpi;
 clngi=(argur-half).*termi+argui.*termr-argui;
 fac=(argur2+argui2).^(-2);
 obasqr=(argur2-argui2).*fac;
 obasqi=-two.*argur.*argui.*fac;
 zfacr=argur;
 zfaci=argui;
 for i=1 : 7;
  termr=zfacr.*obasqr-zfaci.*obasqi;
  termi=zfacr.*obasqi+zfaci.*obasqr;
  fac=fn(i);
  clngr=clngr+termr.*fac;
  clngi=clngi+termi.*fac;
  zfacr=termr;
  zfaci=termi;
 end;
 %
 %      ADD IN THE RELEVANT PIECES FROM THE RECURRENCE FORMULA
 %
 %
 clngr=clngr-half.*log(ovlfr.*ovlfr+ovlfi.*ovlfi);
 clngi=clngi-atan2(ovlfi,ovlfr);
 %
end;
if (lnpfq==1) ;
 cgamma=complex(clngr,clngi);
 return;
end;
%
%     NOW EXPONENTIATE THE COMPLEX LOG GAMMA FUNCTION TO GET
%     THE COMPLEX GAMMA FUNCTION
%
if ((clngr<=expmax) & (clngr>=-expmax)) ;
 fac=exp(clngr);
else;
 ,
 clngr,
 %format (' argument to exponential function (',1p,1d14.7,') out of range.');
 error('stop encountered in original fortran code');
end;
resr=fac.*cos(clngi);
resi=fac.*sin(clngi);
cgamma=complex(resr,resi);
%
return;