genhyper.m

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

sumr(1+2)=one;
numr(1+2)=one;
denomr(1+2)=one;
cnt=sigfig;
temp=complex(0.0,0.0);
oldtemp=temp;
ixcnt=0;
rexp=fix(ibit./2);
x=rexp.*(sumr(l+1+2)-2);
rr10=x.*log2;
ir10=fix(rr10);
rr10=rr10-ir10;
x=rexp.*(sumi(l+1+2)-2);
ri10=x.*log2;
ii10=fix(ri10);
ri10=ri10-ii10;
dum1=(abs(sumr(1+2).*rmax.*rmax+sumr(2+2).*rmax+sumr(3+2)).*sign(sumr(-1+2)));
dum2=(abs(sumi(1+2).*rmax.*rmax+sumi(2+2).*rmax+sumi(3+2)).*sign(sumi(-1+2)));
dum1=dum1.*10.^rr10;
dum2=dum2.*10.^ri10;
cdum1=complex(dum1,dum2);
x=rexp.*(denomr(l+1+2)-2);
rr10=x.*log2;
ir10=fix(rr10);
rr10=rr10-ir10;
x=rexp.*(denomi(l+1+2)-2);
ri10=x.*log2;
ii10=fix(ri10);
ri10=ri10-ii10;
dum1=(abs(denomr(1+2).*rmax.*rmax+denomr(2+2).*rmax+denomr(3+2)).*sign(denomr(-1+2)));
dum2=(abs(denomi(1+2).*rmax.*rmax+denomi(2+2).*rmax+denomi(3+2)).*sign(denomi(-1+2)));
dum1=dum1.*10.^rr10;
dum2=dum2.*10.^ri10;
cdum2=complex(dum1,dum2);
temp=cdum1./cdum2;
%
%     130 IF (IP .GT. 0) THEN
goon1=1;
while (goon1==1);
 goon1=0;
 if (ip<0) ;
  if (sumr(1+2)<half) ;
   mx1=sumi(l+1+2);
  elseif (sumi(1+2)<half);
   mx1=sumr(l+1+2);
  else;
   mx1=max([sumr(l+1+2),sumi(l+1+2)]);
  end;
  if (numr(1+2)<half) ;
   mx2=numi(l+1+2);
  elseif (numi(1+2)<half);
   mx2=numr(l+1+2);
  else;
   mx2=max([numr(l+1+2),numi(l+1+2)]);
  end;
  if (mx1-mx2>2.0) ;
   if (creal>=0.0) ;
    %        write (6,*) ' cdabs(temp1/cnt)=',cdabs(temp1/cnt)
    %
    if (abs(temp1./cnt)<=one) ;
     [sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit]=arydiv(sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit);
     hyper=final;
     return;
    end;
   end;
  end;
 else;
  [sumr,sumi,denomr,denomi,temp,l,lnpfq,rmax,ibit]=arydiv(sumr,sumi,denomr,denomi,temp,l,lnpfq,rmax,ibit);
  %
  %      First, estimate the exponent of the maximum term in the pFq
  %      series.
  %
  expon=0.0;
  xl=real(ixcnt);
  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);
  if (abs(oldtemp-temp)<abs(temp.*accy)) ;
   [sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit]=arydiv(sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit);
   hyper=final;
   return;
  end;
  oldtemp=temp;
 end;
 if (ixcnt~=icount) ;
  ixcnt=ixcnt+1;
  for i1=1 : iq;
   %
   %      TAKE THE CURRENT SUM AND MULTIPLY BY THE DENOMINATOR OF THE NEXT
   %
   %      TERM, FOR BOTH THE MOST SIGNIFICANT HALF (CR,CI) AND THE LEAST
   %
   %      SIGNIFICANT HALF (CR2,CI2).
   %
   %
   [sumr,sumi,cr(i1),ci(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(sumr,sumi,cr(i1),ci(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   [sumr,sumi,cr2(i1),ci2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(sumr,sumi,cr2(i1),ci2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   qr2(l+1+2)=qr2(l+1+2)-1;
   qi2(l+1+2)=qi2(l+1+2)-1;
   %
   %      STORE THIS TEMPORARILY IN THE SUM ARRAYS.
   %
   %
   [qr1,qi1,qr2,qi2,sumr,sumi,wk1,l,rmax]=cmpadd(qr1,qi1,qr2,qi2,sumr,sumi,wk1,l,rmax);
  end;
  %
  %
  %     MULTIPLY BY THE FACTORIAL TERM.
  %
  foo1=sumr;
  foo2=sumr;
  [foo1,cnt,foo2,wk6,l,rmax]=armult(foo1,cnt,foo2,wk6,l,rmax);
  sumr=foo2;
  foo1=sumi;
  foo2=sumi;
  [foo1,cnt,foo2,wk6,l,rmax]=armult(foo1,cnt,foo2,wk6,l,rmax);
  sumi=foo2;
  %
  %     MULTIPLY BY THE SCALING FACTOR, SIGFIG, TO KEEP THE SCALE CORRECT.
  %
  for i1=1 : ip-iq;
   foo1=sumr;
   foo2=sumr;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   sumr=foo2;
   foo1=sumi;
   foo2=sumi;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   sumi=foo2;
  end;
  for i1=1 : iq;
   %
   %      UPDATE THE DENOMINATOR.
   %
   %
   [denomr,denomi,cr(i1),ci(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(denomr,denomi,cr(i1),ci(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   [denomr,denomi,cr2(i1),ci2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(denomr,denomi,cr2(i1),ci2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   qr2(l+1+2)=qr2(l+1+2)-1;
   qi2(l+1+2)=qi2(l+1+2)-1;
   [qr1,qi1,qr2,qi2,denomr,denomi,wk1,l,rmax]=cmpadd(qr1,qi1,qr2,qi2,denomr,denomi,wk1,l,rmax);
  end;
  %
  %
  %     MULTIPLY BY THE FACTORIAL TERM.
  %
  foo1=denomr;
  foo2=denomr;
  [foo1,cnt,foo2,wk6,l,rmax]=armult(foo1,cnt,foo2,wk6,l,rmax);
  denomr=foo2;
  foo1=denomi;
  foo2=denomi;
  [foo1,cnt,foo2,wk6,l,rmax]=armult(foo1,cnt,foo2,wk6,l,rmax);
  denomi=foo2;
  %
  %     MULTIPLY BY THE SCALING FACTOR, SIGFIG, TO KEEP THE SCALE CORRECT.
  %
  for i1=1 : ip-iq;
   foo1=denomr;
   foo2=denomr;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   denomr=foo2;
   foo1=denomi;
   foo2=denomi;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   denomi=foo2;
  end;
  %
  %     FORM THE NEXT NUMERATOR TERM BY MULTIPLYING THE CURRENT
  %     NUMERATOR TERM (AN ARRAY) WITH THE A ARGUMENT (A SCALAR).
  %
  for i1=1 : ip;
   [numr,numi,ar(i1),ai(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(numr,numi,ar(i1),ai(i1),qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   [numr,numi,ar2(i1),ai2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(numr,numi,ar2(i1),ai2(i1),qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
   qr2(l+1+2)=qr2(l+1+2)-1;
   qi2(l+1+2)=qi2(l+1+2)-1;
   [qr1,qi1,qr2,qi2,numr,numi,wk1,l,rmax]=cmpadd(qr1,qi1,qr2,qi2,numr,numi,wk1,l,rmax);
  end;
  %
  %     FINISH THE NEW NUMERATOR TERM BY MULTIPLYING BY THE Z ARGUMENT.
  %
  [numr,numi,xr,xi,qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(numr,numi,xr,xi,qr1,qi1,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
  [numr,numi,xr2,xi2,qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax]=cmpmul(numr,numi,xr2,xi2,qr2,qi2,wk1,wk2,wk3,wk4,wk5,wk6,l,rmax);
  qr2(l+1+2)=qr2(l+1+2)-1;
  qi2(l+1+2)=qi2(l+1+2)-1;
  [qr1,qi1,qr2,qi2,numr,numi,wk1,l,rmax]=cmpadd(qr1,qi1,qr2,qi2,numr,numi,wk1,l,rmax);
  %
  %     MULTIPLY BY THE SCALING FACTOR, SIGFIG, TO KEEP THE SCALE CORRECT.
  %
  for i1=1 : iq-ip;
   foo1=numr;
   foo2=numr;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   numr=foo2;
   foo1=numi;
   foo2=numi;
   [foo1,sigfig,foo2,wk6,l,rmax]=armult(foo1,sigfig,foo2,wk6,l,rmax);
   numi=foo2;
  end;
  %
  %     FINALLY, ADD THE NEW NUMERATOR TERM WITH THE CURRENT RUNNING
  %     SUM OF THE NUMERATOR AND STORE THE NEW RUNNING SUM IN SUMR, SUMI.
  %
  foo1=sumr;
  foo2=sumr;
  bar1=sumi;
  bar2=sumi;
  [foo1,bar1,numr,numi,foo2,bar2,wk1,l,rmax]=cmpadd(foo1,bar1,numr,numi,foo2,bar2,wk1,l,rmax);
  sumi=bar2;
  sumr=foo2;

  %
  %     BECAUSE SIGFIG REPRESENTS "ONE" ON THE NEW SCALE, ADD SIGFIG
  %     TO THE CURRENT COUNT AND, CONSEQUENTLY, TO THE IP ARGUMENTS
  %     IN THE NUMERATOR AND THE IQ ARGUMENTS IN THE DENOMINATOR.
  %
  cnt=cnt+sigfig;
  for i1=1 : ip;
   ar(i1)=ar(i1)+sigfig;
  end;
  for i1=1 : iq;
   cr(i1)=cr(i1)+sigfig;
  end;
  goon1=1;
 end;
end;
[sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit]=arydiv(sumr,sumi,denomr,denomi,final,l,lnpfq,rmax,ibit);
%     write (6,*) 'Number of terms=',ixcnt
hyper=final;
return;


%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ARADD                             *
%     *                                                              *
%     *                                                              *
%     *  Description : Accepts two arrays of numbers and returns     *
%     *    the sum of the array.  Each array is holding the value    *
%     *    of one number in the series.  The parameter L is the      *
%     *    size of the array representing the number and RMAX is     *
%     *    the actual number of digits needed to give the numbers    *
%     *    the desired accuracy.                                     *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [a,b,c,z,l,rmax]=aradd(a,b,c,z,l,rmax);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
ediff=0;i=0;j=0;
%
%
for i=0 : l+1;
 z(i+2)=0.0;
end;
ediff=round(a(l+1+2)-b(l+1+2));
if (abs(a(1+2))<half | ediff<=-l) ;
 for i=-1 : l+1;
  c(i+2)=b(i+2);
 end;
 if (c(1+2)<half) ;
  c(-1+2)=one;
  c(l+1+2)=0.0;
 end;
 return;
else;
 if (abs(b(1+2))<half | ediff>=l) ;
  for i=-1 : l+1;
   c(i+2)=a(i+2);
  end;
  if (c(1+2)<half) ;
   c(-1+2)=one;
   c(l+1+2)=0.0;
  end;
  return;
 else;
  z(-1+2)=a(-1+2);
  goon300=1;
  goon190=1;
  if (abs(a(-1+2)-b(-1+2))>=half) ;
   goon300=0;
   if (ediff>0) ;
    z(l+1+2)=a(l+1+2);
   elseif (ediff<0);
    z(l+1+2)=b(l+1+2);
    z(-1+2)=b(-1+2);
    goon190=0;
   else;
    for i=1 : l;
     if (a(i+2)>b(i+2)) ;
      z(l+1+2)=a(l+1+2);
      break;
     end;
     if (a(i+2)<b(i+2)) ;
      z(l+1+2)=b(l+1+2);
      z(-1+2)=b(-1+2);
      goon190=0;
     end;
    end;
   end;
   %
  elseif (ediff>0);
   z(l+1+2)=a(l+1+2);
   for i=l : -1: 1+ediff ;
    z(i+2)=a(i+2)+b(i-ediff+2)+z(i+2);
    if (z(i+2)>=rmax) ;
     z(i+2)=z(i+2)-rmax;
     z(i-1+2)=one;
    end;
   end;
   for i=ediff : -1: 1 ;
    z(i+2)=a(i+2)+z(i+2);
    if (z(i+2)>=rmax) ;
     z(i+2)=z(i+2)-rmax;
     z(i-1+2)=one;
    end;
   end;
   if (z(0+2)>half) ;
    for i=l : -1: 1 ;
     z(i+2)=z(i-1+2);
    end;
    z(l+1+2)=z(l+1+2)+1;
    z(0+2)=0.0;
   end;
  elseif (ediff<0);
   z(l+1+2)=b(l+1+2);
   for i=l : -1: 1-ediff ;
    z(i+2)=a(i+ediff+2)+b(i+2)+z(i+2);
    if (z(i+2)>=rmax) ;
     z(i+2)=z(i+2)-rmax;
     z(i-1+2)=one;
    end;
   end;
   for i=0-ediff : -1: 1 ;
    z(i+2)=b(i+2)+z(i+2);
    if (z(i+2)>=rmax) ;
     z(i+2)=z(i+2)-rmax;
     z(i-1+2)=one;
    end;
   end;
   if (z(0+2)>half) ;
    for i=l : -1: 1 ;
     z(i+2)=z(i-1+2);
    end;
    z(l+1+2)=z(l+1+2)+one;
    z(0+2)=0.0;
   end;
  else;
   z(l+1+2)=a(l+1+2);
   for i=l : -1: 1 ;
    z(i+2)=a(i+2)+b(i+2)+z(i+2);
    if (z(i+2)>=rmax) ;
     z(i+2)=z(i+2)-rmax;
     z(i-1+2)=one;
    end;
   end;
   if (z(0+2)>half) ;
    for i=l : -1: 1 ;
     z(i+2)=z(i-1+2);
    end;
    z(l+1+2)=z(l+1+2)+one;
    z(0+2)=0.0;
   end;
  end;
  if (goon300==1) ;
   i=i; %here is the line that had a +1 taken from it.
   while (z(i+2)<half & i<l+1);
    i=i+1;
   end;
   if (i==l+1) ;
    z(-1+2)=one;
    z(l+1+2)=0.0;
    for i=-1 : l+1;
     c(i+2)=z(i+2);
    end;
    if (c(1+2)<half) ;
     c(-1+2)=one;
     c(l+1+2)=0.0;
    end;
    return;
   end;
   for j=1 : l+1-i;
    z(j+2)=z(j+i-1+2);
   end;
   for j=l+2-i : l;
    z(j+2)=0.0;
   end;
   z(l+1+2)=z(l+1+2)-i+1;
   for i=-1 : l+1;
    c(i+2)=z(i+2);
   end;
   if (c(1+2)<half) ;
    c(-1+2)=one;
    c(l+1+2)=0.0;
   end;
   return;
  end;
  %
  if (goon190==1) ;
   if (ediff>0) ;
    for i=l : -1: 1+ediff ;
     z(i+2)=a(i+2)-b(i-ediff+2)+z(i+2);
     if (z(i+2)<0.0) ;
      z(i+2)=z(i+2)+rmax;
      z(i-1+2)=-one;
     end;
    end;
    for i=ediff : -1: 1 ;
     z(i+2)=a(i+2)+z(i+2);
     if (z(i+2)<0.0) ;
      z(i+2)=z(i+2)+rmax;
      z(i-1+2)=-one;
     end;
    end;
   else;
    for i=l : -1: 1 ;
     z(i+2)=a(i+2)-b(i+2)+z(i+2);
     if (z(i+2)<0.0) ;
      z(i+2)=z(i+2)+rmax;
      z(i-1+2)=-one;
     end;
    end;
   end;
   if (z(1+2)>half) ;
    for i=-1 : l+1;
     c(i+2)=z(i+2);
    end;
    if (c(1+2)<half) ;
     c(-1+2)=one;
     c(l+1+2)=0.0;
    end;
    return;
   end;
   i=1;
   i=i+1;
   while (z(i+2)<half & i<l+1);
    i=i+1;
   end;
   if (i==l+1) ;
    z(-1+2)=one;
    z(l+1+2)=0.0;
    for i=-1 : l+1;
     c(i+2)=z(i+2);
    end;
    if (c(1+2)<half) ;
     c(-1+2)=one;
     c(l+1+2)=0.0;
    end;
    return;
   end;
   for j=1 : l+1-i;
    z(j+2)=z(j+i-1+2);
   end;
   for j=l+2-i : l;
    z(j+2)=0.0;
   end;
   z(l+1+2)=z(l+1+2)-i+1;
   for i=-1 : l+1;
    c(i+2)=z(i+2);
   end;
   if (c(1+2)<half) ;
    c(-1+2)=one;
    c(l+1+2)=0.0;
   end;
   return;
  end;
 end;
 %
 if (ediff<0) ;
  for i=l : -1: 1-ediff ;
   z(i+2)=b(i+2)-a(i+ediff+2)+z(i+2);
   if (z(i+2)<0.0) ;
    z(i+2)=z(i+2)+rmax;
    z(i-1+2)=-one;
   end;
  end;
  for i=0-ediff : -1: 1 ;
   z(i+2)=b(i+2)+z(i+2);
   if (z(i+2)<0.0) ;
    z(i+2)=z(i+2)+rmax;
    z(i-1+2)=-one;
   end;
  end;
 else;
  for i=l : -1: 1 ;
   z(i+2)=b(i+2)-a(i+2)+z(i+2);
   if (z(i+2)<0.0) ;
    z(i+2)=z(i+2)+rmax;
    z(i-1+2)=-one;
   end;
  end;
 end;
end;
%