回复 3 # syxqq123 的帖子
另外,使用诸福磊等编写《机械故障诊断中的现代信号处理方法》中关于Hilbert谱matlab程序(HHT、InstFreq函数)做出的Hilbert谱图如下图,感觉效果要好一些(相对使用hhspectrum、toimage、disp_hhs函数来生成的Hilbert谱)。
后面是《机械故障诊断中的现代信号处理方法》中关于Hilbert谱matlab程序,现在正在对比这两中方法,看看到底是算法本身的差异,还是作图的差异。望老师也帮忙看看,多谢!
HHT、InstFreq函数做出的HIlbert谱图:
HHT函数:
function [Spectrum,freq]=HHT(Imf)
if (nargin<1)
error('至少需要一个输入!');
end
[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
InstFreq函数:
function [freq,Amp]=InstFreq(Sig);
if (nargin==0)
error('At least one input is required');
end
SigLen=length(Sig);
x=hilbert(real(Sig));
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=find(freq>threshold);
freq(warm)=mean(freq);
threshold=-k*std_freq+mean(freq);
warm=find(freq<threshold);
freq(warm)=mean(freq);
Hilbert谱的绘制:
figure('numbertitle','off','name','Hilbert幅值谱');
mesh(t,freq,Spectrum);
colormap jet;
axis([min(t),max(t),min(freq),max(freq)]);
set(gca,'YLim',[min(freq),max(freq)]);
shading interp;
title('Hilbert-Huang spectrum');
ylabel('Frequency/Fs');
xlabel('Time');
|