hhs.m

最近在看HHT

clc;clear;
t = [0:0.001:0.2]'; 
y1 = zeros(size(t)); 
y2 = zeros(size(t)); 

fs = 100;
N = length(t);
for k = 1:N 
if t(k)>=0 & t(k)<0.1; 
    y1(k) = sin(2*pi*50*t(k))+0.5*sin(2*pi*150*t(k)); 
else t(k)>=0.1 & t(k)<0.2; 
    y2(k) = sin(2*pi*50*t(k))+0.3*sin(2*pi*250*t(k)); 
    end 
end 
y = y1 + y2; 
x = y';
tn=length(t);

t1=min(t(1),t(tn));
t2=max(t(1),t(tn));

% Perform the EMD
imf = semd(x);

ci = imf;

imf_Start = 1;
imf_Stop = size(imf,2);
imf_Last = imf_Stop+1;

for i = imf_Start:imf_Stop
    %%%%%%%%%%%%%%%%%%% Plot the Hilbert Spectrum %%%%%%%%%%%%%%%%%%%

% Perform the Hilbert transform and 
% get instantaneous attributes
[omega,amp,mag,phase,w] = fhilbert(imf(:,i),t1,t2,[]);

% Plot the Hilbert Spectrum
[m,n] = size(mag);
h     = full(sum(mag,2)/n);

% Need time,frequency and amplitude to be column vectors
T = repmat(t(:),size(imf(:,i),2),1);
O = omega(:);
A = amp(:);
    % Sum amp for identical freqs in different IMFs
stoa = sortrows([T O A],[2 1]);
T = stoa(:,1);
O = stoa(:,2);
A = stoa(:,3);
myeps = max(max(abs(T)),max(abs(O)))*eps^(1/3);
ind = [0; ((abs(diff(O)) < myeps) & (abs(diff(T)) < myeps)); 0];
if sum(ind) > 0
  disp(['Summing amplitudes at identical frequencies.']);
  fs = find(ind(1:end-1) == 0 & ind(2:end) == 1);
  fe = find(ind(1:end-1) == 1 & ind(2:end) == 0);
  for n = 1 : length(fs)
    % summing z values
    A(fe(n)) = sum(A(fs(n):fe(n)));
  end
  T = T(~ind(2:end));
  O = O(~ind(2:end));
  A = A(~ind(2:end));
end

figure;
ax1 = subplot(3,1,1);
scatter(T,A,3,'+');%plot(t,sum(imf(:,i),2));
str0 = 'IMF';
str1 = '  幅值';
str2 = [str0 num2str(i) str1];
ylabel(str2);
% ylabel('幅度');
set(gca,'XTickLabel','','xlim',[min(t) max(t)]);% Amplitude
grid on;box on;

ax2 = subplot(3,1,2);
% scatter(T,O,10,20*log10(A),'o','filled');
scatter(T,O,3,'+');
set(ax2,'ydir','normal','xlim',[min(t) max(t)]);
xlabel('时间');
str3 = '  瞬时频率';
str4 = [str0 num2str(i) str3];
ylabel(str4);
% ylabel('频率');% Frequency
%set(gca,'yscale','log');
grid on;box on

% ax4 = colorbar('horiz');axes(ax4);
% xlabel('能量 (dB)')
% grid on;box on;

% ax3 = subplot(3,1,3);plot(h,w);
% xlabel('Amp');set(gca,'YTickLabel','');
% grid on;box on

set(ax1,'Position',[0.1300 0.5822 0.6770 0.3222]);
set(ax2,'Position',[0.1300 0.2100 0.6770 0.3222]);
% set(ax3,'Position',[0.8480 0.2100 0.0770 0.3222]);
% set(ax4,'Position',[0.1300 0.1000 0.6770 0.0100]);
end