genhyper.m

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

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;


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ARSUB                             *
%     *                                                              *
%     *                                                              *
%     *  Description : Accepts two arrays and subtracts each element *
%     *    in the second array from the element in the first array   *
%     *    and returns the solution.  The parameters L and RMAX are  *
%     *    the size of the array and the number of digits needed for *
%     *    the accuracy, respectively.                               *
%     *                                                              *
%     *  Subprograms called: ARADD                                   *
%     *                                                              *
%     ****************************************************************
%
function [a,b,c,wk1,wk2,l,rmax]=arsub(a,b,c,wk1,wk2,l,rmax);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
i=0;
%
%
for i=-1 : l+1;
 wk2(i+2)=b(i+2);
end;
wk2(-1+2)=(-one).*wk2(-1+2);
[a,wk2,c,wk1,l,rmax]=aradd(a,wk2,c,wk1,l,rmax);


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ARMULT                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Accepts two arrays and returns the product.   *
%     *    L and RMAX are the size of the arrays and the number of   *
%     *    digits needed to represent the numbers with the required  *
%     *    accuracy.                                                 *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [a,b,c,z,l,rmax]=armult(a,b,c,z,l,rmax);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
b2=0;carry=0;
i=0;
%
%
z(-1+2)=(abs(one).*sign(b)).*a(-1+2);
b2=abs(b);
z(l+1+2)=a(l+1+2);
for i=0 : l;
 z(i+2)=0.0;
end;
if (b2<=eps | a(1+2)<=eps) ;
 z(-1+2)=one;
 z(l+1+2)=0.0;
else;
 for i=l : -1: 1 ;
  z(i+2)=a(i+2).*b2+z(i+2);
  if (z(i+2)>=rmax) ;
   carry=fix(z(i+2)./rmax);
   z(i+2)=z(i+2)-carry.*rmax;
   z(i-1+2)=carry;
  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;
  if (z(1+2)>=rmax) ;
   for i=l : -1: 1 ;
    z(i+2)=z(i-1+2);
   end;
   carry=fix(z(1+2)./rmax);
   z(2+2)=z(2+2)-carry.*rmax;
   z(1+2)=carry;
   z(l+1+2)=z(l+1+2)+one;
  end;
  z(0+2)=0.0;
 end;
end;
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;


%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE CMPADD                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Takes two arrays representing one real and    *
%     *    one imaginary part, and adds two arrays representing      *
%     *    another complex number and returns two array holding the  *
%     *    complex sum.                                              *
%     *              (CR,CI) = (AR+BR, AI+BI)                        *
%     *                                                              *
%     *  Subprograms called: ARADD                                   *
%     *                                                              *
%     ****************************************************************
%
function [ar,ai,br,bi,cr,ci,wk1,l,rmax]=cmpadd(ar,ai,br,bi,cr,ci,wk1,l,rmax);
%
%
%
%
%
%
[ar,br,cr,wk1,l,rmax]=aradd(ar,br,cr,wk1,l,rmax);
[ai,bi,ci,wk1,l,rmax]=aradd(ai,bi,ci,wk1,l,rmax);


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE CMPSUB                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Takes two arrays representing one real and    *
%     *    one imaginary part, and subtracts two arrays representing *
%     *    another complex number and returns two array holding the  *
%     *    complex sum.                                              *
%     *              (CR,CI) = (AR+BR, AI+BI)                        *
%     *                                                              *
%     *  Subprograms called: ARADD                                   *
%     *                                                              *
%     ****************************************************************
%
function [ar,ai,br,bi,cr,ci,wk1,wk2,l,rmax]=cmpsub(ar,ai,br,bi,cr,ci,wk1,wk2,l,rmax);
%
%
%
%
%
%
[ar,br,cr,wk1,wk2,l,rmax]=arsub(ar,br,cr,wk1,wk2,l,rmax);
[ai,bi,ci,wk1,wk2,l,rmax]=arsub(ai,bi,ci,wk1,wk2,l,rmax);


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE CMPMUL                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Takes two arrays representing one real and    *
%     *    one imaginary part, and multiplies it with two arrays     *
%     *    representing another complex number and returns the       *
%     *    complex product.                                          *
%     *                                                              *
%     *  Subprograms called: ARMULT, ARSUB, ARADD                    *
%     *                                                              *
%     ****************************************************************
%
function [ar,ai,br,bi,cr,ci,wk1,wk2,cr2,d1,d2,wk6,l,rmax]=cmpmul(ar,ai,br,bi,cr,ci,wk1,wk2,cr2,d1,d2,wk6,l,rmax);
%
%
%
%
i=0;
%
%
[ar,br,d1,wk6,l,rmax]=armult(ar,br,d1,wk6,l,rmax);
[ai,bi,d2,wk6,l,rmax]=armult(ai,bi,d2,wk6,l,rmax);
[d1,d2,cr2,wk1,wk2,l,rmax]=arsub(d1,d2,cr2,wk1,wk2,l,rmax);
[ar,bi,d1,wk6,l,rmax]=armult(ar,bi,d1,wk6,l,rmax);
[ai,br,d2,wk6,l,rmax]=armult(ai,br,d2,wk6,l,rmax);
[d1,d2,ci,wk1,l,rmax]=aradd(d1,d2,ci,wk1,l,rmax);
for i=-1 : l+1;
 cr(i+2)=cr2(i+2);
