希尔伯特变换的实现——数据分析漫谈

在利用希尔伯特变换提求地震资料的地震子波时，其实际的困难在于相位谱的求取，附加希尔伯特变换的程序，望以后能够认真分析，认真学习，积累点点滴滴的知识。

#define PI 3.1415926      
#define PI2 6.2831853
#include "stdio.h"
#include "math.h"
/*  inv=1 forward transform; inv=-1 inverse transform  */
fft(float sr[],float sx[],int m0,int inv)
{
 int i,j,lm,li,k,lmx,np,lix,mm2;
 float t1,t2,c,s,cv,st,ct;
 if(m0<0)
  return;
 inv=-inv;
 lmx=1;
 for(i=1;i<=m0;++i)
        lmx+=lmx;   //get 2**m0
    cv=PI2/(double)lmx;
 ct=cos(cv); st=sin(cv)*inv;
 np=lmx;mm2=m0-2;
 /* fft butterfly numeration */
 for(i=1;i<=mm2;++i)
 {
  lix=lmx;lmx/=2;
  c=ct;s=st;
  for(li=0;li<np;li+=lix)
  {
   j=li;k=j+lmx;
   t1=sr[j]-sr[k];t2=sx[j]-sx[k];
   sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=t1;sx[k]=t2;
   ++j;++k;
   t1=sr[j]-sr[k];t2=sx[j]-sx[k];
   sr[j]+=sr[k];sx[j]+=sx[k];
   sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;
  }
  for(lm=2;lm<lmx;++lm)
  {
   cv=c;c=ct*c-st*s;s=st*cv+ct*s;
   for(li=0;li<np;li+=lix)
   {
    j=li+lm;k=lmx+j;
    t1=sr[j]-sr[k];t2=sx[j]-sx[k];
    sr[j]+=sr[k];sx[j]+=sx[k];
    sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;
   }
  }
  cv=ct;ct=2.0*ct*ct-1.0;st=2.0*st*cv;
 }
 /* 4 points DFT */
 if(m0>=2)
 for(li=0;li<np;li+=4)
 {
  j=li;k=j+2;
  t1=sr[j]-sr[k];t2=sx[j]-sx[k];
  sr[j]+=sr[k];
  sx[j]+=sx[k];
  sr[k]=t1;sx[k]=t2;
  ++j;++k;
  t1=sr[j]-sr[k];t2=sx[j]-sx[k];
  sr[j]+=sr[k];sx[j]+=sx[k];
  sr[k]=-inv*t2;sx[k]=inv*t1;
 }
 /* 2 points DFT */
 for(li=0;li<np;li+=2)
 {
  j=li;k=j+1;
  t1=sr[j]-sr[k];t2=sx[j]-sx[k];
  sr[j]+=sr[k];sx[j]+=sx[k];
  sr[k]=t1;sx[k]=t2;
 }
 /* sort according to bit reversal */
 lmx=np/2;j=0;
 for(i=1;i<np-1;++i)
 {
  k=lmx;
  while(k<=j)
  {
   j-=k;k/=2;
  }
  j+=k;
  if(i<j)
  {
   t1=sr[j];sr[j]=sr[i];sr[i]=t1;
   t1=sx[j];sx[j]=sx[i];sx[i]=t1;
  }
 }
 /*  if Inverse FFT, multiply 1.0/np  */
 if(inv!=1.0)
  return;
 t1=1.0/np;
 for(i=0;i<np;++i)
 {
  sr[i]*=t1;sx[i]*=t1;
 }
}

main()
{
 float xr[2048],xi[2048],xt[2048];
  int i,np,nfft,k,nf;
  float t,dt,df,f,hf;
  float f0=10,f1=20,f2=30;
  FILE *fp1,*fp2;
  char fil1[60],fil2[60];
  printf("ENTER TIME SIGNAL FILE\n");
  scanf("%s",fil1);
  printf("ENTER AMPLITUDE SPECTRUM FILE\n");
  scanf("%s",fil2);
  printf("SAMPLE POINTS?\n");
  scanf("%d",&np);
  printf("SAMPLE INTERVAL(S)?\n");
  scanf("%f",&dt);
  printf("HIGH CUTOFF FREQUENCY(Hz)?\n");
  scanf("%f",&hf);
  fp1=fopen(fil1,"w");
  fp2=fopen(fil2,"w");
  for(i=0;i<np;i++)
 { t=(i-np/2)*dt;
   if(t!=0.0)xt[i]=1.0/(PI*t);
   else xt[i]=0.0;
   xr[i]=xt[i];
 }
  // calculate fft point
   k=log(np*1.0)/log(2.0);
 if(np>pow(2.0,k*1.0))k=k+1;
 nfft=pow(2.0,k*1.0);
 df=1.0/(nfft*dt);
 nf=hf/df+1;
 printf("nfft=%d  k=%d\n",nfft,k);
 printf("dt=%f  df=%f  nf=%d\n",dt,df,nf);
  // fill zero
 if(np<nfft)
  {for(i=np;i<nfft;i++)
   xr[i]=0;
  }
 for(i=0;i<nfft;i++)
  xi[i]=0.0;
  // calculate fft
 fft(xr,xi,k,1);
  // output amplitude and argument
 for(i=0;i<nf;i++)
  { f=i*df;
    fprintf(fp2,"%d %8.2f %12.4f %12.4f\n",i,f,atan(xr[i]/xi[i]),sqrt(xr[i]*xr[i]+xi[i]*xi[i]));
  }
 fft(xr,xi,k,-1);
 for(i=0;i<np;i++)
  { t=i*dt;
    fprintf(fp1,"%10d %10.4f %10.4f %10.4f\n",i+1,t,xt[i],xr[i]);
  }
 fclose(fp1);
 fclose(fp2);
}

 

希尔伯特变换c程序例子

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"FFT1.h"
const int L=4;//信号长度
//x初始序列，也存变换结果序列
void HT(float x[],int num)
{
    float *xi,*x1r, *x1i, *sgn;
    int i,k;
    k=log(num)/log(2);
                     xi=(float *)calloc(num,sizeof(float));
    x1r=(float *)calloc(num,sizeof(float));
    x1i=(float *)calloc(num,sizeof(float));
    sgn=(float *)calloc(num,sizeof(float));
    //apply FFT to x[]
    fft(x,xi,k,1);
                     //construct sgn[]
    for( i=1; i<num/2; i++)
        sgn=1.0;
    for(i=num/2+1; i<num; i++)
        sgn=-1.0;
    //X(k)sgn(k)
    for(i=0; i<num; i++)
    {
        x1r=x*sgn;
        x1i=xi*sgn;
    }
    //apply IFFT to x1
    fft(x1r,x1i,k,-1);
    for(i=0; i<num; i++)
          x= x1i;
                      free(x1r);
    free(x1i);
    free(sgn);
    free(xi);
}
void main()
{
    float x[L]={1,2,3,4};
    float xx[L];
    int i;
    for(i=0; i<L; i++)
            xx=x;
    HT(x,L);
    printf("Hilbert变换得到的解析信号为:\n");
    for(i=0; i<L; i++)
             printf("%f+i(%f)\n",xx,x);
}