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

#define PI 3.1415926

#define PI2 6.2831853

#include "stdio.h"

#include "math.h"

/* inv=1 forward transforminv=-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=1i<=m0++i)

lmx+=lmx //get 2**m0

cv=PI2/(double)lmx

ct=cos(cv)st=sin(cv)*inv

np=lmxmm2=m0-2

/* fft butterfly numeration */

for(i=1i<=mm2++i)

{

lix=lmxlmx/=2

c=cts=st

for(li=0li<npli+=lix)

{

j=lik=j+lmx

t1=sr[j]-sr[k]t2=sx[j]-sx[k]

sr[j]+=sr[k]sx[j]+=sx[k]sr[k]=t1sx[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*t2sx[k]=s*t1+c*t2

}

for(lm=2lm<lmx++lm)

{

cv=cc=ct*c-st*ss=st*cv+ct*s

for(li=0li<npli+=lix)

{

j=li+lmk=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*t2sx[k]=s*t1+c*t2

}

}

cv=ctct=2.0*ct*ct-1.0st=2.0*st*cv

}

/* 4 points DFT */

if(m0>=2)

for(li=0li<npli+=4)

{

j=lik=j+2

t1=sr[j]-sr[k]t2=sx[j]-sx[k]

sr[j]+=sr[k]

sx[j]+=sx[k]

sr[k]=t1sx[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*t2sx[k]=inv*t1

}

/* 2 points DFT */

for(li=0li<npli+=2)

{

j=lik=j+1

t1=sr[j]-sr[k]t2=sx[j]-sx[k]

sr[j]+=sr[k]sx[j]+=sx[k]

sr[k]=t1sx[k]=t2

}

/* sort according to bit reversal */

lmx=np/2j=0

for(i=1i<np-1++i)

{

k=lmx

while(k<=j)

{

j-=kk/=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=0i<np++i)

{

sr[i]*=t1sx[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=0i<npi++)

{ 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=npi<nffti++)

xr[i]=0

}

for(i=0i<nffti++)

xi[i]=0.0

// calculate fft

fft(xr,xi,k,1)

// output amplitude and argument

for(i=0i<nfi++)

{ 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=0i<npi++)

{ t=i*dt

fprintf(fp1,"%10d %10.4f %10.4f %10.4f\n",i+1,t,xt[i],xr[i])

}

fclose(fp1)

fclose(fp2)

}

#include <math.h>

#include <stdio.h>

#define N 8

void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il)

void main()

{

double xr[N],xi[N],Yr[N],Yi[N],l=0,il=0

int i,j,n=N,k=3

for(i=0i<Ni++)

{

xr[i]=i

xi[i]=0

}

printf("------FFT------\n")

l=0

kkfft(xr,xi,n,k,Yr,Yi,l,il)

for(i=0i<Ni++)

{

printf("%-11lf + j* %-11lf\n",Yr[i],Yi[i])

}

printf("-----DFFT-------\n")

l=1

kkfft(Yr,Yi,n,k,xr,xi,l,il)

for(i=0i<Ni++)

{

printf("%-11lf + j* %-11lf\n",xr[i],xi[i])

}

getch()

}

void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il)

