小波与hilbert变换进行信号包络谱分析
这是我自己搞的一个仿真信号用小波与hilbert变换进行信号包络谱分析的程序大家看一下,给点意见啊。我是结合书上的编的,有几个地方一直不懂,
1,许多论文上说hilbert变换是这样进行的:原始信号X(t)带通滤波,进行hilbert变换得到x^(t),作为虚部, 然后用w(t)=x(t)+jx^(t)作为解析信号。包络谱的求法是:a(t)=X(t)和x^(t)的平方和再开放.但是这个程序里面怎么没有这个功能??
2,plot((0:nfft/2-1)/nfft*fs,p(1:nfft/2));这句程序有点难懂,麻烦哪位编程高手解释一下
t=0:0.005:1*pi;
fs=10000;
s=4*sin(2*200*pi*t).*(sin(2*4500*pi*t))+25*(sin(2*4500*pi*t));
subplot(411);plot(t,s)
=wavedec(s,1,'db10');
d1=wrcoef('d',c,l,'db10');
a1=0;
subplot(412);plot(d1);title('重构高频信号');
y=hilbert(d1);
y1=abs(y);
ydata=y-mean(y);
nfft=1024;
p=abs(fft(y1,nfft));
figure(2);
plot((0:nfft/2-1)/nfft*fs,p(1:nfft/2));
xlabel('频率');
ylabel('功率谱'); 1,许多论文上说hilbert变换是这样进行的:原始信号X(t)带通滤波,进行hilbert变换得到x^(t),作为虚部,然后用w(t)=x(t)+jx^(t)作为解析信号。包络谱的求法是:a(t)=X(t)和x^(t)的平方和再开放.但是这个程序里面怎么没有这个功能??
你的程序中小波分解实现了上面说的滤波功能,而hilbert这个函数就是求上面说的w(t).
2,plot((0:nfft/2-1)/nfft*fs,p(1:nfft/2));这句程序有点难懂,麻烦哪位编程高手解释一下
这个可以看一下有关“如何正确得到Fourier分析幅值”的帖子。 我觉得程序应该是这样的
s=data;
plot(s);
hold on
=wavedec(s,2,'db10');
d1=wrcoef('d',c,l,'db10');
y=hilbert(d1);
y1=imag(y).^2;
w=sqrt(d1.^2+y1)
y2=abs(w);
plot(w,'r-')
figure(2);
plot(y2);
p=abs(fft(y2,10000));
figure(3);
plot(p);
xlabel('频率');
ylabel('功率谱');
不知大家有何看法??
是对还是错?? 本帖最后由 VibInfo 于 2016-10-21 15:24 编辑
原帖由 wxk8000 于 2008-1-9 16:05 发表
原始信号X(t)带通滤波,进行hilbert变换得到x^(t),作为虚部,
然后用w(t)=x(t)+jx^(t)作为解析信号。包络谱的求法是:a(t)=X(t)和x^(t)的平方和再开方。
楼主程序中的这两语句
y=hilbert(d1);
y1=abs(y);
便是完成把滤波后的信号进行hilbert变换得到x^(t),作为虚部,然后用w(t)=x(t)+jx^(t)作为解析信号。包络谱的求法是:a(t)=X(t)和x^(t)的平方和再开方。
原帖由 wxk8000 于 2008-1-9 16:05 发表
plot((0:nfft/2-1)/nfft*fs,p(1:nfft/2));
((0:nfft/2-1)/nfft*fs是计算FFT后频谱各谱线对应的(正)频率值,而p(1:nfft/2)是包络谱的谱线值。 刚才看了下matlab的帮助文件,谈谈我的看法:
1、matlab中的y=hilbert(x)命令虽然字面上是hilbert变换,实际上得到的复数序列y的虚部才是真正的原信号x的hilbert变换,实部就是原信号x,所以经y=hilbert(x)得到的就是原序列的解析信号,没必要再求其虚部等等。帮助文件原文:
Description:x=hilbert(xr) returns a complex helical sequence, sometimes called the analytic signal, from a real data sequence. The analytic signal x = xr + i*xi has a real part, xr, which is the original data, and an imaginary part, xi, which contains the Hilbert transform.
2、本人对于包络谱的算法也一直弄不明白,但看到楼上几位的说法:“用w(t)=x(t)+jx^(t)作为解析信号。包络谱的求法是:a(t)=X(t)和x^(t)的平方和再开方。”,感觉还是有问题。如果只是求解析信号的实部和虚部平方和,得到的仍然只是时域的变换信号,又谈的上什么“谱”的概念呢?我的理解是:对解析信号再做FFT得到的也许才是包络谱(实际上wxk8000兄的程序里已经这样做了)。不过这样做和直接对原信号做FFT有何区别,本人也不清楚,期待高人指点。
[ 本帖最后由 liyaohua522 于 2008-3-23 13:24 编辑 ] 我是初学者,看大家都很牛的样子,有点基本的问题想讨教,边际谱和功率谱是一回事吗?还有能量谱这些都有什么区别,有点晕了? 边际谱和功率谱是不一样的,边际谱是Hilbert变换以后再求得的,功率谱一般是FFT求得 楼上的问题我也困扰很久了,希望有高人能够解答 我也想问个 问题。我在做轴承的故障诊断,运用包络解调方法,我要怎样确定带通滤波器,如中心频率··· y=hilbert(d1);
y1=abs(y);
我的理解,到这,y1求的是a(t)=X(t)和x^(t)的平方和再开方,而这时的a(t)叫包络信号,不叫包络谱,要求包络谱在对a(t)做傅里叶变换,而不是对解析信号做傅里叶变化。 我将信号做hilbert变换后求的实部和虚部的平方和再开方,但是得到的并不是信号的包络。是不是因为我没滤波呢 用w(t)=x(t)+jx^(t)作为解析信号。包络谱的求法是:a(t)=X(t)和x^(t)的平方和再开方得到信号b(t)。然后对b(t)进行fft变换。
b(t)=abs(hilbert(a(t)))
p=fft(b(t))
谢谢
学习中!:@D 看了10楼和12楼的帖子就非常明白包络谱的概念以及用Matlab怎么实现了!非常感谢!
页:
[1]