这是我的hht谱图程序，请大家帮忙看一下对不对

hht123

本来觉得挺对的，继续往下做的时候怎么也不对，就回头来看了

hht123

function [Spectrum,freq]=Book_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);                      %用hilbert变换估计每个本征模分量的瞬时频率和功率

     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');

%  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;

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(~isempty(IndZero)),

    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

function[freq,Amp]=InstFreq(Sig)                                            %计算信号的瞬时频率和功率

    %  freq:信号的瞬时频率

    %  Amp :信号的瞬时功率

if(nargin==0),

    error('At least one input is required');

end;

SigLen=length(Sig);

x=hilbert(real(Sig));                                                       %计算信号的Hilbert变换

Amp=abs(x);                                                                 %计算信号的瞬时功率

freq=(angle(-x(2:SigLen).*conj(x(1:SigLen-1)))+pi)/(2*pi);                  %计算信号的瞬时功率

freq=[freq(1) freq];

std_freq=std(freq);

k=3;

threshold=k*std_freq+mean(freq);

warm= freq>threshold;

freq(warm)=mean(freq);

threshold=-k*std_freq+mean(freq);

warm= find(freq<threshold);

freq(warm)=mean(freq);     

hht123

上面是用到的程序，下面是主程序，需要用到emd.m
SampFreq=500000;
t=0:1/SampFreq:0.014;
Sig1=(t>=0&t<=0.014).*(1+sin(2*pi*15000*t)).*cos(2*pi*60000*t+sin(2*pi*15000*t));
Sig2=(t>=0&t<=0.028).*(1+sin(2*pi*20000*t)).*cos(2*pi*150000*t+sin(2*pi*20000*t));
Sig3=(t>=0.1128&t<=0.0084).*cos(2*pi*150000*t.*(1+sin(2*pi*20000*t))).*cos(2*pi*150000*t+sin(2*pi*20000*t));
x=Sig1+Sig2+Sig3;
Sig=awgn(x,0.1);                          %信噪比为0.5dB

figure(1);
[Spectrum,freq]=Book_HHT(Sig);               %谱图

feiyu12345

你好！不知道你的HHT程序处理好了没？能否发给我一份研究学习一下，我现在编的程序在调用toimage.m后，幅值和瞬时频率个数不相等了，不知道咋办，我的qq：894881296

leaukai

看起来好复杂的样子，mark一下，学习学习 。

cufflink

很好的参考，研究下。