end;


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ARYDIV                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Returns the double precision complex number   *
%     *    resulting from the division of four arrays, representing  *
%     *    two complex numbers.  The number returned will be in one  *
%     *    of two different forms:  either standard scientific or as *
%     *    the log (base 10) of the number.                          *
%     *                                                              *
%     *  Subprograms called: CONV21, CONV12, EADD, ECPDIV, EMULT.    *
%     *                                                              *
%     ****************************************************************
%
function [ar,ai,br,bi,c,l,lnpfq,rmax,ibit]=arydiv(ar,ai,br,bi,c,l,lnpfq,rmax,ibit);
%
%
cdum=[];ae=[];be=[];ce=[];n1=[];e1=[];n2=[];e2=[];n3=[];e3=[];
global  zero   half   one   two   ten   eps;
%
%
%
%
ae=zeros(2,2);be=zeros(2,2);ce=zeros(2,2);dum1=0;dum2=0;e1=0;e2=0;e3=0;n1=0;n2=0;n3=0;phi=0;ri10=0;rr10=0;tenmax=0;x=0;x1=0;x2=0;
cdum=0;
%
dnum=0;
ii10=0;ir10=0;itnmax=0;rexp=0;
%
%
%
rexp=fix(ibit./2);
x=rexp.*(ar(l+1+2)-2);
rr10=x.*log10(two)./log10(ten);
ir10=fix(rr10);
rr10=rr10-ir10;
x=rexp.*(ai(l+1+2)-2);
ri10=x.*log10(two)./log10(ten);
ii10=fix(ri10);
ri10=ri10-ii10;
dum1=(abs(ar(1+2).*rmax.*rmax+ar(2+2).*rmax+ar(3+2)).*sign(ar(-1+2)));
dum2=(abs(ai(1+2).*rmax.*rmax+ai(2+2).*rmax+ai(3+2)).*sign(ai(-1+2)));
dum1=dum1.*10.^rr10;
dum2=dum2.*10.^ri10;
cdum=complex(dum1,dum2);
[cdum,ae]=conv12(cdum,ae);
ae(1,2)=ae(1,2)+ir10;
ae(2,2)=ae(2,2)+ii10;
x=rexp.*(br(l+1+2)-2);
rr10=x.*log10(two)./log10(ten);
ir10=fix(rr10);
rr10=rr10-ir10;
x=rexp.*(bi(l+1+2)-2);
ri10=x.*log10(two)./log10(ten);
ii10=fix(ri10);
ri10=ri10-ii10;
dum1=(abs(br(1+2).*rmax.*rmax+br(2+2).*rmax+br(3+2)).*sign(br(-1+2)));
dum2=(abs(bi(1+2).*rmax.*rmax+bi(2+2).*rmax+bi(3+2)).*sign(bi(-1+2)));
dum1=dum1.*10.^rr10;
dum2=dum2.*10.^ri10;
cdum=complex(dum1,dum2);
[cdum,be]=conv12(cdum,be);
be(1,2)=be(1,2)+ir10;
be(2,2)=be(2,2)+ii10;
[ae,be,ce]=ecpdiv(ae,be,ce);
if (lnpfq==0) ;
 [ce,c]=conv21(ce,c);
