在利用希尔伯特变换提求地震资料的地震子波时,其实际的困难在于相位谱的求取,附加希尔伯特变换的程序,望以后能够认真分析,认真学习,积累点点滴滴的知识。
#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); } |