{

int it,m,is,i,j,nv,l0

double p,q,s,vr,vi,poddr,poddi

for (it=0it<=n-1it++)

{

m = it

is = 0

for(i=0i<=k-1i++)

{

j = m/2

is = 2*is+(m-2*j)

m = j

}

fr[it] = pr[is]

fi[it] = pi[is]

}

pr[0] = 1.0

pi[0] = 0.0

p = 6.283185306/(1.0*n)

pr[1] = cos(p)

pi[1] = -sin(p)

if (l!=0)

pi[1]=-pi[1]

for (i=2i<=n-1i++)

{

p = pr[i-1]*pr[1]

q = pi[i-1]*pi[1]

s = (pr[i-1]+pi[i-1])*(pr[1]+pi[1])

pr[i] = p-q

pi[i] = s-p-q

}

for (it=0it<=n-2it=it+2)

{

vr = fr[it]

vi = fi[it]

fr[it] = vr+fr[it+1]

fi[it] = vi+fi[it+1]

fr[it+1] = vr-fr[it+1]

fi[it+1] = vi-fi[it+1]

}

m = n/2

nv = 2

for (l0=k-2l0>=0l0--)

{

m = m/2

nv = 2*nv

for(it=0it<=(m-1)*nvit=it+nv)

for (j=0j<=(nv/2)-1j++)

{

p = pr[m*j]*fr[it+j+nv/2]

q = pi[m*j]*fi[it+j+nv/2]

s = pr[m*j]+pi[m*j]

s = s*(fr[it+j+nv/2]+fi[it+j+nv/2])

poddr = p-q

poddi = s-p-q

fr[it+j+nv/2] = fr[it+j]-poddr

fi[it+j+nv/2] = fi[it+j]-poddi

fr[it+j] = fr[it+j]+poddr

fi[it+j] = fi[it+j]+poddi

}

}

/*逆傅立叶变换*/

if(l!=0)

{

for(i=0i<=n-1i++)

{

fr[i] = fr[i]/(1.0*n)

fi[i] = fi[i]/(1.0*n)

}

}

/*是否计算模和相角*/

if(il!=0)

{

for(i=0i<=n-1i++)

{

pr[i] = sqrt(fr[i]*fr[i]+fi[i]*fi[i])

if(fabs(fr[i])<0.000001*fabs(fi[i]))

{

if ((fi[i]*fr[i])>0)

pi[i] = 90.0

else

pi[i] = -90.0

}

else

pi[i] = atan(fi[i]/fr[i])*360.0/6.283185306

}

}

return

}

希尔伯特变换的意义本文不提，本文的目标是举例说明到底离散傅里叶变换是如何实现的，并编写程序与matlab自带的hilbert函数结果对比，验证我们的实现过程是否正确。

原始信号为 ，其离散希尔伯特变换的定义公式为：

说明：先让原始信号与 信号做卷积，然后一起合并成一个 复数 信号 。

问题：看上去很简单，但这里的卷积不是一般意义上的卷积的 *** 作！

所以：实际中得到 的方法是通过借助 离散傅里叶变换DFT 来实现的！

因此：本文就用 离散傅里叶变换 来实现 离散希尔伯特变换 。

设原始信号为 ，总长度一般 (下标n是从 1 开始到N)，总体实现步骤可分为3步：

相应的matlab程序：

结果：

手动实现正确！

希尔伯特变换的结果是给原始信号 提供了一个幅值、频率不变，但相位平移90°的信号 。

所以，希尔伯特变换是从"时域"到"时域"的变换！只改变了相位，所以又叫90°移相滤波器；

所以，原始信号 与它的希尔伯特变换 构成正交副。

现在，我们换回到最初的记法：原始信号和它对应的希尔伯特变换信号分别用 和 表示，那么对应的" 解析信号 "就可以用这两个东西组成：

对于这个解析信号，我们可以得它的3瞬属性：瞬时振幅、瞬时相位、瞬时频率。

可以看出，3瞬属性是相互关联的！

瞬时振幅、瞬时相位可以直接求，有意义；

但是瞬时频率 直接 根据解析信号这么按公式 求 ，是 没有物理意义 的！

并且 离散信号 ，它的瞬时频率求导只能按照" 差分 "来近似，即：

给出matlab的实现程序：

效果：

3瞬属性中的瞬时频率，很明显可以看出它有很多的" 负频率 "！这很明显是错误的。

所以，直接根据" 解析信号 "算瞬时频率是无意义的！

所以，真正做 3瞬属性 的分析，做原信号的" 时频谱 "分析，我们用的：

—— 希尔伯特-黄变换。

---