pascal（delphi7）写的一个thht类，能实现emd分解、hilbert变换算法.txt

Pascal（Delphi7）写的一个THHT类

  for i:=low(indmax) to high(indmax) do
  begin
    tmax[i+2]:=indmax[i];
    mmax[i+2]:=m[indmax[i]];
  end;
  //能不能换成Copy


  BoundaryExpand(m,tmin,tmax,mmin,mmax);       //*

  setlength(y2Min,nMin);
  setlength(y2Max,nMax);
  p1Min:=(mmin[1]-mmin[0])/(tmin[1]-tmin[0]+0.0000001);
  pnMin:=(mmin[nMin-1]-mmin[nMin-2])/(tmin[nMin-1]-tmin[nMin-2]+0.0000001);
  p1Max:=(mmax[1]-mmax[0])/(tmax[1]-tmax[0]+0.0000001);
  pnMax:=(mmax[nMax-1]-mmax[nMax-2])/(tmax[nMax-1]-tmax[nMax-2]+0.0000001);

  spline(tmin,mmin,nMin-1,p1Min,pnMin,y2Min);
  spline(tmax,mmax,nMax-1,p1Max,pnMax,y2Max);

  setlength(tempEnvMin,(tmin[nMin-1]-tmin[0]+1));
  setlength(tempEnvMax,(tmax[nMax-1]-tmax[0]+1));
  ind1:=abs(tmin[0]);
  ind2:=abs(tmax[0]);
  for i:=0 to (tmin[nMin-1]-tmin[0]) do
    splint(tmin,mmin,y2Min,nMin-1,i-ind1,tempEnvMin[i]);
  for i:=0 to (tmax[nMax-1]-tmax[0]) do
    splint(tmax,mmax,y2Max,nMax-1,i-ind2,tempEnvMax[i]);

  n:=length(m);
  setlength(envmoy,n);
  setlength(amp,n);
  for i:=0 to n-1 do
  begin
    envmoy[i]:=(tempEnvMin[i+ind1]+tempEnvMax[i+ind2])/2;
    amp[i]:=2/abs(tempEnvMax[i+ind2]-tempEnvMin[i+ind1]+0.0000001);         //小心除数为零
  end;

end;


////////////////////////////////////////////////////////////////////////////////
procedure TEmd.extr(const datas:myArray;var indmin,indmax,indzero:myIntArray);
var
  m:integer;     //m(信号长度)
  i:integer;
  nMin,nMax,nZero:integer;        //极值点个数
  diffData:myArray;
  bad:myIntArray;           //导数为零的点
  nBad,nDebs,nFins:integer;
  dd,debs,fins:myIntArray;
  iz,diffIz,zer,dz,debz,finz:myIntArray;
  nIz,nz,nDebz,nFinz:integer;
