hilbert变换怪现象！

tangaoming

我做了一个简单的例子
fs=100;N=1024;n=0:1:N-1;dt=1/fs;t=n*dt;
x=cos(2*pi*3*t);x=hilbert(x');[fnor,t]=instfreq(x);
plot(t,fnor*fs);
发现得出来的频率在3周围波动，尤其是在端点处波动很厉害
但是如果把fs改成fs=102.4，得出来的结果就是一条值为3的水平线
然后总结出一个规律，那就是T=N/fs必须是一个整数才不会在端点处出现波动，谁能解释这是为什么吗？

yangzj

规律应该是频率3是分辨率fs/N的整数倍,也就是整周期采样.
还应该是泄露的原因

[ 本帖最后由 yangzj 于 2007-4-25 10:41 编辑 ]

zhangnan3509

如果在处理实际数据的时候，我们为了取得很好的时频图，那么就的注意这个问题了。保证要分析的频率是分辨率的整数倍？:@) 斑竹是这个意思吗？

yangzj

这个一般都没法保证吧,因为你不能确定实际的频率到底是多少.
但如果在采样的时候做一些措施,尽量做到整周期的话,也可以.像旋转类机械用键相信号来保证起止采样点.

zhangnan3509

严重同意！但是在实验室的实验当中做一些模拟故障实验，频率是多少还是基本能够保证的。这样的话实验效果会好很多。拿一些典型的旋转机械故障来说，多出现一些倍频或者分频，这样话如果我们注意到这些问题效果会很好的。对我们实验还是有作用的。你说是不是？:lol

yangzj

呵呵,那是,一般来说像我们自己做这类的实验的话都会在FFT后加校正处理的.

tangaoming

不对，当我fs=100，N=1100，3就不是fs/N的整数倍，但还是得到值为3的水平线，好像只要要求N/fs是整数就行了，那位知道为什么吗？

yangzj

呵呵,不好意思,我臆测了,对时频分析不熟,期待有关高手来解决

zhlong

3是fs/N=1/11的33倍，所以是整数倍。
再举例如下：

一、N/fs 为整数，但是有波动

