关于Hilbert谱的一些问题

圣经里的菜刀

大家好，我刚开始接触HHT，在做些练习，用HHT画出来的Hilbert谱的频率范围是0到1，这个结果肯定不对，但是不知道什么原因，也不知道该怎么修改，希望大神指点，谢谢。

shuihai707

程序贴出来，图贴出来，用事实说话。

圣经里的菜刀

shuihai707 发表于 2014-6-15 20:46
程序贴出来，图贴出来，用事实说话。

嗯嗯好的，谢谢指点，我试着把程序和图贴出来，用的褚福磊的  机械故障诊断中的现代信号处理方法  这本书里的程序

圣经里的菜刀

shuihai707 发表于 2014-6-15 20:46
程序贴出来，图贴出来，用事实说话。

function [Spectrum, freq] = HHT(Sig),  %计算HHT的主程序
if(nargin <1),
    error('At least one input is needed!');
end
Imf = Emd(Sig'); %对信号进行经验模式分解
[nImf, SigLen] = size(Imf);
[nImf, SigLen] = 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,:));
    [freq,Amp] = 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([min(t) max(t) min(freq) max(freq)]);
shading interp;
title('Hilbert-Huang spectrum');
ylabel('Frequency');
xlabel('Time');

圣经里的菜刀

shuihai707 发表于 2014-6-15 20:46
程序贴出来，图贴出来，用事实说话。

function [Imf] = 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;
[IndMax, IndMin] = ExtrPoint(c);       %提取信号的极值点：极大值点和极小值点
[Index_zer] = 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([LMaxEtd t(IndMax) RMaxEtd], [hlMaxEtd h(IndMax) hrMaxEtd], t, 'spline');
  %计算信号的下包络线
    Min_Env = interp1([LMinEtd t(IndMin) RMinEtd], [hlMinEtd h(IndMin) hrMinEtd], 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));
    [IndMax, IndMin] = ExtrPoint(h);     %提取信号的极值点
    [Index_zer] = ZeroPoint(h);          %提取信号的零点
    NExtr = length(IndMin) + length(IndMax);
    NZero = length(Index_zer);
    nLoop = nLoop+1;
end   %提取到一个本征模分量
IMFNum = IMFNum + 1;
disp([num2str(IMFNum), 'IMFs have been obtained.']);
Imf = [Imf; h];
c = c - h;
[IndMax, IndMin] = ExtrPoint(c);           %提取信号的极值点
[Index_zer] = 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 = [Imf;c];
[n m ] = 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 = [xCor max(abs(Coef))];
end
k = find(xCor > max(xCor)/8);
disp([num2str(length(k)),'IMFs have been selected.']);
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 [Pos_max,Pos_min] = 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_max]);
Pos_min = sort([Pos_min])

圣经里的菜刀

function [PZero] = 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([0 zero 0]);
    LZero = find(DZero == 1);
    RZero = find(DZero ==-1) - 1;
    IndZero = round((LZero + RZero) / 2);
    PZero = sort([PZero IndZero]);
end

shuihai707

看程序太累，论坛有现成的主流程序，法国人的EMD，自己用用试试，有问题再提问吧。

圣经里的菜刀

shuihai707 发表于 2014-6-17 18:28
看程序太累，论坛有现成的主流程序，法国人的EMD，自己用用试试，有问题再提问吧。

嗯嗯，非常感谢，我找找那个程序试试