begin
  m:=length(datas);
  setlength(diffData,m-1);
  nMin:=0;
  nMax:=0;
  for i:=0 to m-2 do diffData[i]:=datas[i+1]-datas[i];
  for i:=0 to m-3 do
  begin
    if (diffData[i]*diffData[i+1]<0) and (diffData[i]<0) then
    begin
      inc(nMin);
      setlength(indmin,nMin);
      indmin[nMin-1]:=i+1;
    end;
    if (diffData[i]*diffData[i+1]<0) and (diffData[i]>0) then
    begin
      inc(nMax);
      setlength(indmax,nMax);
      indmax[nMax-1]:=i+1;
    end;
  end;
  {如果有多个值同为极值点,求得bad}
  nBad:=0;
  setlength(bad,m-1);
  for i:=0 to m-2 do
    if diffData[i]=0 then
    begin
      inc(nBad);
      bad[i]:=1;
    end
    else
      bad[i]:=0;

  if nBad>0 then                  //下面的代码必须包含在这个if结构中
  begin
    setlength(dd,m);
    dd[0]:=bad[0]-0;
    for i:=0 to m-3 do dd[i+1]:=bad[i+1]-bad[i];
    dd[m-1]:=0-bad[m-2];                             //dd

    nDebs:=0;
    nFins:=0;
    for i:=0 to m-1 do
    begin
      if dd[i]=1 then       //debs
      begin
        inc(nDebs);
        setlength(debs,nDebs);
        debs[nDebs-1]:=i;
      end;
      if dd[i]=-1 then      //fins
      begin
        inc(nFins);
        setlength(fins,nFins);
        fins[nFins-1]:=i;
      end;
    end;

    if debs[0]=1 then        //
    begin
      if length(debs)>1 then
      begin
        for i:=0 to nDebs-2 do debs[i]:=debs[i+1];
        dec(nDebs);
        setlength(debs,nDebs);
        for i:=0 to nFins-2 do fins[i]:=fins[i+1];
        dec(nFins);
        setlength(fins,nFins);
      end
      else
      begin
        nDebs:=0;
        nFins:=0;
        setlength(debs,nDebs);
        setlength(fins,nFins);
      end;
    end;

    if length(debs)>0 then
    begin
      if fins[nFins-1]=m then
      begin
        if length(debs)>1 then
        begin
          dec(nDebs);
          dec(nFins);
          setlength(debs,nDebs);
          setlength(fins,nFins);
        end
        else
        begin
          nDebs:=0;
          nFins:=0;
          setlength(debs,nDebs);
          setlength(fins,nFins);
        end;
      end;
    end;

    if nDebs>0 then
    begin
      for i:=0 to nDebs-1 do
      begin
        if diffData[debs[i]-1]>0 then
        begin
          if diffData[fins[i]]<0 then
          begin
            inc(nMax);
            setlength(indmax,nMax);
            indmax[nMax-1]:=round((fins[i]+debs[i])/2);
          end;
        end
        else
        begin
          if diffData[fins[i]]>0 then
          begin
            inc(nMin);
            setlength(indmin,nMin);
            indmin[nMin-1]:=round((fins[i]+debs[i])/2);
          end;
        end;
      end;      //for i:=0 to nDebs-1
    end;    //if nDebs>0

    sort(indmax);
    sort(indmin);

  end;   //if nBad>0

  {下面是indzero}
  nZero:=0;
  for i:=0 to m-2 do
    if datas[i]*datas[i+1]<0 then
    begin
      inc(nZero);
      setlength(indzero,nZero);
      indzero[nZero-1]:=i;
    end;
  nIz:=0;
  for i:=0 to m-1 do
    if datas[i]=0 then
    begin
      inc(nIz);
      setlength(iz,nIz);
      iz[nIz-1]:=i;
    end;

  if nIz>0 then
  begin
    setlength(diffIz,nIz-1);
    nz:=0;
    for i:=0 to nIz-2 do
    begin
      diffIz[i]:=iz[i+1]-iz[i];
      if diffIz[i]=0 then inc(nz);
    end;
    if nz>0 then
    begin
      setlength(zer,m);
      for i:=0 to m-1 do
        if datas[i]=0 then zer[i]:=1 else zer[i]:=0;      //zer

      setlength(dz,m+1);
      dz[0]:=zer[0]-0;
      for i:=0 to m-2 do dz[i+1]:=zer[i+1]-zer[i];
      dz[m]:=0-zer[m-1];                 //dz

      nDebz:=0;
      nFinz:=0;
      for i:=0 to m do
      begin
        if dz[i]=1 then
        begin
          inc(nDebz);
          setlength(debz,nDebz);
          debz[nDebz-1]:=i;
        end;
        if dz[i]=-1 then
        begin
          inc(nFinz);
          setlength(finz,nFinz);
          finz[nFinz-1]:=i-1;
        end;
      end; // i:=0 to m

      for i:=0 to nDebz-1 do
      begin
        inc(nZero);
        setlength(indzero,nZero);
        indzero[nZero-1]:=round((debz[i]+finz[i])/2);
      end;
    end
    else
    begin
      for i:=0 to nIz-1 do
      begin
        inc(nZero);
        setlength(indzero,nZero);
        indzero[nZero-1]:=iz[i];
      end;
    end;    //if nz>0
    sort(indzero);
  end;      //if nIz>0

end;

///////////////////////////////////////////////////////////////////////////////////////////
procedure TEmd.Splint(XA:myIntArray; YA, Y2A:myarray; N:integer; X:integer;var Y:double);
label 1;
var
    K,KLO,KHI:integer;
    H,A,B,AAA,BBB,Q,QQ:double;
begin
    KLO:=1;
    KHI:=N;