fs=100;N=600;n=0:1:N-1;dt=1/fs;t=n*dt;
x=cos(2*pi*3.3*t);x=hilbert(x');[fnor,t]=instfreq(x);
plot(t,fnor*fs);

N/fs为整数，但3.3/(fs/N)=3.3/(100/600)=19.8 不是整数，这个时候出现波动。
显然不满足楼主下面这个猜想。

好像只要要求N/fs是整数就行了

二、N/fs不为整数，却没有波动

fs=100;N=1025;n=0:1:N-1;dt=1/fs;t=n*dt;
x=cos(2*pi*4*t);x=hilbert(x');[fnor,t]=instfreq(x);
plot(t,fnor*fs)

N/fs不为整数，但它的瞬时频率线却近似为直线，所以这又是楼主假设的一个反例。这里4/(fs/N)=4/(100/1025)=41是整数，满足整周期采样。


所以可以得出结论：基于信号hilbert变换和函数instfreq求取信号瞬时频率时为了避免波动，必须对信号进行整周期采样。由于HHT也是基于这一方法求信号瞬时频率，由此可知在利用EMD、HHT画信号时频图时，其中时频线的波动一部分可能是由于没有对信号进行整周期采样造成的。

三、验证非整周期采样对HHT时频图的影响

1、整周期采样

t=[0.001:0.001:1];
x=sin(2*pi*200*t);
y=sin(2*pi*80*t);
imf=[x;y];%直接把x、y当作emd分解结果
[A,f,tt]=hhspectrum(imf);
[im,tt]=toimage(A,f);
disp_hhs(im)
colormap(flipud(gray))
title('整周期采样时的HHT时频图')

fs=1000,N=1000，所以x和y都是整周期采样，得出hht时频图如下，两条时频线都是比较好的直线。


2、非整周期采样

令x=sin(2*pi*200.5*t); 此时x不是整周期采样。得出HHT时频图如下：

可以看出归一化频率0.2处的时频线两端有微小波动，将局部放大如下图：

由上可知：非整周期采样确实会造成HHT时频图时频线两端的波动，不过这种波动相对于采样频率较小，所以在频率轴显示范围为[0 , fs/2]的HHT时频图上，由非整周期采样造成的波动视觉上还是较小的。

四、待解决问题

非整周期采样在FFT分析中造成的泄漏容易理解，非整周期采样是如何造成瞬时频率曲线的波动的？理论上如何解释？

[ 本帖最后由 zhlong 于 2007-8-21 20:01 编辑 ]

zhlong

非整周期采样在FFT分析中造成的泄漏容易理解，非整周期采样是如何造成瞬时频率曲线的波动的？理论上如何解释？

想了一下这个问题，其实还是由于FFT造成的，因为hilbert()函数是调用FFT的。

http://forum.vibunion.com/forum/viewthread.php?tid=47303&extra=&highlight=%2B%D3%AA%C9%FA&page=1  
15楼  xray的回贴摘录如下：

至于Hilbert变换在matlab中的实现，其实很简单，就是把实信号fft变换，然后把负频率部分强制赋0，然后ifft，就得到了y=hilbert(x)。从这个意义上说，两端差距较大，是因为fft以后存在能量扩散，正频率部分存在从负频率扩散而来的能量，因此即使把负频率部分强制赋0，也不能完全保证消除所有的负频率能量。从滤波器的角度来说，hilbert变换实际上是通过一个冲击响应为1/t的滤波器，在前端由于最前面认为是全0，所以滤波存在误差，若采用某些方法对前端数据进行预测（例如镜像方法），可能可以消除一部分误差。

所以非整周期采样会造成hilbert求解的不正确，我们可以用下面的程序进行验证。
1. 整周期采样情况

t=[0.001:0.001:1];
x = cos(20*pi*t)+cos(40*pi*t); %两信号分量都是整周期采样
y = hilbert(x);  %hilbert函数得到的是信号的解析信号形式
x1 = imag(y);   %取hilbert的虚部，即信号的希尔伯特变换
x2 = sin(20*pi*t)+sin(40*pi*t); %这是信号的理论希尔伯特变换
figure(1);
plot(t,x1,t,x2)   % 画图比较hilbert函数计算得到的和理论上的信号希尔伯特变换

运行上面程序，得到下图，可以看到整周期采样时hilbert函数误差极小，很难观察到hilbert函数计算得到的希尔伯特变换与理论值有什么区别。

2. 非整周期采样情况

t=[0.001:0.001:1];
x = cos(20.5*pi*t)+cos(40.5*pi*t); %两信号分量都是非整周期采样
y = hilbert(x);  %hilbert函数得到的是信号的解析信号形式
x1 = imag(y);   %取hilbert的虚部，即信号的希尔伯特变换
x2 = sin(20.5*pi*t)+sin(40.5*pi*t); %这是信号的理论希尔伯特变换
figure(1);
plot(t,x1,t,x2)   % 画图比较hilbert函数计算得到的和理论上的信号希尔伯特变换

运行上面的程序，得到下图，可以看到两端差别较大，由此也可以理解非整周期采样时，楼主计算的瞬时频率两端波动较大了。


[ 本帖最后由 zhlong 于 2007-8-10 14:03 编辑 ]

tangaoming

那这个问题如何解决呢，实际采集到的数据很难做到整周期采样呀

zhlong

hilbert变换两端存在波动的现象，统称为hilbert端点效应，造成原因还有能量泄漏、gibbs现象等。可以参照emd去除端点效应的办法，将信号两端延长后再hilbert变化，之后再去掉两端部分。

form

现在想想对信号是否整周期采样有些疑问：
x=sin(2*pi*200*t);
y=sin(2*pi*80*t);
对于这两个信号，采样频率fs=1000；怎么说他们是整周期采样呢？
这两个信号的频率不是倍数关系，满足一个整周期采样，另一个就不是了吧；
由此引出满足一个信号整周期采样，另一个不满足，出现结果是一个两端不波动，一个波动；两个信号同时满足整周期采样，两端都不波动；两个信号都不满足整周期采样，两端都波动。

如果对于非周期信号来讲，是不是就不用考虑这样的影响
不知道这样理解是否正确

zhlong

现在想想对信号是否整周期采样有些疑问：
x=sin(2*pi*200*t);
y=sin(2*pi*80*t);
对于这两个信号，采样频率fs=1000；怎么说他们是整周期采样呢？
这两个信号的频率不是倍数关系，满足一个整周期采样，另一个就不是了吧；

t=[0.001:0.001:1];
x=sin(2*pi*200*t);
y=sin(2*pi*80*t);

采样点数1000，采样频率1000，所以频率分辨率为1Hz，因此200和80都是1Hz的整数倍。

zhangnan3509

N/fs=L*T,是这样计算吧 其中L取整数