这是我的hht谱图程序，请大家帮忙看一下对不对
本来觉得挺对的，继续往下做的时候怎么也不对，就回头来看了 function =Book_HHT(Sig) %计算HHT的主程序if(nargin<1), error('At least one input is needed');endImf=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); %用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(); shading interp; title('Hilbert-Huang spectrum'); ylabel('Frequency');xlabel('Time');%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; function=ExtrPoint(Sig) %提取信号的极大值和极小值点 if(nargin==0) error('At least one input is required');endSigLen=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');endSigLen=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(); LZero=find(DZero==1); RZero=find(DZero==-1)-1; IndZero=round((LZero+RZero)/2); PZero=sort();endfunction=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=;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); 上面是用到的程序，下面是主程序，需要用到emd.mSampFreq=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);
=Book_HHT(Sig); %谱图
你好！不知道你的HHT程序处理好了没？能否发给我一份研究学习一下，我现在编的程序在调用toimage.m后，幅值和瞬时频率个数不相等了，不知道咋办，我的qq：894881296 看起来好复杂的样子，mark一下，学习学习 。 很好的参考，研究下。
页:
[1]