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

Pascal（Delphi7）写的一个THHT类

{Author:Wang qiang
 Current version:1.0 (20/4/2008)


 E-Mail:[e-mail]
 This code is translated from Matlat code by Rilling
 The TEmd include Empirical Mode Decompsition and Hilbert Transform
 }


unit EMD;

interface
uses
  math,dialogs,series;

const                   //常数
  sd=0.05;
  sd2=0.5;
  tol=0.05;
  MAXMODES=20;
  MAXITERATIONS=2000;

type
  myDArray=array of array of double;
  myArray=array of double;
  myIntArray=array of integer;

type
  TEmd=class (TObject)
  private

    data:myArray;
    nDatas:integer;         //信号长度
    imf:myDArray;
    nImf:integer;
    hReal:myDArray;
    hImaginary:myDArray;
    a:myDArray;
    omg:myDArray;
    spotMatrix:myDArray;
    marginalMatrix:myArray;

    function TheMax(datas:myArray):double;                                               //返回信号中最大值
    procedure BoundaryExpand(const datas:myArray;var mins,maxes:myIntArray;var mins_h,maxes_h:myArray);
    procedure sort(var datas:myIntArray);
    procedure Spline(X:myIntArray;Y:myArray;N:integer;YP1,YPN:double;var Y2:myArray);   //样条插值函数
    procedure Splint(XA:myIntArray;YA,Y2A:myarray;N:integer;X:integer;var Y:double);//样条插值函数
    procedure extr(const datas:myArray;var indmin,indmax,indzero:myIntArray);           //得到数据极值点和过零点的索引
    procedure MeanAndAmplitude(const m:myarray;var envmoy,amp:myarray;var nem,nzm:integer);   //得到包络线和均值
    procedure stopSifting(const m:myarray;var stop:boolean;var envmoy:myarray);//;var s:double);     //IMF的判定函数
    function stopEMD(const datas:myArray):boolean;   //是否停止筛选

    procedure Four1(var Datas:myArray;Isign:integer);  //傅立叶变换
    procedure unwrap(var Datas:myArray);

    procedure GetSpectrum(const deltaT:double);
    procedure hilbert();
    procedure HilbertSpectrum(const deltaT:double);

  public
    property sigal:myArray read data;       //sigal属性，读出信号
    property IMFS:myDArray read imf;        //IMFS属性，读出imf
    property NImfs:integer read nImf;
    property hilbertSpotMatrix:myDArray read SpotMatrix;  //读出Hilbert谱
    property marginarySpotMatrix:myArray read marginalMatrix;

    constructor Create(const datas:myArray);
    destructor Destroy;override;  //?

    function emd():myDArray;               //进行emd筛选计算
    function SpotSpectrum(const XMesh,YMesh:integer;const deltaT:double):myDArray;
    function MarginalSpectrum():myArray;

  end;


implementation




///////////////////////////////////////////////////////////////////////
constructor TEmd.Create(const datas:myArray);
begin
  data:=copy(datas);
  
end;



////////////////////////////////////////////////////////////////////////
destructor TEmd.Destroy;
begin
  inherited;
end;
////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////
function TEmd.MarginalSpectrum():myArray;
var
  i,j:integer;
  m,n:integer;
  s:double;
begin
  m:=length(SpotMatrix);
  n:=length(SpotMatrix[0]);
  setlength(MarginalMatrix,n);
  s:=0;
  for i:=0 to n-1 do
  begin
    for j:=0 to m-1 do
      s:=s+sqr(SpotMatrix[j,i]);
    s:=s/m;
    MarginalMatrix[i]:=s;
  end;//for i...

  result:=copy(MarginalMatrix);
end;
///////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////////////////
function TEmd.SpotSpectrum(const XMesh,YMesh:integer;const deltaT:double):myDArray;
var
  //m,n:integer; //m-imf个数；n-信号长度
  omgMax,omgMin:double;
  i,j:integer;
  tempTime,tempFreq:integer;
  //a,omg:myDArray;
