Getting Started
Examples on how to use the FECGSYN
Examples on how to use the FECGSYN
The code will produce simulation files, make sure to select the desired destination path on the Matlab root (cd savedir). An exemplary dataset can be simply generated and plotted by running the following lines:
param.fs = 1000; % data generation at 1kHz
param.n = 15*param.fs;
param.SNRfm = -3; % fetal maternal SNR
param.SNRmn = 8; % maternal to noise SNR
param.ntype = {'MA'}; % adding muscular artifact type noise
param.noise_fct = {1}; % noise at constant SNR level
debug = 5; % debug level (5 = max number of plots generated)
out = run_ecg_generator(param,debug);
FECGSYN_plotmix(out) % plots random channels of one of the created datasets
Please consider that various parameters (e.g. location of fetal heart, location of noise sources) are not explicitly defined in the previous example. Therefore, to reproduce physiological conditions, running the code twice will produce considerably different signals. In the following examples, more complex simulations are performed, where further parameters are controlled.
This example make use of the exp_datagen1.m
function to generate the dataset used in (Behar et al 2014). The code will produce simulation files, so be sure to cd dirname the desired destination path.
pathdir = pwd;
exp_datagen1(pathdir,5,false) % generate dataset (may take a few minutes)
load('fecgsyn_c1.mat')
FECGSYN_plotmix(out) % plots random channels of one of the created datasets
This example uses the asymmetric volume conductor modeling process as proposed in (Keenan et al 2018) to project the maternal and fetal source models onto each sensor. This code requires that you have installed the pre-processed anatomic models into your fecgsyn toolbox. This code converts these anatomic models into a finite element mesh with specified compartment mesh size and conductivity, then computes the lead field matrices for a pre-defined set of source and sensor positions as follows:
models = get_anatomic_models(); % check available anatomic models
model = read_anatomic_model(models.femonum_32WGA); % read selected model
model.compartments.vernix.mesh = model.variants.vernix.back_3mm; % set vernix model as 3mm in the back region
model.compartments.vernix.conductivity = 1e-5; % set vernix compartment conductivity to 1e-5 S/m
model.compartments.vernix.meshsize = 2; % set maximum vernix compartment mesh size to 2mm^3
model = generate_fe_mesh(model); % generate finite element mesh and assign to model
plot_fe_mesh(model); % plot finite element mesh
model = calculate_lead_field(model); % calculate lead field matrices and assign to model
This code may take several minutes to run depending on your system configuration and will generate some large files during this process. Please note that the selected mesh size will influence the final solution quality and computational cost (see Fig 8. in Keenan et al 2018 for further details). Once you have generated the lead field matrices and assigned them to a model, you can simulate abdominal potentials using the following code:
param.fs = 1000; % data generation at 1kHz
param.n = 15*param.fs;
param.SNRfm = -3; % fetal maternal SNR
param.SNRmn = 8; % maternal to noise SNR
param.ntype = {'MA'}; % adding muscular artifact type noise
param.noise_fct = {1}; % noise at constant SNR level
debug = 5; % debug level (5 = max number of plots generated)
out = run_ecg_generator(param,debug);
plot_asymmetric_mixture(out,model) % plots mixture using vectorcardiogram in out and lead field matrices in model
Similarly to the exp_datagen1.m
example, the FECGSYN includes three further pre-built routines to data generation. For performance purpose, these functions DO NOT include wfdb support. Nevertheless, this can be simply expanded by using the fecgsyn2wfdb.m
function (see this example). Below is a list of further examples included and a short description. For mor information, see :
FECGSYNDB_datagen.m
: This function generates a dateset similar to the one used in (Andreotti et al 2016).exp_datagen2.m
: This example performs some additional maternal ECG amplitude changes (after generation) and SNR modulation (during generation). It is useful to stress testing the ability of Template Subtraction techniques in adapting their average template to an ECG with varying amplitude. More information is available in the Documentation. exp_datagen3.m
: This example include highly oscilatting (e.g sinusoidal) changes in the noise SNR level and ECG amplitude changes. The aim of the dataset generate is to train algorithms that detect variations in the MECG/FECG signal quality, i.e. signal quality indexes (SQIs). More information is available in the Documentation.FECGSYNDB_datagen.m
will NOT be identical. For comparison purposes with (Andreotti et al 2016), make sure to download Physionet's FECGSYNDB ( available here). Reference:
Behar J., Andreotti F., Zaunseder S., Li Q., Oster J. and Clifford G D., An ECG model for simulating maternal-foetal activity mixtures on abdominal ECG recordings. Physiol Meas 35(8), pp.1537-50, 2014.
Andreotti F., Behar J., Zaunseder S.,Oster J. and Clifford G D., An Open-Source Framework for Stress-Testing Non-Invasive Foetal ECG Extraction Algorithms. Physiol Meas 37(5), pp. 627-648, 2016.
Keenan E., Karmakar C K. and Palaniswami M., The effects of asymmetric volume conductor modeling on non-invasive fetal ECG extraction. Physiol Meas 39(10), pp. 105013, 2018
This minimalistic example makes use of the template subtraction (TS) method by Cerutti et al 1986, using the FECGSYN_ts_extraction.m
function, using the FECGSYN_test.mat examplary set. As input the function expects the maternal QRS locations (mref
), abdominal signal (mixture
- rows represent different channels), a predetermined number of cycles to use on the average maternal template (NbCycles
) and sampling frequency used (fs
). See the example:
load('FECGSYN_test.mat')
NbCycles = 20; % number of cycles used
ch = 3;
residual = FECGSYN_ts_extraction(out.mqrs,out.mixture(ch,:),'TS-CERUTTI',0,...
NbCycles,'',out.param.fs);
% Plotting results
plot(out.mixture(ch,:),'Color',[0.8 0.8 0.8])
hold on
plot(residual)
plot(out.mqrs,2.*ones(1,length(out.mqrs)),'gd')
plot(out.fqrs{1},1,'xb')
legend('abdm mixture','extracted signal','MQRS','FQRS')
hold off
As output, the function returns the residual
variable, which has the same dimensions as mixture
.
Reference:
Cerutti, S., Baselli, G., Civardi, S., Ferrazzi, E., Marconi, A. M., Pagani, M., & Pardi, G. (1986). Variability analysis of fetal heart rate signals as obtained from abdominal electrocardiographic recordings. Journal of Perinatal Medicine, 14(6), 445–452.
FECGSYN_ts_extraction.m
function has several pre-built TS methods that can be used by simply substituting 'TS-CERUTTI'
for 'TS'
, 'TS-SUZANNA'
, 'TS-LP'
or 'TS-PCA'
. Make sure to see the Documentation. Aside from Template Subtraction methods (FECGSYN_ts_extraction.m
function), the FECGSYN toolbox includes Blind Source Separation (FECGSYN_bss_extraction.m
function), Kalman Filter (FECGSYN_kf_extraction.m
function) and Adaptive Filters (FECGSYN_adaptfilt_extraction.m
function). For more information, checkout Andreotti et al 2016 and the Documentation.
In this example, extraction using every method available is performed by applying the FECGSYN_main_extract.m
script. The script accepts WFDB formatted data, as well as .mat
files (following FECGSYN's structure). In order to inform the function which format is used, the bool wfdb
should be toggled to true or false. The function comes with to preset pre-processing bands (3-100 Hz or 0.5-100 Hz) which can be set with the argument narrowband
. The function may be applied directly on the FECGSYNDB, as follows:
pathdir = fullfile(which('FECGSYN_test.mat'));
pathdir = pathdir(1:end-16);
wfdb = false; % true/false, data format in WFDB?
narrowband = false; % preprocessing filtering band
ch = 3; % extract channel 3 of test data
refchs = 1; % use channel 1 as reference
debug = 1; % plot results
FECGSYN_main_extract(pathdir,narrowband,wfdb,ch,refchs,debug)
This function saves each extracted signal in an separate .mat
file under the directory 'path/data/ext/'
. In order to reduce the filenames, a coding is introduce consisting of 'rec' (for recording) and a number, followed by the extraction method used, for example rec1_tsc.mat
or rec5_PCA.mat
. An index.mat
file is also saved on the folder to allow the backwards mapping of simulated and extracted datasets.
This example demonstrates the signal quality estimation proposed in Andreotti et al 2017. This can be used to select/merge multichannel extracted FECG information or e.g. in BSS techniques component selection.
First, we modify the initial example by simulating a minute of data and modulating noise functions to a sinusoidal function and a constant function by setting the variable, as follows:
param.fs = 1000; % data generation at 1kHz
param.n = 60*param.fs;
param.SNRfm = -3; % fetal maternal SNR
param.SNRmn = 8; % maternal to noise SNR
param.ntype = {'MA'}; % adding muscular artifact type noise
param.noise_fct = {arrayfun(@(x) 1+sin(5*pi*x),linspace(0,1,param.n)) 0.5}; % variable noise
debug = 0; % debug level (skip it)
out = run_ecg_generator(param,debug);
Then extract a few channels of the produced simulation using for example the TSpca algorithm (Kanjilal et al 1997) using:
NbCycles = 20;
mixture = out.mixture([1 15 30],:);
for j = 1:size(mixture,1)
residual(j,:) = FECGSYN_ts_extraction(out.mqrs,mixture(j,:),'TS-PCA',0,...
NbCycles,'',out.param.fs);
end
This example illustrates how to produce a single signal quality estimate for a single channel segment using the bsqi algorithm (Li et al 2008), which compares how similar the result of two different QRS detectors are.
qrs1 = OSET_MaxSearch(residual(1,:),1.3/out.param.fs); % 2.3 Hz is the expected FHR parameter
qrs2 = jqrs_mod(residual(1,:),0.25,0.3,out.param.fs,''); % refractory period and threshold as parameters
sqi = bsqi(qrs1,qrs2,0.05,out.param.fs);
sprintf('SQI for segment was equal %2.2f',sqi)
Obviously, the previous example is suboptimal since we are estimating SQI for a complete 1 minute segment, during which we simulated a varying noise modulation. In the following example, the several SQI metrics used in Andreotti et al 2017 are applied. The example makes use of 5sec windows in calculating the different fetal SQI metrics for each channel. The results are saved in a table sqi
. Further, the different metrics are then merged using the pre-trained Naive Bayes Classifier into one metric as follows:
sqi = main_sqi_analysis(out.mqrs,out.mecg,mixture,residual,out.param.fs,'testrec');
load('nbmodel.mat') % loading classifier
for ch = 1:max(sqi.channel)
idx = sqi.channel == ch;
features = sqi{idx,5:end}; % obtaining features
[estimate,scores] = nbmodel.predict(features);
final_sqi(:,ch) =scores*nbmodel.ClassNames; % get continues value based on posterior probabilites
end
Where the columns of final_sqi
represent the window-based SQI values for each of the available channels. The results for one of the channels (e.g. channel 2) can be plotted as:
ch = 2;
ax(1)=subplot(3,1,1);
plot(mixture(ch,:),'Color',[0.8 0.8 0.8])
hold on
plot(residual(ch,:),'b')
plot(qrs1,residual(ch,qrs1),'xr')
legend('abdmECG','FECG','FQRS')
ax(2)=subplot(3,1,2);
sqiidx = [10 30 35 40]; % see variable "sqi" for more info
plot(table2array(sqi(sqi.channel==ch,4)),table2array(sqi(sqi.channel==ch,sqiidx)))
legend(sqi.Properties.VariableNames(sqiidx))
ax(3)=subplot(3,1,3);
plot(table2array(sqi(sqi.channel==ch,4)),final_sqi(:,ch))
legend('Final SQI')
linkaxes(ax,'x');
gqrs()
and wqrs()
detectors. For more information see FECGSYN dependencies. In this section, the multichannel Kalman filter (KF) fetal heart rate (FHR) estimation is presented. Following the previous steps shown in the previous section ( Multichannel SQI estimatione), the FHR can be estimated by applying the function runKFHR
as:
runKFHR(residual,final_sqi,out.fqrs{1},out.param.fs)
The algorithm should plot the results for window-based FHR estimation using the method proposed by Li et al 2008.
References:
Andreotti et al. 2017 Non-Invasive Fetal ECG Signal Quality Assessment for Multichannel Heart Rate Estimation. IEEE Trans. Biomed. Eng. (in press).
Kanjilal, P. P., Palit, S., & Saha, G. (1997). Fetal ECG extraction from single-channel maternal ECG using singular value decomposition. IEEE Trans. Biomed. Eng., 44(1), 51–9.
Li, Q., Mark, R. G., & Clifford, G. D. (2008). Robust heart rate estimation from multiple asynchronous noisy sources using signal quality indices and a Kalman filter. Physiol. Meas., 29(1), 15-32.
slashchar = char('/'*isunix + '\'*(~isunix)); % find out OS
old_path=which('FECGSYNDB_datagen');
[~,old_path] = strtok(old_path(end:-1:1),slashchar);
old_path = fliplr(old_path);
load('test.mat')% load test set
NbCycles = 20; % number of cycles used
ch = 3;
residual = FECGSYN_ts_extraction(out.mqrs,out.mixture(ch,:),'TS-CERUTTI',0,...
NbCycles,'',out.param.fs);
Single channel FQRS detection can be performed on extracted signals using the jqrs_mod.m
function by Behar et al 2014. The detector is a simple offline implementation of the Pan & Tompkins algorithm, adapted for FECG (i.e. the refractory period is shorter than for adult ECG). The function returns the FQRS locations, given a signal
, threshold level (thr
, same unit as signal amplitude), refractory period (refrac
, in s) and sampling frequency (fs
, in Hz). It is called as fqrs = jqrs_mod(signal,thr,refrac,fs)
. The following example can be used subsequently to the Simple TSc Extraction Example:
fqrs = jqrs_mod(residual,0.25,0.100,250,[],[],1); % performs single channel FQRS detection
The function Bxb_compare.m
, based on Physionet's bxb()
, enables the assessment on FQRS detection accuracy, in terms of F1-measure (F1
), mean absolute error (MAE
- a distance measure), positive predictive value (PPV
), sensitivity (SE
), number of true positives (TP
), false negatives (FN
) and false positives (FP
). As input the reference (out.fqrs{1}
) and test FQRS locations (fqrs
) sould be provided, along with an acceptance interval (acceptint
- in samples).
acceptint = 50*out.param.fs/1000; % 50 ms acceptance interval interval
[F1,MAE,PPV,SE,TP,FN,FP] = Bxb_compare(out.fqrs{1}, fqrs, acceptint)
The function FECGSYN_QRSmincompare.m
makes use of the Bxb_compare.m
on a minute basis, as in Andreotti et al 2016. A broad benchmark across different methods is also available, however a slightly more complex test set is required. This dataset is generated and extracted using the following functions:
slashchar = char('/'*isunix + '\'*(~isunix)); % find out OS
narrowband = false; % preprocessing filtering band
ch = 3; % extract channel 3 of test data
refchs = 1; % use channel 1 as reference
wfdb = false; % true/false, data format in WFDB?
debug = 0; % plot results
mkdir testset
cd testset
FECGSYNDBmini_datagen(pwd,0) % generates a miniture version of the FECGSYNDB set
FECGSYN_main_extract([pwd slashchar],narrowband,wfdb,ch,refchs,debug)
FECGSYN_benchFQRS.m
is then used to benchmark the different methods applied on FQRS detection as follows:
pathdir = pwd;% for original and extracted data
ch = []; % ch = optional input, case data in WFDB format
debug = 5; % enable plotting
stats=FECGSYN_benchFQRS(pathdir,ch,debug)
References:
Pan, J., & Tompkins, W. J. (1985). A Real-Time QRS Detection Algorithm. IEEE.J.BME, (3), 230–236.
Behar, J., Oster, J., & Clifford, G. D. (2014). Combining and Benchmarking Methods of Foetal ECG Extraction Without Maternal or Scalp Electrode Data. Physiological Measurement, 35(8), 1569–1589.
Analogous to the FQRS detection benchmark, morphological features can be evaluated using the FECGSYN_benchMorph.m
function. The function can be summarized into the following steps:
The function may be called by using the following code (after obtaining the variables pathdir
, ch
and debug
as previously shown):
ch = 1:8; % using channels from testset
morph=FECGSYN_benchMorph(pathdir,ch,debug)
pathdir
where extracted signals have a broader band than by the one used in FECGSYN_benchFQRS.m
. This because a greater high-pass cut frequency will inevitably distort the FECG T-wave.The script may take a while to run since it makes use of the ECGPUWAVE function, where several WFDB - Matlab transformations take place. If you are experiencing errors while running this function or no segmentation has taken place, chances are that your installation of the WFDB App for Matlab is not working properly. In nominal function, some figures demonstrating successful segmentations should appear and files with the extension *.ecgpu
should be created on the current directory, which indicate that segmenation has taken place. In case of errors check the existence of these files and recheck your WFDB App installation.
Reference:
Jane, R et al. 1997 Evaluation of An Automatic Threshold Based Detector of Waveform Limits in Holter ECG with the QT Database. Computers in Cardiology.
WFDB is an open-source and widely used data format from Physionet ( www.physionet.org). It is also the format adopted by the FECGSYNDB ( available here).
The following subsections aim at demonstrating how to convert from FECGSYN Matlab/Octave format to WFDB format, and vice-versa. In order to use these functions, the WFDB App should be installed, as described in Installing the FECGSYN.
This is done using the fecgsyn2wfdb.m
function. The function can be used to convert multiple .mat
saved files using FECGSYN internal structure, such as:
fecgsyn2wfdb('/path_to_mats/')
% outstruct = structure containing dataset (intern from FECGSYN)
% recName = string containing name which will be used to save WFDB dataset
% dest_path = string containing destination folder for WFDB dataset
fecgsyn2wfdb(dest_path,outstruct,recName)
outstruct
before exporting it to WFDB format. This can be done using the clean_compress.m
function, see the Documentation. exp_datagen1.m
function for a simple use of the fecgsyn2wfdb.m
, see this example. Similarly, the wfdb2fecgsyn.m
function is able to convert WFDB datasets into the .mat
internal structure used by FECGSYN FECGSYN_main_extract.m
script. This function is particularly useful when converting the recordings from the FECGSYNDB to Matlab/Octave. An examplary code follows:
% pathdir = string containing origin folder for WFDB dataset
wfdb2fecgsyn(pathdir)
pathdir
must be an unique identifier for a WFDB dataset, which excludes the extension of the files. For instance, for the FECGSYNDB 'sub01_snr00dB_l1_c0'
represents an unique identifier for a dataset consisting of multiple files: _mecg
,_fecg1
,_noise1
with .dat
, .hea
, .mqrs
and .fqrs
extensions.main_extract_data.m
function for a simple use of the wfdb2fecgsyn.m
, see this example.