1:  If KHI - KLO > 1 Then
    begin
        K:=(KHI + KLO) Div 2;
        If XA[K] > X Then
            KHI:=K
        Else
            KLO:=K;
        GoTo 1;
    end;
    H:=XA[KHI] - XA[KLO];
    If H = 0 Then
    begin
        ShowMessage('  PAUSE  "BAD  XA  INPUT"');
        Exit;
    end;
    A:=(XA[KHI] - X) / H;
    B:=(X - XA[KLO]) / H;
    AAA:=A * YA[KLO] + B * YA[KHI];
    if A = 0 then
      Q:= 0
    else
      if A > 0 then
        Q:= exp(3*ln(A))
      else
        Q:= -EXP(3*LN(-A));
    if B = 0 then
      QQ:= 0
    else
      if B > 0 then
        QQ:= exp(3*ln(B))
      else
        QQ:= -EXP(3*LN(-B));
    BBB:=(Q - A) * Y2A[KLO] + (QQ - B) * Y2A[KHI];
    Y:=AAA + BBB * (H * H) / 6;
end;


///////////////////////////////////////////////////////////////////////////////////////
procedure TEmd.Spline(X:myIntArray;Y:myArray;N:integer;YP1,YPN:double;var Y2:myArray);
var
    U:array of double;
    AAA,SIG,BBB,CCC,P,QN,UN:double;
    I,K:integer;
begin
    setlength(U,N);
    If YP1 > 9.9E+29 Then
    begin
        Y2[1]:=0;
        U[1]:=0;
    end
    Else
    begin
        Y2[1]:=-0.5;
        AAA:=(Y[2] - Y[1]) / (X[2] - X[1]);
        U[1]:=(3 / (X[2] - X[1])) * (AAA - YP1);
    end;
    For I:=2 To N - 1 do
    begin
        SIG:=(X[I] - X[I - 1]) / (X[I + 1] - X[I - 1]);
        P:=SIG * Y2[I - 1] + 2;
        Y2[I]:=(SIG - 1) / P;

        AAA:=(Y[I + 1] - Y[I]) / (X[I + 1] - X[I]);
        BBB:=(Y[I] - Y[I - 1]) / (X[I] - X[I - 1]);    //改过了
        CCC:=X[I + 1] - X[I - 1];
        U[I]:=(6 * (AAA - BBB) / CCC - SIG * U[I - 1]) / P;
    end;
    If YPN > 9.9E+29 Then
    begin
        QN:=0;
        UN:=0;
    end
    Else
    begin
        QN:=0.5;
        AAA:=YPN - (Y[N] - Y[N - 1]) / (X[N] - X[N - 1]);
        UN:=(3 / (X[N] - X[N - 1])) * AAA;
    end;
    Y2[N]:=(UN - QN * U[N - 1]) / (QN * Y2[N - 1] + 1);
    For K:=N-1 DownTo 1 do
        Y2[K]:=Y2[K]*Y2[K + 1] + U[K];
end;


///////////////////////////////////////////////////////////////////////////////////////////////////////
procedure TEmd.BoundaryExpand(const datas:myArray;var mins,maxes:myIntArray;var mins_h,maxes_h:myArray);
var
  h_maxes,h_mins,h_datas:integer;
  num:integer;