begin
  //time:=nil;freq:=nil;
  SpotMatrix:=nil;
  getSpectrum(deltaT);       //计算a,omg

  //setlength(time,XMESH+1);setlength(freq,YMESH+1);
  setlength(SpotMatrix,XMESH+1,YMESH+1);

  //m:=length(a);
  //n:=length(a[0]);

  omgMax:=omg[0,0];
  omgMin:=omg[0,0];

  for i:=0 to nImf-1 do
    for j:=0 to nDatas-1 do
    begin
      //if aMax<a[i,j] then aMax:=a[i,j];
      //if aMin>a[i,j] then aMin:=a[i,j];
      if omgMax<omg[i,j] then omgMax:=omg[i,j];
      if omgMin>omg[i,j] then omgMin:=omg[i,j];
    end;

  for i:=0 to nImf-1 do
    for j:=0 to nDatas-1 do
    begin
       tempTime:=round(XMESH*j/(nDatas-1));
       tempFreq:=round(YMESH*(omg[i,j]-omgMin)/(omgMax-omgMin));
       SpotMatrix[tempTime,tempFreq]:=a[i,j];
    end;

  result:=copy(spotMatrix);

  //for i:=0 to XMESH do time[i]:=round(i/XMESH*(n-1));
  //for i:=0 to YMESH do freq[i]:=round(i/YMESH*(omgMax-omgMin)+omgMin);
end;
///////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////
procedure TEmd.GetSpectrum(const deltaT:double);
begin
  hilbert();
  Hilbertspectrum(DeltaT);
end;
////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////
procedure TEmd.Hilbertspectrum(const deltaT:double);
var
  ang:myarray;
  n,i,k:integer;
begin
  a:=nil;
  omg:=nil;
  setlength(a,length(hReal));
  setlength(omg,length(hImaginary));

for k:=low(hReal) to high(hImaginary) do
begin
  n:=length(hReal[k]);
  setlength(ang,n);
  setlength(a[k],n);
  setlength(omg[k],n);
  for i:= 0 to n-1 do
  begin
    a[k,i]:=sqrt(sqr(hReal[k,i])+sqr(hImaginary[k,i]));
    ang[i]:=arctan2(hReal[k,i],hImaginary[k,i]);
  end;
  unwrap(ang);
  for i:=0 to n-2 do omg[k,i]:=(ang[i+1]-ang[i])/(2*pi*DeltaT);       //最后一个数为0
  omg[k,n-1]:=omg[k,n-2];

end;
end;
////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////
procedure TEmd.unwrap(var Datas:myArray);
var
  DataOut:myArray;
  n,i:integer;
begin
  n:=length(Datas);
  setlength(DataOut,n);
  DataOut[0]:=Datas[0];
  for i:=1 to n-1 do
  if Datas[i]>Datas[i-1] then
    DataOut[i]:=DataOut[i-1]+(Datas[i]-Datas[i-1])
  else
    DataOut[i]:=DataOut[i-1]+2*pi+(Datas[i]-Datas[i-1]);
  //for i:=0 to n-1 do Datas[i]:=DataOut[i];
  Datas:=copy(DataOut);
end;
////////////////////////////////////////////////////////////////////////


////////////////////////////////////////////////////////////////////////
procedure TEmd.Four1(var Datas:myArray;Isign:integer);
var
    N,J,JJ,I,M,II,ISTEP,MMAX:integer;
    THETA,WR,WI,WPR,WTEMP,TEMPR,TEMPI,WPI:real;
