Test_emd_IF.m

% signal -> x
%begin
fs=600;
Ts=1/fs;
t=1:length(x);
t=(t-1).*Ts;
%figure('name','original signal');plot(t,x)
%IMF=eemd_test(x); %try to use eemd to calculate IOF
IMF=emd(x);
% original signal
figure('name','original signal')
plot(t,x)
xlabel('t (s)')
ylabel('Amplitude')
%IMF components
figure('name','IMF')
for n=1:size(IMF,1)
    subplot(size(IMF,1),1,n)
    plot(t,IMF(n,:))
    xlabel('t (s)')
    ylabel('Amplitude')
    %hold on
    %plot(t,IMF(n,:),'rx')
end
%fft for each IMF components
figure('name','fft')
maxf_fft=[];
for n=2:size(IMF,1)-1
    subplot(size(IMF,1)-2,1,n-1)
    [frequency,fft_result]=fft_plot(IMF(n,:),fs);
    plot(frequency,abs(fft_result(1:length(frequency))))
    xlabel('f (Hz)')
    ylabel('Amplitude')
    axis([0 100 0 max(abs(fft_result(1:length(frequency))))])
    maxf_fft_index=find(abs(fft_result(1:length(frequency)))==max(abs(fft_result(1:length(frequency)))));
    maxf_fftvalue=frequency(maxf_fft_index);
    maxf_fft(n-1)=maxf_fftvalue(1);
end
%IF not using hilbert transform
%IF property
figure('name','IF')
xlabel('frequency')
ylabel('probility (%)')
meanf=[];
maxf_emd=[];
for n=2:size(IMF,1)-1
    subplot(size(IMF,1)-2,1,n-1)
    [f,tf,f_GZC]=Instaneous_Frequency_Test(IMF(n,:),600);
    meanf(n-1)=f_GZC;
    f_range=0:1:100;%frequency range
    f_pro=hist(f,f_range);
    plot(f_range,f_pro./length(f)*100);
    xlabel('frequency')
    ylabel(strcat('IMF',num2str(n),'(%)'))
    maxf_emd_pro=max(f_pro);
    maxf_emd_pro_index=find(f_pro==maxf_emd_pro);
    maxf_emdvalue=f_range(maxf_emd_pro_index);
    maxf_emd(n-1)=maxf_emdvalue(1);
end
%IF in time domain
figure('name','IF value')
for h=2:size(IMF,1)-1
    subplot(size(IMF,1)-2,1,h-1)
    xlabel('t (s)')
    ylabel('f (Hz)')
    [f,tf,f_GZC]=Instaneous_Frequency_Test(IMF(h,:),600);
    hold on
    for n=1:length(f)
        plot([tf(n,1),tf(n,2)],[f(n),f(n)],'k')
        if n~=length(f)
            plot([tf(n,2),tf(n,2)],[f(n+1),f(n)],'k')
        end
    end
end
%IF using hilbert transform
%IF property
figure('name','IF property hilbert')
for n=2:size(IMF,1)-1
    subplot(size(IMF,1)-2,1,n-1)
    xlabel('f (Hz)')
    ylabel('Number')
     [ XPhase,f,Imag,t,f_OGZ,maxf_emdvalue ]=IF_hilbert(IMF(n,:),fs);
     maxf_emd2(n-1)=maxf_emdvalue;
     f_range=0:1:100;%frequency range
    f_pro=hist(f,f_range);
    plot(f_range,f_pro./length(f)*100);
    xlabel('frequency')
    ylabel(strcat('IMF',num2str(n),'(%)'))
end
%IF in time domain
figure('name','IF value hilbert')
for h=2:size(IMF,1)-1
    subplot(size(IMF,1)-2,1,h-1)
    xlabel('t (s)')
    ylabel('f (Hz)')
    [ XPhase,f,Imag,t,f_OGZ,maxf_emdvalue ]=IF_hilbert(IMF(h,:),fs);
    plot(t,f);
end
disp(strcat('frequency in fft from IMF2 to IMF',num2str(size(IMF,1)-1),':'))
disp(maxf_fft)
disp(strcat('frequency in emd not uing heilbert from IMF2 to IMF',num2str(size(IMF,1)-1),':'))
disp(maxf_emd)
disp(strcat('frequency in emd using hilbert from IMF2 to IMF',num2str(size(IMF,1)-1),':'))
disp(maxf_emd2)