begin
  num:=high(datas);
  if maxes[2]<mins[2] then
  begin
     if datas[0]>mins_h[2] then
     begin
       maxes[0]:=2*maxes[2]-maxes[4];
       maxes[1]:=2*maxes[2]-maxes[3];
       maxes_h[0]:=maxes_h[4];
       maxes_h[1]:=maxes_h[3];
       mins[0]:=2*maxes[2]-mins[3];
       mins[1]:=2*maxes[2]-mins[2];
       mins_h[0]:=mins_h[3];
       mins_h[1]:=mins_h[2];
     end
     else
     begin
       maxes[1]:=0-maxes[2];
       maxes[0]:=0-maxes[3];
       maxes_h[1]:=maxes_h[2];
       maxes_h[0]:=maxes_h[3];
       mins[1]:=0;
       mins[0]:=0-mins[2];
       mins_h[1]:=datas[0];
       mins_h[0]:=mins_h[2];
     end;
  end
  else
  begin
    if datas[0]<maxes_h[2] then
    begin
      maxes[1]:=2*mins[2]-maxes[2];
      maxes[0]:=2*mins[2]-maxes[3];
      maxes_h[1]:=maxes_h[2];
      maxes_h[0]:=maxes_h[3];
      mins[1]:=2*mins[2]-mins[3];
      mins[0]:=2*mins[2]-mins[4];
      mins_h[1]:=mins_h[3];
      mins_h[0]:=mins_h[4];
    end
    else
    begin                           //否则以第一个点作为对称中心
      maxes[1]:=0;
      maxes[0]:=0-maxes[2];
      maxes_h[1]:=datas[0];
      maxes_h[0]:=maxes_h[2];
      mins[1]:=0-mins[2];
      mins[0]:=0-mins[3];
      mins_h[1]:=mins_h[2];
      mins_h[0]:=mins_h[3];
    end;
  end;

  h_maxes:=high(maxes)-2;
  h_mins:=high(mins)-2;
  h_datas:=high(datas);
  //右延拓
  if maxes[h_maxes]<mins[h_mins] then
  begin
    if datas[h_datas]<maxes_h[h_maxes] then   //如果端点的数值比最后一个极大值小
    begin
      maxes[h_maxes+1]:=2*mins[h_mins]-maxes[h_maxes];
      maxes[h_maxes+2]:=2*mins[h_mins]-maxes[h_maxes-1];
      maxes_h[h_maxes+1]:=maxes_h[h_maxes];
      maxes_h[h_maxes+2]:=maxes_h[h_maxes-1];
      mins[h_mins+1]:=2*mins[h_mins]-mins[h_mins-1];
      mins[h_mins+2]:=2*mins[h_mins]-mins[h_mins-2];
      mins_h[h_mins+1]:=mins_h[h_mins-1];
      mins_h[h_mins+2]:=mins_h[h_mins-2];
    end
    else
    begin
      maxes[h_maxes+1]:=num;
      maxes[h_maxes+2]:=2*num-maxes[h_maxes];
      maxes_h[h_maxes+1]:=datas[h_datas];
      maxes_h[h_maxes+2]:=maxes_h[h_maxes];
      mins[h_mins+1]:=2*num-mins[h_mins];
      mins[h_mins+2]:=2*num-mins[h_mins-1];
      mins_h[h_mins+1]:=mins_h[h_mins];
      mins_h[h_mins+2]:=mins_h[h_mins-1];
    end;
  end
  else
  begin
    if datas[h_datas]>mins_h[h_mins] then
    begin
      maxes[h_maxes+1]:=2*maxes[h_maxes]-maxes[h_maxes-1];
      maxes[h_maxes+2]:=2*maxes[h_maxes]-maxes[h_maxes-2];
      maxes_h[h_maxes+1]:=maxes_h[h_maxes-1];
      maxes_h[h_maxes+2]:=maxes_h[h_maxes-2];
      mins[h_mins+1]:=2*maxes[h_maxes]-mins[h_mins];
      mins[h_mins+2]:=2*maxes[h_maxes]-mins[h_mins-1];
      mins_h[h_mins+1]:=mins_h[h_mins];
      mins_h[h_mins+2]:=mins_h[h_mins-1];
    end
    else
    begin
      maxes[h_maxes+1]:=2*num-maxes[h_maxes];
      maxes[h_maxes+2]:=2*num-maxes[h_maxes-1];
      maxes_h[h_maxes+1]:=maxes_h[h_maxes];
      maxes_h[h_maxes+2]:=maxes_h[h_maxes-1];
      mins[h_mins+1]:=num;
      mins[h_mins+2]:=2*num-mins[h_mins];
      mins_h[h_mins+1]:=datas[h_datas];
      mins_h[h_mins+2]:=mins_h[h_mins];
    end;
  end;    //右延拓的end
end;


/////////////////////////////////////////////////////////////////////////////
procedure TEmd.sort(var datas:myIntArray);
var
  i,j,t:integer;
begin
  for i := low(datas) to high(datas) - 1 do
    for j := high(datas) downto i + 1 do
      if datas[i] > datas[j] then
      begin
        t := datas[i];
        datas[i] := datas[j];
        datas[j] := t;
      end;
end;


/////////////////////////////////////////////////////////////////////////////



/////////////////////////////////////////////////////////////////////////////
function TEmd.TheMax(datas:myArray):double;
var
  i:integer;
  temp:double;
begin
  temp:=datas[low(datas)];
  for i:=low(datas) to high(datas) do
    if temp < datas[i] then temp:=datas[i];
  result:=temp;
end;

end.