begin
    N:=2 * nDatas;
    J:=1;
    For II:=1 To nDatas do
    begin
        I:= II * 2 - 1;
        If J > I Then
        begin
            TEMPR:=DATAS[J];
            TEMPI:=DATAS[J + 1];
            DATAS[J]:=DATAS[I];
            DATAS[J + 1]:=DATAS[I + 1];
            DATAS[I]:=TEMPR;
            DATAS[I + 1]:=TEMPI;
        end; 
        M:=N div 2;
        While (M >= 2) And (J > M) do
        begin
            J:=J - M;
            M:=M div 2;
        end;
        J:=J + M;
    end;
    MMAX:=2;
    While N > MMAX do
    begin
        ISTEP:=2 * MMAX;
        THETA:=6.28318530717959 / (ISIGN * MMAX);
        WPR:=-2 * Sqr(Sin(0.5 * THETA));
        WPI:=Sin(THETA);
        WR:=1;
        WI:=0;
        For II:=1 To (MMAX div 2) do
        begin
            M:= II * 2 - 1;
            For JJ:=0 To ((N-M) div ISTEP) do
            begin
                I:= M + JJ*ISTEP ;
                J:=I + MMAX;
                TEMPR:=WR * DATAS[J] - WI * DATAS[J + 1];
                TEMPI:=WR * DATAS[J + 1] + WI * DATAS[J];
                DATAS[J]:=DATAS[I] - TEMPR;
                DATAS[J + 1]:=DATAS[I + 1] - TEMPI;
                DATAS[I]:=DATAS[I] + TEMPR;
                DATAS[I + 1]:=DATAS[I + 1] + TEMPI;
            end; 
            WTEMP:=WR;
            WR:=WR * WPR - WI * WPI + WR;
            WI:=WI * WPR + WTEMP * WPI + WI;
        end; 
        MMAX:=ISTEP;
    end;
end;
////////////////////////////////////////////////////////////////////////




////////////////////////////////////////////////////////////////////////
procedure TEmd.hilbert();
var
  x,h:myarray;
  i,j,k,isign:integer;
begin
  hReal:=nil;
  hImaginary:=nil;//初始化
  setlength(hReal,length(imf));
  setlength(hImaginary,length(imf));


//for k:=low(imf) to high(imf) do
for k:=0 to high(imf) do
begin
  //n:=length(imf[k]);
  setlength(hReal[k],nDatas);
  setlength(hImaginary[k],nDatas);
  setlength(h,nDatas);
  for i:=0 to nDatas-1 do h[i]:=0;
  setlength(x,2*nDatas+1);
  x[0]:=0;
  for i:=0 to nDatas-1 do
  begin
      x[2*i+1]:=imf[k,i];
      x[2*i+2]:=0;
  end;
  isign:=1;
  four1(x,isign);
  if (nDatas mod 2)=0 then
  begin
    j:=nDatas div 2;
    h[0]:=1;h[j]:=1;
    for i:=1 to j-1 do h[i]:=2;
  end
  else
  begin
    j:=((nDatas+1) div 2)-1;
    h[0]:=1;
    for i:=1 to j do h[i]:=2;
  end;
  for i:=0 to nDatas-1 do
  begin
   x[2*i+1]:=h[i]*x[2*i+1];
   x[2*i+2]:=h[i]*x[2*i+2];
  end;
  isign:=-1;
  four1(x,isign);
  for i:=0 to nDatas-1 do
  begin
   hReal[k,i]:=x[2*i+1]/nDatas;
   hImaginary[k,i]:=x[2*i+2]/nDatas;
  end;

end;//for k:=
end;
////////////////////////////////////////////////////////////////////////


////////////////////////////////////////////////////////////////////////
function TEmd.emd():myDArray;
var
  m,r,mp,absM,absDatas:myarray;
 // nDatas:integer;         //原始数据长度
  i,n:integer;
  stopsift:boolean;
  moyenne:myarray;
  //s:double;
