Segmentation label/Custom

andbiz commented 7 years ago

An additional window (end part of the signal, label=2) is expected.

# import libraries
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

# for windowed plots
get_ipython().magic(u'matplotlib qt')

# import all pyphysio classes and methods
import pyphysio as ph 

# In[2]:

# import data and creating a signal

ecg_data = ph.TestData.ecg()

fsamp = 2048
ecg = ph.EvenlySignal(values = ecg_data, sampling_freq = fsamp, signal_nature = 'ecg')

# In[3]:

ecg.plot()

# ** Step 1: Filtering and preprocessing **

# In[4]:

# (optional) IIR filtering : remove high frequency noise
ecg = ph.IIRFilter(fp=45, fs = 50, ftype='ellip')(ecg)

# In[5]:

# normalization : normalize data
ecg = ph.Normalize(norm_method='standard')(ecg)

# In[6]:

# resampling : increase the sampling frequency by cubic interpolation
ecg = ecg.resample(fout=4096, kind='cubic')
fsamp = 4096

# In[7]:

ecg.plot()

# ** Step 2: Information Extraction **
# 
# The information we want to extract from the ECG signal is the position of the heartbeats and the Inter Beat Interval signal.

# In[8]:

ibi = ph.BeatFromECG()(ecg)

# In[9]:

# check results so far
ax1 = plt.subplot(211)
ecg.plot()
plt.vlines(ibi.get_times(), np.min(ecg), np.max(ecg))

plt.subplot(212, sharex = ax1)
ibi.plot('o-')
plt.vlines(ibi.get_times(), np.min(ibi), np.max(ibi))
plt.show()

# ** Step 3: Physiological Indicators **

# In[10]:

# define a list of indicators we want to compute
hrv_indicators = [ph.Mean(name='RRmean'), ph.StDev(name='RRstd'), ph.RMSSD(name='rmsSD'), 
              ph.PowerInBand(name = 'HF', interp_freq=4, freq_max=0.4, freq_min=0.15, method = 'ar'),
              ph.PowerInBand(name = 'LF', interp_freq=4, freq_max=0.15, freq_min=0.04, method = 'ar')
             ]

# In[11]:

# create fake label
label = np.zeros(1200)

label[300:600] = 1

label[900:1200] = 2

label = ph.EvenlySignal(label, sampling_freq = 10, signal_nature = 'label')

# In[25]:

# check label
ax1 = plt.subplot(211)
ibi.plot('.-')

plt.subplot(212, sharex = ax1)
label.plot('.-')
plt.show()

# In[16]:

#fixed length windowing
fixed_length = ph.FixedSegments(step = 5, width = 20, labels = label)

indicators, col_names = ph.fmap(fixed_length(ibi), hrv_indicators)

# In[19]:

# extract column with the labels for each window
label_w = indicators[:, np.where(col_names == 'label')[0]]

# extract column with the RRmean values computed from each window
rrmean_w = indicators[:, np.where(col_names == 'RRmean')[0]]

# create a box and whisker plot
plt.boxplot([rrmean_w[label_w==1], rrmean_w[label_w==2]],  
            labels=['image1', 'image2'])
plt.show()

# In[21]:

#label based windowing
label_based = ph.LabelSegments(labels = label)

indicators, col_names = ph.fmap(label_based(ibi), hrv_indicators)

# In[23]:

indicators[:,0:2]

# In[32]:

custom_based = ph.CustomSegments(begins = [0, 30, 60, 90], ends = [30, 60, 90, label.get_duration()], labels=label, drop_shorter=False, drop_mixed=False)

indicators, col_names = ph.fmap(custom_based, hrv_indicators, ibi)

indicators[:,0:2]

# In[29]:

len(label)

# In[ ]:

# ### 2.2 EDA processing pipeline

# ** Step 0: Import data **

# In[ ]:

# import data and creating a signal

eda_data = ph.TestData.eda()

fsamp = 2048
eda = ph.EvenlySignal(values = eda_data, sampling_freq = fsamp, signal_nature = 'eda')

eda.plot()

# ** Step 1: Filtering and preprocessing **

# In[ ]:

# resampling : decrease the sampling frequency by cubic interpolation
eda = eda.resample(fout=8, kind='cubic')

# In[ ]:

# IIR filtering : remove high frequency noise
eda = ph.IIRFilter(fp=0.8, fs = 1.1, ftype='ellip')(eda)

# In[ ]:

eda.plot()

# ** Step 2: Information Extraction **
# 
# The information we want to extract from the EDA signal is the phasic component associated to the sympathetic activity.

# In[ ]:

# estimate the driver function
driver = ph.DriverEstim(delta=0.02)(eda)

# In[ ]:

# compute the phasic component
phasic, tonic, _ = ph.PhasicEstim(delta=0.02)(driver)
phasic.plot()

# In[ ]:

# check results so far
ax1 = plt.subplot(211)
eda.plot()

plt.subplot(212, sharex = ax1)
driver.plot()
phasic.plot()
plt.grid()
plt.show()

# ** Step 3: Physiological Indicators **

# In[ ]:

# define a list of indicators we want to compute
indicators = [ph.Mean(), ph.StDev(), ph.AUC(), 
              ph.PeaksMean(delta=0.02),
              ph.DurationMean(delta=0.02)
             ]

# In[ ]:

# define the windowing method

andbiz commented 7 years ago

Add duration= default = time last sample - start time

Change segmentation (drop shorter if duration < width) stop segmentation if start > end time.