% This script demonstrates the use of MA-EMD on real-life central body
% temperature (CBT)
% This cript is best run under Matlab R2016a and tested on R2017a
% Code writer: Hui-Wen Yang 2019 @GICE, NTU and LIBSA, NCU
% Please cite:
% Hui-Wen Yang ; Shyh-Kang Jeng ; Hsu-Wen V. Young ; Chen Lin ; Yung-Hung Wang ; 
% Kun Hu ; Men-Tzung Lo. (2019). A Minimum Arclength Method for Removing Spikes 
% in Empirical Mode Decomposition. IEEE Access PP(99):1-1 ?January 2019
clear; close all;
load('cbt_data.mat');   % data, fs = 10 per hour


L = size(data,1);
t = (1:L)/fs;

sig1 = data(:,2);

% find extrema
cd ./MA_EMD_code    

%% detect spike
winsize = 20; thrld = 1.2;
[spikeIDMx, spikeIDMn] = detecspike(sig1,thrld,winsize,0);

 %%  EMD - minArc
allmode_minArc = MA_EMD(sig1,spikeIDMx,spikeIDMn);

z = sig1-allmode_minArc(:,2);
allmode_MA = eemd(z,0,1,6);
allmode_MA = allmode_MA';
IMF_ma = [allmode_minArc(:,2), allmode_MA];

%% EMD
allmode_EMD = eemd(sig1,0,1,7);
allmode_EMD = allmode_EMD';
%% Figure 1: raw data and detected spike
fig1 = figure;
plot(t,sig1,'Color',[0 0 0],'LineWidth',1.5); hold on;
plot(t(spikeIDMx),sig1(spikeIDMx),'o');
plot(t(spikeIDMn),sig1(spikeIDMn),'*');
legend('raw data','detected spike(max)','detected spike(min)')
title('core body temperature data')
ylabel('Temperature (^oC)')
xlabel('Time (hour)')

%% Figure 2: IMFs from EMD and MA-EMD
fig2 = figure;

ax = zeros(6,2);
for i=1:6
    % left panel: EMD
    ax(i,1) = subplot(6,2,i*2-1);
    plot(t,allmode_EMD(:,i),'Color',[0 0 0],'LineWidth',1);
    ylabel(['IMF ' num2str(i)])
    % right panel: MA-EMD
    ax(i,2) = subplot(6,2,i*2);
    plot(t,IMF_ma(:,i),'Color',[0 0 0],'LineWidth',1);
    ylabel(['IMF ' num2str(i)])
end
title(ax(1,1),'EMD')
title(ax(1,2),'MA-EMD')

