最近毕设需要用到希尔伯特变换的知识,今天做完之后决定还是记录下思路。
当然是数字信号的希尔伯特变换

上面是连续信号的希尔伯特变换,离散的应该也能根据上面写(没现成的图片,懒得编辑公式了)。
这里打算采用使用卷积的方法来计算。由于希尔伯特变换的传输函数的傅里叶变换是 H(w)= -j w>=0
j w<0
所以我们可以先求原始信号的离散傅里叶变换E(K),然后按照下面的公式就可以求出希尔伯特变换之后的信号的离散傅里叶变换Eh(K)。
然后对Eh(K)求反傅里叶变换就可得到我们需要的信号的希尔伯特变换信号。
下面贴代码思路
先建立一个复信号的结构体:
typedef struct {
Float64 r; /* 实部 */
Float64 i; /* 虚部 */
} CPX;
接着是离散傅里叶变换的函数 第一个参数dir代表正变换和反变换。
void DFT(int dir,int framelen,CPX *signal,CPX *dft_s)
{
int i,k;
double arg;
double cosarg,sinarg;
for(i=0;i<framelen;i++)
{
arg=-dir*2.0*3.141592654*(double)i/(double)framelen;
for(k=0;k<framelen;k++)
{
cosarg=cos(k*arg);
sinarg=sin(k*arg);
dft_s[i].r+=(signal[k].r*cosarg-signal[k].i*sinarg);
dft_s[i].i+=(signal[k].r*sinarg+signal[k].i*cosarg);
}
}
/*返回数据*/
if(dir==-1)
{
for(i=0;i<framelen;i++)
{
dft_s[i].r=dft_s[i].r/(double)framelen;
dft_s[i].i=dft_s[i].i/(double)framelen;
}
}
}
上主函数
int main(void)
{
CPX *signal;
CPX *dft_s;
CPX *hdft_s;
CPX *hsignal;
signal=calloc(framelen,sizeof(CPX)); // 原始信号
dft_s=calloc(framelen,sizeof(CPX)); // 原始信号的傅里叶变换
hdft_s=calloc(framelen,sizeof(CPX)); // 希尔伯特变换的离散 傅里叶变换
hsignal=calloc(framelen,sizeof(CPX)); // 希尔伯特变换后信号
DFT(1,framelen,signal, dft_s); //求原始信号 傅里叶变换
for(i=0;i<framelen;i++) //求出希尔伯特变换信号的傅里叶变换
{
if(i<=framelen/2)
{
hdft_s[i].r=dft_s[i].i;
hdft_s[i].i=-dft_s[i].r;
}
else
{
hdft_s[i].r=-dft_s[i].i;
hdft_s[i].i=dft_s[i].r;
}
}
DFT(-1,framelen, hdft_s,hsignal); //利用反傅里叶变换求出希尔伯特变换信号
}
最后这个hsignal就是我们要的希尔伯特变换信号了。
|