hhtspe.m

HHT谱计算程序

function [S,freq] = HHTspe(imf,N,fs);
%%N--------分辨率
%fs-------采样频率
%s--------时间-频率-幅值矩阵
if(nargin<3)
    fs=1;
end
L = size(imf,1); 
S = []; 
     
clear x z m p freq
   
x = imf';  
z = hilbtmf(x); % 用改进的Hilbert变换，效果好一些。
m = abs(z);
p = angle(z); 

for i = 1
      freq(:,i) = instfreq(z(:,i))*fs; %乘以了采样频率
                  ceilfreq(:,i) = ceil(freq(:,i)*N);
            for j = 1:length(x)-2
            S(ceilfreq(j,i),j+1) = m(j+1,i);
         end
     end
%plot S
figure;
t=(1:length(x))/fs; 
f=(1:size(S,1))/N; %
imagesc(t,f,S); 
colorbar;
set(gca,'YDir','normal');
xlabel('Time t/s');
ylabel(' Frequency/HZ');

return