请问如何用DFT实现希尔伯特变换？Matlab里面是怎么实现的？

请问如何用DFT实现希尔伯特变换？Matlab里面是怎么实现的？

#define PI 3.1415926
#define PI2 6.2831853
#include "stdio.h"
#include "math.h"

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;

for(i=1;i<=mm2;++i)
{
lix=lmx;lmx/=2;
c=ct;s=st;
for(li=0;li{
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{
cv=c;c=ct*c-st*s;s=st*cv+ct*s;
for(li=0;li {
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;
}

if(m0>=2)
for(li=0;li {
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;
}

for(li=0;li {
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;
}

lmx=np/2;j=0;
for(i=1;i {
k=lmx;
while(k<=j)
{
j-=k;k/=2;
}
j+=k;
if(i{
t1=sr[j];sr[j]=sr[i];sr[i]=t1;
t1=sx[j];sx[j]=sx[i];sx[i]=t1;
}
}

if(inv!=1.0)
return;
t1=1.0/np;
for(i=0;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 { 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=%dk=%d\n",nfft,k);
printf("dt=%fdf=%fnf=%d\n",dt,df,nf);
// fill zero
if(np{for(i=np;i xr[i]=0;
}
for(i=0;ixi[i]=0.0;
// calculate fft
fft(xr,xi,k,1);
// output amplitude and argument
for(i=0;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{ t=i*dt;
fprintf(fp1,"%10d %10.4f %10.4f %10.4f\n",i+1,t,xt[i],xr[i]);
}
fclose(fp1);
fclose(fp2);
}