%脈動信息emd分析,斜率
%用HT分析法分析瞬時頻律
clear
clc
load 'k.txt';

%選取觀察窗
window=201:2200;
sig2=k(window,1);
%采樣頻率
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)))-1;
% for i=1:li;
%     subplot(li,1,i);
%     plot((1:length(f(1,:)))/fs,f(i,:)*fs,'LineWidth',1);
%     xlabel('time(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
%     ylabel(strimf);
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%各分量的功率譜%%%%%%%%%%%%%%%%%%%%%%%%%
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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%功率譜估計比較
% Fs=100;
% nfft=1024;
% x=sig2;
% 
% figure;
% window1=boxcar(100);
% window2=hamming(100);
% window3=blackman(100);
% noverlap=40;
% [pxx1,f1]=pwelch(x,window1,noverlap,nfft,Fs);
% [pxx2,f2]=pwelch(x,window2,noverlap,nfft,Fs);
% [pxx3,f3]=pwelch(x,window3,noverlap,nfft,Fs);
% % pxx1=10*log10(pxx1);
% % pxx2=10*log10(pxx2);
% % pxx3=10*log10(pxx3);
% subplot(311);
% plot(f1,pxx1);
% title('矩形窗');
% subplot(312);
% plot(f2,pxx2);
% title('海明窗');
% ylabel('功率');
% subplot(313);
% plot(f3,pxx3);
% xlabel('頻率(Hz)');
% title('布萊克曼窗');
% 
% figure;
% [pxx1,f1]=pmtm(x,2,nfft,Fs);
% [pxx2,f2]=pmtm(x,4,nfft,Fs);
% [pxx3,f3]=pmtm(x,10,nfft,Fs);
% % pxx1=10*log10(pxx1);
% % pxx2=10*log10(pxx2);
% % pxx3=10*log10(pxx3);
% subplot(311);
% plot(f1,pxx1);
% title('Nw=2');
% subplot(312);
% plot(f2,pxx2);
% title('Nw=4');
% ylabel('功率');
% subplot(313);
% plot(f3,pxx3);
% xlabel('頻率(Hz)');
% title('Nw=10');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%