HHT边际谱的定义以及理解
(1)我对HHT边际谱的定义理解为就是将时频分布对时间求积分,具体地说就是将频率相同的时刻所对应的幅值进行累加。这样理解正确吗?(2) 如果我的理解是正确的,那么写程序时关键的一步就是判断在各种IMF下哪些时刻的频率是相同的。在实际计算中,计算机的精度导致了各种IMF下各时刻的频率值几乎不可能是完全一致的,那么我要用什么样的算法来做边际谱呢?哪位高手能提供一下具体的算法或者有现成的code可以下载?
[ 本帖最后由 yejet 于 2006-9-17 12:57 编辑 ] Marginal spectrum是对Hilbert spectrum即H(A,w,t)进行时间积分得到的。Hilbert spectrum是指H(A,w,t)时、频、幅三维图。 Fourier变换是将任何信号分解为正弦信号的加权和,而每一个正弦信号对应着一个固定的频率(Fourier频率)和固定的幅值,因此,用Fourier 变换分析频率不随时间变化的平稳信号是十分有效的。但对于频率随时间变化的非平稳信号,Fourier 变换就无能为力了。
HHT是历史上首次对Fourier变换的基本信号和频率定义作的创造性的改进。他们不再认为组成信号的基本信号是正弦信号,而是一种称为固有模态函数的信号,也就是满足以下两个条件的信号: (1) 整个信号中,零点数与极点数相等或至多相差1 ; (2) 信号上任意一点,由局部极大值点确定的包络线和由局部极小值点确定的包络线的均值均为零,即信号关于时间轴局部对称。
无论Hilbert谱中的频率还是边际谱中的频率(即瞬时频率) ,其意义都与Fourier分析中的频率(即Fourier 频率) 完全不同,但在Fourier分析中,某一频率处能量的存在,代表一个正弦或余弦波在整个时间轴上的存在,而边际谱h中某一频率处能量的存在仅代表在整个时间轴上可能有这样一个频率的振动波在局部出现过,h越大,代表该频率出现的可能性越大。
[ 本帖最后由 yejet 于 2006-9-17 12:57 编辑 ] 目前最好的贴
感谢 meil32的回贴
我正急需emd分解的matlab的程序呢,meil32的回贴借了我的燃眉之急。 HHT问题很多,下面是别人提出来的(1)因为分解的不稳定,分解所得的固有模态函数并非是系统的固有(一定)模态。
(2)且正因为分解的不稳定,分解结果一般也就很难用确切的物理过程来解释了。
(3)此外,再对一般情况下不稳定分解所得的本来就无明确物理意义的分量进行诸如波内调频(intra-wave frequency modulation)的细致分析,无疑亦是在做多此一举的无用功了。
[ 本帖最后由 yejet 于 2006-9-17 12:57 编辑 ] to:meil32,你的源程序里好像缺一个函数吧,extr(m),缺了这个不能正常运行啊,麻烦再传一下了,谢谢!
to mail32
对你的程序里却两个函数,是不能正常运行的,你要是有的话请传上来,或发到邮箱wugeren@126.com[ 本帖最后由 yejet 于 2006-9-17 12:57 编辑 ]
HHT著作已经出版,
HHT著作已经出版,那位高人去搞个电子版的 !http://www.worldscibooks.com/mathematics/5862.html
http://www.amazon.com/gp/product ... s&v=glance&n=283155
[ 本帖最后由 yejet 于 2006-9-17 12:57 编辑 ] 这两本书那位高手能介绍买的地方,谢谢了!!! 非常感谢你的程序
[ 本帖最后由 zhlong 于 2007-10-15 22:29 编辑 ] 谢谢,我正在求这个程序,谢谢!
回复:(qlisikefed)[转帖]HHT边际谱
其中缺少的extr试试下面的代码%extr
function =extr(x)
=uigetfile('*.txt','输入信号');
watchon;
FILENAME=;
x=load(FILENAME);
numzer=1;
nummin=1;
nummax=1;
for g=2:(length(x)-1)
if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
indzer(numzer)=g-1;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
indzer(numzer)=g;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
indzer(numzer)=(2*g-1)/2;
numzer=numzer+1;
end
if (x(g-1)>x(g))&(x(g)>x(g+1))
continue;
elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
nummin=nummin+1;
elseif (x(g-1) continue;
elseif (x(g-1)x(g+1))
indmax(nummax)=g;
nummax=nummax+1;
end
end
%subplot(2,3,1)
plot(indzer,x(indzer),'r-+')
title('indzer')
pause
% subplot(2,3,2)
plot(indmin,x(indmin),'b-o')
title('indmin')
pause
%subplot(2,3,3)
plot(indmax,x(indmax),'g-*')
title('indmax')
pause
%subplot(2,3,4)
plot(1:length(x),x,'k-s')
title('x')
pause
envmax = interp1(indmax,x(indmax),1:length(x),'spline');
envmin = interp1(indmin,x(indmin),1:length(x),'spline');
%subplot(2,3,5)
plot(1:length(x),envmin,'y-d')
title('envmin')
pause
%subplot(2,3,6)
plot(1:length(x),envmax,'c-h')
title('envmax')
pause
plot(1:length(x),mean((envmax+envmin)/2),'m-p')
title('mean((envmax+envmin)/2)')%extr
function =extr(x)
=uigetfile('*.txt','输入信号');
watchon;
FILENAME=;
x=load(FILENAME);
numzer=1;
nummin=1;
nummax=1;
for g=2:(length(x)-1)
if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
indzer(numzer)=g-1;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
indzer(numzer)=g;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
indzer(numzer)=(2*g-1)/2;
numzer=numzer+1;
end
if (x(g-1)>x(g))&(x(g)>x(g+1))
continue;
elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
nummin=nummin+1;
elseif (x(g-1) continue;
elseif (x(g-1)x(g+1))
indmax(nummax)=g;
nummax=nummax+1;
end
end
%subplot(2,3,1)
plot(indzer,x(indzer),'r-+')
title('indzer')
pause
% subplot(2,3,2)
plot(indmin,x(indmin),'b-o')
title('indmin')
pause
%subplot(2,3,3)
plot(indmax,x(indmax),'g-*')
title('indmax')
pause
%subplot(2,3,4)
plot(1:length(x),x,'k-s')
title('x')
pause
envmax = interp1(indmax,x(indmax),1:length(x),'spline');
envmin = interp1(indmin,x(indmin),1:length(x),'spline');
%subplot(2,3,5)
plot(1:length(x),envmin,'y-d')
title('envmin')
pause
%subplot(2,3,6)
plot(1:length(x),envmax,'c-h')
title('envmax')
pause
plot(1:length(x),mean((envmax+envmin)/2),'m-p')
title('mean((envmax+envmin)/2)')%extr
function =extr(x)
=uigetfile('*.txt','输入信号');
watchon;
FILENAME=;
x=load(FILENAME);
numzer=1;
nummin=1;
nummax=1;
for g=2:(length(x)-1)
if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
indzer(numzer)=g-1;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
indzer(numzer)=g;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
indzer(numzer)=(2*g-1)/2;
numzer=numzer+1;
end
if (x(g-1)>x(g))&(x(g)>x(g+1))
continue;
elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
nummin=nummin+1;
elseif (x(g-1) continue;
elseif (x(g-1)x(g+1))
indmax(nummax)=g;
nummax=nummax+1;
end
end
%subplot(2,3,1)
plot(indzer,x(indzer),'r-+')
title('indzer')
pause
% subplot(2,3,2)
plot(indmin,x(indmin),'b-o')
title('indmin')
pause
%subplot(2,3,3)
plot(indmax,x(indmax),'g-*')
title('indmax')
pause
%subplot(2,3,4)
plot(1:length(x),x,'k-s')
title('x')
pause
envmax = interp1(indmax,x(indmax),1:length(x),'spline');
envmin = interp1(indmin,x(indmin),1:length(x),'spline');
%subplot(2,3,5)
plot(1:length(x),envmin,'y-d')
title('envmin')
pause
%subplot(2,3,6)
plot(1:length(x),envmax,'c-h')
title('envmax')
pause
plot(1:length(x),mean((envmax+envmin)/2),'m-p')
title('mean((envmax+envmin)/2)')%extr
function =extr(x)
=uigetfile('*.txt','输入信号');
watchon;
FILENAME=;
x=load(FILENAME);
numzer=1;
nummin=1;
nummax=1;
for g=2:(length(x)-1)
if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
indzer(numzer)=g-1;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
indzer(numzer)=g;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
indzer(numzer)=(2*g-1)/2;
numzer=numzer+1;
end
if (x(g-1)>x(g))&(x(g)>x(g+1))
continue;
elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
nummin=nummin+1;
elseif (x(g-1) continue;
elseif (x(g-1)x(g+1))
indmax(nummax)=g;
nummax=nummax+1;
end
end
%subplot(2,3,1)
plot(indzer,x(indzer),'r-+')
title('indzer')
pause
% subplot(2,3,2)
plot(indmin,x(indmin),'b-o')
title('indmin')
pause
%subplot(2,3,3)
plot(indmax,x(indmax),'g-*')
title('indmax')
pause
%subplot(2,3,4)
plot(1:length(x),x,'k-s')
title('x')
pause
envmax = interp1(indmax,x(indmax),1:length(x),'spline');
envmin = interp1(indmin,x(indmin),1:length(x),'spline');
%subplot(2,3,5)
plot(1:length(x),envmin,'y-d')
title('envmin')
pause
%subplot(2,3,6)
plot(1:length(x),envmax,'c-h')
title('envmax')
pause
plot(1:length(x),mean((envmax+envmin)/2),'m-p')
title('mean((envmax+envmin)/2)')%extr
function =extr(x)
=uigetfile('*.txt','输入信号');
watchon;
FILENAME=;
x=load(FILENAME);
numzer=1;
nummin=1;
nummax=1;
for g=2:(length(x)-1)
if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
indzer(numzer)=g-1;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
indzer(numzer)=g;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
indzer(numzer)=(2*g-1)/2;
numzer=numzer+1;
end
if (x(g-1)>x(g))&(x(g)>x(g+1))
continue;
elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
nummin=nummin+1;
elseif (x(g-1) continue;
elseif (x(g-1)x(g+1))
indmax(nummax)=g;
nummax=nummax+1;
end
end
%subplot(2,3,1)
plot(indzer,x(indzer),'r-+')
title('indzer')
pause
% subplot(2,3,2)
plot(indmin,x(indmin),'b-o')
title('indmin')
pause
%subplot(2,3,3)
plot(indmax,x(indmax),'g-*')
title('indmax')
pause
%subplot(2,3,4)
plot(1:length(x),x,'k-s')
title('x')
pause
envmax = interp1(indmax,x(indmax),1:length(x),'spline');
envmin = interp1(indmin,x(indmin),1:length(x),'spline');
%subplot(2,3,5)
plot(1:length(x),envmin,'y-d')
title('envmin')
pause
%subplot(2,3,6)
plot(1:length(x),envmax,'c-h')
title('envmax')
pause
plot(1:length(x),mean((envmax+envmin)/2),'m-p')
title('mean((envmax+envmin)/2)')%extr
function =extr(x)
=uigetfile('*.txt','输入信号');
watchon;
FILENAME=;
x=load(FILENAME);
numzer=1;
nummin=1;
nummax=1;
for g=2:(length(x)-1)
if (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)==0)
indzer(numzer)=g-1;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g)==0)
indzer(numzer)=g;
numzer=numzer+1;
elseif (x(g-1)*x(g)<0)|(x(g-1)*x(g)==0)&(x(g-1)~=0)&(x(g)~=0)
indzer(numzer)=(2*g-1)/2;
numzer=numzer+1;
end
if (x(g-1)>x(g))&(x(g)>x(g+1))
continue;
elseif (x(g-1)>x(g))&(x(g) indmin(nummin)=g;
nummin=nummin+1;
elseif (x(g-1) continue;
elseif (x(g-1)x(g+1))
indmax(nummax)=g;
nummax=nummax+1;
end
end
%subplot(2,3,1)
plot(indzer,x(indzer),'r-+')
title('indzer')
pause
% subplot(2,3,2)
plot(indmin,x(indmin),'b-o')
title('indmin')
pause
%subplot(2,3,3)
plot(indmax,x(indmax),'g-*')
title('indmax')
pause
%subplot(2,3,4)
plot(1:length(x),x,'k-s')
title('x')
pause
envmax = interp1(indmax,x(indmax),1:length(x),'spline');
envmin = interp1(indmin,x(indmin),1:length(x),'spline');
%subplot(2,3,5)
plot(1:length(x),envmin,'y-d')
title('envmin')
pause
%subplot(2,3,6)
plot(1:length(x),envmax,'c-h')
title('envmax')
pause
plot(1:length(x),mean((envmax+envmin)/2),'m-p')
title('mean((envmax+envmin)/2)')
[ 本帖最后由 yejet 于 2006-9-17 12:58 编辑 ]
我也正准备研究一下
本人东南大学无线电系大三学生,正准备做个SRTP,各位大侠多多指教.MEIL32的程序我收下了,缺的两个函数有时间发到syyu514@126.com.谢谢各位。[ 本帖最后由 yejet 于 2006-9-17 12:58 编辑 ] 我刚刚下载
多谢了!