begin

  nDatas:=length(data);
  setlength(r,nDatas);setlength(m,nDatas);setlength(mp,nDatas);
  setlength(absM,nDatas);setlength(absDatas,nDatas);        //初始化

  r:=copy(data);
  for i:=0 to nDatas-1 do absDatas[i]:=abs(data[i]);


  nImf:=0;//imf的个数
  while not (stopEMD(r)) do
  begin


    for i:=0 to nDatas-1 do      //m,mp,
    begin
      m[i]:=r[i];
      mp[i]:=m[i];
      absM[i]:=abs(m[i]);
      //absDatas[i]:=abs(datas[i]);
    end;


    stopsifting(m,stopsift,moyenne);//,s);              //**

    if theMax(absM)<(1e-10)*theMax(absDatas) then
      break;


    n:=0;
    while not stopsift do       //找m
    begin

      for i:=low(m) to high(m) do m[i]:=m[i]-moyenne[i];


      stopsifting(m,stopsift,moyenne);//,s);
      inc(n);
      if n>MAXITERATIONS then break;

    end;
    inc(nImf);
    setlength(imf,nImf,length(m));
    for i:=low(m) to high(m) do imf[nImf-1,i]:=m[i];

    for i:=low(r) to high(r) do r[i]:=r[i]-m[i];

    if nImf>MAXMODES then break;
  end; // while not (stopEMD(datas)) 大循环

  result:=copy(imf);
end;


////////////////////////////////////////////////////////////////////////////
function TEmd.stopEMD(const datas:myArray):boolean;
var
  indmax,indmin,indzero:myIntArray;
begin
  extr(datas,indmin,indmax,indzero);
  if (length(indmin)+length(indmax))<20 then stopEMD:=true else stopEMD:=false;
end;


//////////////////////////////////////////////////////////////////////////////////////////
procedure TEmd.StopSifting(const m:myarray;var stop:boolean;var envmoy:myarray);//;var s:double);
var
  amp:myarray;
  nem,nzm:integer;
  sx:myarray;
  n:integer;
  temp,temp1:double;
  keepon1,keepon2,keepon3,keepon4:boolean;
  i:integer;
begin
try
  envmoy:=nil;
  //amp:=nil;
  MeanAndAmplitude(m,envmoy,amp,nem,nzm);       //*
  n:=length(envmoy);
  setlength(sx,n);
  //s:=0;
  temp:=0;
  temp1:=0;
  for i:=0 to n-1 do
  begin
    sx[i]:=abs(envmoy[i])/amp[i];
    if sx[i]>sd then temp:=temp+1;
    if sx[i]>sd2 then temp1:=temp1+1;
    //s:=s+sx[i];
  end;
  //s:=s/n;
  temp:=temp/n;

  if temp>tol then keepon1:=true else keepon1:=false;
  if temp1>0 then keepon2:=true else keepon2:=false;
  if nem>2 then keepon3:=true else keepon3:=false;
  if abs(nzm-nem)>1 then keepon4:=true else keepon4:=false;
  stop:= not (keepon1 or keepon2 and keepon3 and keepon4);
except
  stop:=true;
  n:=length(m);
  setlength(envmoy,n);
  for i:=0 to n-1 do envmoy[i]:=0;
end;//try

end;

/////////////////////////////////////////////////////////////////////////////////////////////
procedure TEmd.MeanAndAmplitude(const m:myarray;var envmoy,amp:myarray;var nem,nzm:integer);
var
  indmin,indmax,indzero:myIntArray;
  tmin,tmax:myIntArray;
  mmin,mmax:myarray;
  nMin,nMax:integer;
  i,n:integer;
  p1Min,pnMin,p1Max,pnMax:double;
  y2Min,y2Max:myarray;
  tempEnvMin,tempEnvMax:myarray;
  ind1,ind2:integer;
begin
  envmoy:=nil;
  amp:=nil;         //初始化

  extr(m,indmin,indmax,indzero);                 //**
  nzm:=length(indzero);                          //过零点个数
  nMin:=length(indmin);
  nMax:=length(indmax);
  nem:=nMin+nMax;            //极值点个数

  nMin:=nMin+4;
  nMax:=nMax+4;
  setlength(tmin,nMin);
  setlength(mmin,nMin);
  setlength(tmax,nMax);
  setlength(mmax,nMax);

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