maihht.m

hilbert-huang 变换分析一段脉搏波 运行maihht即可

function maihht()
%pp是件名，window是选取的观察窗
%脉动信息emd分析，包含振幅、斜率和两者乘积
%用HT分析法分析瞬时频律
clc
vw = get(gca,'View' );
sig2=load ('pp.txt');

%选取观察窗
window=101:2100;
sig2=sig2(window,1);
sig2=sig2/4096;
%采样频率
fs=100;
N=length(sig2);
T=1:N;

%EMD分解IMFs
imf_a=emd(sig2,1:length(sig2),0.2);
imf_a=imf_a';

%显示EMD结果%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%A
figure(1);
subplot(9,1,1);
plot(T/fs,sig2,'LineWidth',1);
title('各imf分量');
ylabel('A(kPa)');
li=length((imf_a(1,:)));

for i=1:li;
    subplot(li+1,1,i+1);
    plot(T/fs,imf_a(1:length(T),i),'LineWidth',1);
    xlabel('时间(s)');
    if i==1
        strimf='c1';
    end
    if i==2
        strimf='c2';
    end
    if i==3
        strimf='c3';
    end
    if i==4
        strimf='c4';
    end
    if i==5
        strimf='c5';
    end
    if i==6
        strimf='c6';
    end
    if i==7
        strimf='c7';
    end
    if i==8
        strimf='c8';
    end
    ylabel(strimf);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %估计各分量的瞬时频率并作图
DEFSPL=200;
[A,f,tt] = hhspectrum(imf_a');
% [im,tt] = toimage(A,f);
% disp_hhs(im);
% psd=hhmspectrum(A,f);
% disp_hhms(psd,DEFSPL);
psd=hhmspectrum_new(A,f,200);
disp_hhms_new(psd,DEFSPL);

figure
li=length((f(:,1)));
for i=1:li;
    subplot(li,1,i);
    psd=hhmspectrum1(A(i,:),f(i,:),200);
    interval=0.5*fs/DEFSPL;
    xvalue(1)=interval;
    yvalue=fliplr(psd);
        for j=2:200
            xvalue(j)=xvalue(j-1)+interval;
        end
    plot(xvalue,yvalue);
    xlabel('time(s)');
    if i==1
        strimf='c1';
        title('各IMF分量的频谱');
        end
    if i==2
        strimf='c2';
    end
    if i==3
        strimf='c3';
    end
    if i==4
        strimf='c4';
    end
    if i==5
        strimf='c5';
    end
    if i==6
        strimf='c6';
    end
    if i==7
        strimf='c7';
    end
    ylabel(strimf);
end