关于Hilbert谱的一些问题
大家好，我刚开始接触HHT，在做些练习，用HHT画出来的Hilbert谱的频率范围是0到1，这个结果肯定不对，但是不知道什么原因，也不知道该怎么修改，希望大神指点，谢谢。程序贴出来，图贴出来，用事实说话。 shuihai707 发表于 2014-6-15 20:46
程序贴出来，图贴出来，用事实说话。
嗯嗯好的，谢谢指点，我试着把程序和图贴出来，用的褚福磊的机械故障诊断中的现代信号处理方法这本书里的程序 shuihai707 发表于 2014-6-15 20:46
程序贴出来，图贴出来，用事实说话。
function = HHT(Sig),%计算HHT的主程序
if(nargin <1),
error('At least one input is needed!');
end
Imf = Emd(Sig'); %对信号进行经验模式分解
= size(Imf);
= size(Imf);
t = 1: SigLen;
HHSpectrum = zeros(max(1,nImf-1), SigLen);
AmpSpectrum = zeros(max(1, nImf - 1), SigLen);
for k = 1:nImf - 1,
Imf0 = real(Imf(k,:));
= InstFreq(Imf0);
HHSpectrum(k,:) = freq';
AmpSpectrum(k,:) = Amp';
end
%组合所有本征模分量的瞬时频率和瞬时功率为Hilbert谱
MaxFreq = max(max(HHSpectrum));
MaxFreq = ceil(MaxFreq/0.5)*0.5;
NumFreq = 128;
freq = linspace(0,MaxFreq,NumFreq);
Spectrum = zeros(NumFreq,SigLen);
Temp = HHSpectrum;
Temp = min(round((Temp/MaxFreq)* NumFreq)+1, NumFreq);
for k = 1:nImf - 1,
for n = 1:SigLen,
Spectrum(Temp(k,n),n) = Spectrum(Temp(k,n),n) + AmpSpectrum(k,n);
end
end
%显示信号的Hilbert谱
clf;
mesh(t, freq, Spectrum);
colormap jet;
axis();
shading interp;
title('Hilbert-Huang spectrum');
ylabel('Frequency');
xlabel('Time'); shuihai707 发表于 2014-6-15 20:46
程序贴出来，图贴出来，用事实说话。
function = Emd(Sig);%对信号进行经验模式分解
if(nargin < 1),
error('At least one input is needed!');
end
Sig = Sig';
SigLen = length(Sig);
t = 1:SigLen;
c = Sig;
Imf = [];
nLoop = 1;
MaxLoop = 1000;
= ExtrPoint(c); %提取信号的极值点：极大值点和极小值点
= ZeroPoint(c); %提取信号的零点
NExtr = length(IndMin) + length(IndMax); %极值点个数
NZero = length(Index_zer); %零点个数
Bool = 1;
NEtd = 2;
IMFNum = 0;
while( NExtr >2 & Bool >0),
h = c;
nLoop = 1;
SD = 1;
while((SD>0.3|(abs(NZero-NExtr)>1))&(NExtr>2)&nLoop<MaxLoop)%镜像拓展序列
LMaxEtd = fliplr(IndMax(1:min(end,NEtd)));
LMinEtd = fliplr(IndMin(1:min(end,NEtd)));
hlMaxEtd = h(LMaxEtd);
hlMinEtd = h(LMinEtd);
LMaxEtd = -LMaxEtd +1;
LMinEtd = -LMinEtd +1;
RMaxEtd = fliplr(IndMax(max(end-NEtd+1,1):end));
RMinEtd = fliplr(IndMin(max(end-NEtd+1,1):end));
hrMaxEtd = h(RMaxEtd);
hrMinEtd = h(RMinEtd);
RMaxEtd = 2*SigLen - RMaxEtd;
RMinEtd = 2*SigLen - RMinEtd;
%计算信号的上包络线
Max_Env = interp1(, , t, 'spline');
%计算信号的下包络线
Min_Env = interp1(, , t, 'spline');
Mean_Env = (Max_Env +Min_Env)/2; %计算上下包络线的平均值
Preh = h;
h = h - Mean_Env;
alarm = 1.0e-08;
SD = sum (((Preh - h).^2)./(Preh.^2 + alarm));
= ExtrPoint(h); %提取信号的极值点
= ZeroPoint(h); %提取信号的零点
NExtr = length(IndMin) + length(IndMax);
NZero = length(Index_zer);
nLoop = nLoop+1;
end %提取到一个本征模分量
IMFNum = IMFNum + 1;
disp();
Imf = ;
c = c - h;
= ExtrPoint(c); %提取信号的极值点
= ZeroPoint(c); %提取信号的零点
NExtr = length(IndMin) + length(IndMax);
NZero = length(Index_zer);
Bool = 1;
if((range(c)/range(Sig)) < 2e-4),
Bool = -1;
end
end
Imf = ;
= size(Imf);
xCor = [];
%选择本征模分量
norm_Sig = (Sig - mean(Sig)) / std(Sig);
for i = 1:n-1,
norm_Imf = (Imf(i,:)-mean(Imf(i,:))) / std(Imf(i,:));
Coef = xcorr(norm_Sig,norm_Imf) / SigLen;
xCor = ;
end
k = find(xCor > max(xCor)/8);
disp();
Temp = zeros(length(k)+1, SigLen);
Temp1 = Sig;
for i = 1:length(k),
Temp(i,:) = Imf(k(i),:);
Temp1 = Temp1 -Imf(k(i),:);
end
Temp(length(k)+1,:) = Temp1; %最后一个为残差分量
Imf = Temp; shuihai707 发表于 2014-6-15 20:46
程序贴出来，图贴出来，用事实说话。
function = ExtrPoint(Sig); %提取信号的极大值和极小值点
if(nargin == 0)
error('At least one input is required');
end
SigLen = length(Sig);
DSig = diff(Sig);
DSig1 = DSig(1:end - 1);
DSig2 = DSig(2: end);
Pos_min = find(DSig1.*DSig2<0 & DSig1<0) + 1;
Pos_max = find(DSig1.*DSig2<0 &DSig1>0) + 1;
Pos_max = sort();
Pos_min = sort() function = ZeroPoint(Sig); %提取信号的过零点
if(nargin == 0)
error('At least one input is required');
end
SigLen = length(Sig);
s1 = Sig(1:SigLen-1);
s2 = Sig(2:SigLen);
PZero = find(s1.*s2 <0);
IndZero = find( Sig == 0);
if(length(IndZero) > 0),
zero = find(Sig ==0);
DZero = diff();
LZero = find(DZero == 1);
RZero = find(DZero ==-1) - 1;
IndZero = round((LZero + RZero) / 2);
PZero = sort();
end
看程序太累，论坛有现成的主流程序，法国人的EMD，自己用用试试，有问题再提问吧。 shuihai707 发表于 2014-6-17 18:28
看程序太累，论坛有现成的主流程序，法国人的EMD，自己用用试试，有问题再提问吧。
嗯嗯，非常感谢，我找找那个程序试试
页:
[1]