pl.m

hibert变换程序

function pl(Y,F,A,Fs)
N=240;
close all;
s=size(Y);
pp=N:length(Y)-N;
t=pp/Fs;  %  Time
dec=round(s(1)/150);

c=['k- '; 'b- '; 'r- '; 'm  ';  'g- ';  'c- ';  'y- ' ];
if s(2)>7 s(2)=7; end

F=F*Fs;  %  Frequency, [Hz]

for k=1:s(2),
figure(1); 
subplot(s(2),1,k,'align');
plot(t,Y(pp,k),c(k,:));grid on;drawnow;
ylabel([int2str(k) '^Y']);
axis([min(t) max(t) min(Y(pp,k)) max(Y(pp,k))]);
xlabel('Time, s');

figure(2);
subplot(211);
plot(t,F(pp,k),c(k,:),'Linewidth',2);drawnow;hold on;
grid on;
axis([min(t) max(t) 0 1.2*max(max(F(pp,:)))]);
ylabel('Frequency, Hz');

subplot(212);
plot(t,A(pp,k),c(k,:),'Linewidth',2);drawnow;hold on;
grid on;
axis([min(t) max(t) 0 1.2*max(max(A(pp,:))) ]);
ylabel('Amplitude')
xlabel('Time, s');

figure(3);
subplot(211)
plot(t,Y(pp,k),c(k,:));drawnow;hold on;
grid on;
axis([min(t) max(t) min(min(Y(pp,:))) max(max(Y(pp,:)))]);
ylabel('Y');

end

figure(1); subplot(s(2),1,1,'align'); title('Components');
figure(2); subplot(211); title('Component instantaneous frequency');
subplot(212); title('Component envelope');

figure(3); subplot(211); title('Components');
subplot(212); plot(t,sum(Y(pp,:)')); drawnow; grid on;
axis([min(t) max(t) min(sum(Y(pp,:)')) max(sum(Y(pp,:)'))]);
xlabel('Time, s'); ylabel('Y'); title('Sum of components')

figure(4);
for k=1:s(2),
stem3(t(1:dec:length(pp)),F(pp(1:dec:length(pp)),k),(A(pp(1:dec:length(pp)),k)),c(k),'.'); hold on;drawnow;
end
xlabel('Time, s');ylabel('Frequency, Hz');zlabel('Amplitude')
axis([min(t) max(t) min(min(F(pp,:))) max(max(F(pp,:))) min(min(A(pp,:))) max(max(A(pp,:)))]);
title('3D plot of the instantaneous frequency and the envelope of each component')
%view(-50,50);
view(-50,70);

tilefigs([2 2],10)