else;
 [ce(1,1),ce(1,2),ce(1,1),ce(1,2),n1,e1]=emult(ce(1,1),ce(1,2),ce(1,1),ce(1,2),n1,e1);
 [ce(2,1),ce(2,2),ce(2,1),ce(2,2),n2,e2]=emult(ce(2,1),ce(2,2),ce(2,1),ce(2,2),n2,e2);
 [n1,e1,n2,e2,n3,e3]=eadd(n1,e1,n2,e2,n3,e3);
 n1=ce(1,1);
 e1=ce(1,2)-ce(2,2);
 x2=ce(2,1);
 %
 %      TENMAX - MAXIMUM SIZE OF EXPONENT OF 10
 %
 %      THE FOLLOWING CODE CAN BE USED TO DETERMINE TENMAX, BUT IT
 %
 %      WILL LIKELY GENERATE AN IEEE FLOATING POINT UNDERFLOW ERROR
 %
 %      ON A SUN WORKSTATION.  REPLACE TENMAX WITH THE VALUE APPROPRIATE
 %
 %      FOR YOUR MACHINE.
 %
 %
 tenmax=320;
 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-1;
 tenmax=real(itnmax);
 %
 if (e1>tenmax) ;
  x1=tenmax;
 elseif (e1<-tenmax);
  x1=0.0;
 else;
  x1=n1.*(ten.^e1);
 end;
 if (x2~=0.0) ;
  phi=atan2(x2,x1);
 else;
  phi=0.0;
 end;
 c=complex(half.*(log(n3)+e3.*log(ten)),phi);
end;


%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE EMULT                             *
%     *                                                              *
%     *                                                              *
%     *  Description : Takes one base and exponent and multiplies it *
%     *    by another numbers base and exponent to give the product  *
%     *    in the form of base and exponent.                         *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [n1,e1,n2,e2,nf,ef]=emult(n1,e1,n2,e2,nf,ef);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
nf=n1.*n2;
ef=e1+e2;
if (abs(nf)>=ten) ;
 nf=nf./ten;
 ef=ef+one;
end;


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE EDIV                              *
%     *                                                              *
%     *                                                              *
%     *  Description : returns the solution in the form of base and  *
%     *    exponent of the division of two exponential numbers.      *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [n1,e1,n2,e2,nf,ef]=ediv(n1,e1,n2,e2,nf,ef);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
nf=n1./n2;
ef=e1-e2;
if ((abs(nf)<one) & (nf~=zero)) ;
 nf=nf.*ten;
 ef=ef-one;
end;


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE EADD                              *
%     *                                                              *
%     *                                                              *
%     *  Description : Returns the sum of two numbers in the form    *
%     *    of a base and an exponent.                                *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [n1,e1,n2,e2,nf,ef]=eadd(n1,e1,n2,e2,nf,ef);
%
%
global  zero   half   one   two   ten   eps;
%
ediff=0;
%
%
ediff=e1-e2;
if (ediff>36.0d0) ;
 nf=n1;
 ef=e1;
elseif (ediff<-36.0d0);
 nf=n2;
 ef=e2;
else;
 nf=n1.*(ten.^ediff)+n2;
 ef=e2;
 while (1);
  if (abs(nf)<ten) ;
   while ((abs(nf)<one) & (nf~=0.0));
    nf=nf.*ten;
    ef=ef-one;
   end;
   break;
  else;
   nf=nf./ten;
   ef=ef+one;
  end;
 end;
end;


%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE ESUB                              *
%     *                                                              *
%     *                                                              *
%     *  Description : Returns the solution to the subtraction of    *
%     *    two numbers in the form of base and exponent.             *
%     *                                                              *
%     *  Subprograms called: EADD                                    *
%     *                                                              *
%     ****************************************************************
%
function [n1,e1,n2,e2,nf,ef]=esub(n1,e1,n2,e2,nf,ef);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
[n1,e1,dumvar3,e2,nf,ef]=eadd(n1,e1,n2.*(-one),e2,nf,ef);


%
%
%     ****************************************************************
%     *                                                              *
%     *                 SUBROUTINE CONV12                            *
%     *                                                              *
%     *                                                              *
%     *  Description : Converts a number from complex notation to a  *
%     *    form of a 2x2 real array.                                 *
%     *                                                              *
%     *  Subprograms called: none                                    *
%     *                                                              *
%     ****************************************************************
%
function [cn,cae]=conv12(cn,cae);
%
%
global  zero   half   one   two   ten   eps;
%
%
%
%