请问如何用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)

}

n = 0:40

a = 2b = -3

x1 = cos(2*pi*0.1*n)

x2 = cos(2*pi*0.4*n)

x = a*x1 + b*x2

num = [2.2403 2.4908 2.2403]

den = [1 -0.4 0.75]

ic = [0 0]

y1 = filter(num,den,x1,ic)% Compute the output y1[n]

y2 = filter(num,den,x2,ic)% Compute the output y2[n]

y = filter(num,den,x,ic)% Compute the output y[n]

yt = a*y1 + b*y2

d = y - yt

subplot(3,1,1)

stem(n,y)

ylabel('Amplitude')

title('Output Due to Weighted Input: a \cdot x_{1}[n] + b \cdot x_{2}[n]')

subplot(3,1,2)

stem(n,yt)

ylabel('Amplitude')

title('Weighted Output: a \cdot y_{1}[n] + b \cdot y_{2}[n]')

subplot(3,1,3)

stem(n,d)

xlabel('Time index n')ylabel('Amplitude')

title('Difference Signal')

作为学习HHT的的第二部分，第一部分emd分解可参见链接：

用matlab进行振动波形的emd分解

希尔伯特变换的物理意义十分简单： 把信号的所有频率分量的相位推迟90度。 因此又叫90°移相器，所以原始信号与它的希尔伯特变换构成正交副。

当然，我知道大家最感兴趣的是：把相位推迟90度有什么用？

答案是： 希尔伯特变换可以用来做解调器，调幅、调频都能解。

我们构造一个信号 z(t)=x(t)+i*y(t)，将该图像在三维空间中画出来，如图所示

补充：为什么通过瞬时相位求导可以定义为瞬时频率：从信号投影来看可以建立时间t和一个角度的极坐标方程，所以单位时间角度的变化就是角速度，而角速度与频率成倒数。

这样，我们就利用希尔伯特变换从一个幅度、频率均被调制的调制波中把幅度、频率都解调了出来。

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

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

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

—— 希尔伯特-黄变换（HHT）。HHT变换先将信号进行EMD分解，得到的是各个不同尺度的分量，对每一个分量进行Hilbert变换后得到的是有实际意义的瞬时频率。

举例如下：

希尔伯特-黄变换最初的理论是采用emd的经验模态分解，目前已经改进到采用ceemdan的模态分解方式

【1】https://www.zhihu.com/question/30372795   希尔伯特变换将信号表示为复解析信号的物理意义是什么？

【2】https://www.jianshu.com/p/b591d95ae80b    一维离散希尔伯特变换实现与3瞬属性

【3】https://www.jianshu.com/p/3363abb64f32    离散数据希尔伯特-黄变换

【4】https://blog.csdn.net/yrlgg/article/details/79595859    傅里叶变换与希尔伯特变换

---