Parametric modelling of cardiac system multiple measurement signals: an open-source computer framework for performance evaluation of ECG, PCG and ABP event detectors

Abstract

The major focus of this study is to present a performance accuracy assessment framework based on mathematical modelling of cardiac system multiple measurement signals. Three mathematical algebraic subroutines with simple structural functions for synthetic generation of the synchronously triggered electrocardiogram (ECG), phonocardiogram (PCG) and arterial blood pressure (ABP) signals are described. In the case of ECG signals, normal and abnormal PQRST cycles in complicated conditions such as fascicular ventricular tachycardia, rate dependent conduction block and acute Q-wave infarctions of inferior and anterolateral walls can be simulated. Also, continuous ABP waveform with corresponding individual events such as systolic, diastolic and dicrotic pressures with normal or abnormal morphologies can be generated by another part of the model. In addition, the mathematical synthetic PCG framework is able to generate the S4–S1–S2–S3 cycles in normal and in cardiac disorder conditions such as stenosis, insufficiency, regurgitation and gallop. In the PCG model, the amplitude and frequency content (5–700 Hz) of each sound and variation patterns can be specified. The three proposed models were implemented to generate artificial signals with varies abnormality types and signal-to-noise ratios (SNR), for quantitative detection-delineation performance assessment of several ECG, PCG and ABP individual event detectors designed based on the Hilbert transform, discrete wavelet transform, geometric features such as area curve length (ACLM), the multiple higher order moments (MHOM) metric, and the principal components analysed geometric index (PCAGI). For each method the detection-delineation operating characteristics were obtained automatically in terms of sensitivity, positive predictivity and delineation (segmentation) error rms and checked by the cardiologist. The Matlab m-file script of the synthetic ECG, ABP and PCG signal generators are available in the Appendix.

Keywords::

• Electrocardiogram
• Phonocardiogram
• Arterial blood pressure
• Macroscopic modelling
• Event detector
• Performance evaluation
• Area–curve length
• Multiple higher-order moments
• Principal components analysed geometric index

Delcaration of interest: The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.

Appendix: Matlab m-file script of the mathematical model

clc;clear all; close all;

epsval = 0.0004; % Truncation Accuracy Value

Fs = 1000; % Sampling Frequency

MeanRR = 0.8; % Mean Value of RR Intervals (Mean Heart Rate)

BeatsNumber = 500; % Number of Beats to be Generated

SNR_ECG = 30; %Signal to noise ratio for ECG

SNR_ABP = 2; %Signal to noise ratio for ABP

SNR_PCG = 35; %Signal to noise ratio for PCG

%========================================================%

% Synthetic ECG Parameters

% P–Q–R–S–T Events

Events_Location = [10 25 30 35 65]/100;

Events_Duration = [0.30, 0.05, 0.14, 0.05, 0.50];

Events_Amplitude = 10*[0.035 -0.025 0.21 -0.02 0.06];

%========================================================%

Events_Number = length(Events_Amplitude);

SumOfRR = 0;

TotalTime = MeanRR *BeatsNumber;

Allx = (0: (1/Fs): TotalTime)’;

SizeOfAllx = length(Allx);

SyntheticECG = zeros(SizeOfAllx, 1);

Lambda=zeros(1,length(Events_Amplitude));

%========================================================%

for ii = 1: BeatsNumber,

  BeatStartTime = SumOfRR;

  SumOfRR = SumOfRR + MeanRR;

  for jj = 1: Events_Number,

     Lambda(jj) = - 4*log(epsval/abs(Events_Amplitude(jj)))*(1/(Events_ Duration(jj)*MeanRR)^2);

end,

KthBeatEventsLocation = BeatStartTime + Events_ Location*MeanRR;

for jj = 1: Events_Number,

SyntheticECG = SyntheticECG + …

Events_Amplitude(jj)*exp((-(Lambda(jj)))*…

(Allx - KthBeatEventsLocation(jj)).^2);

end,

end,

%========================================================%

% Synthetic ABP Parameters

% Diastolic Pressure, Systolic Pressure, Dicrotic Notch

%========================================================%

Events_Location = [0 14 37 39 42 60 65 67 77 84 90 98 ]/100;

Events_Duration = [1, 1,5,1, 0.9, 1, 1,1,1,1,1,1];

Events_Amplitude = 10*[-5.5 7 -.5 -4 3.5 -.5 1.2 -0.1 -2.5.2 0.5 1.2];

c1=92;tt=0.1;s=.4;r=0.45;m=-100;ss=-0.2450;

%========================================================%

Events_Number = length(Events_Amplitude);

SumOfRR = 0;

TotalTime = MeanRR *BeatsNumber;

Allx = (0: (1/Fs): TotalTime)’;

SizeOfAllx = length(Allx);

SyntheticABP = zeros(SizeOfAllx, 1);

[T,Allx] = Twave(c1,tt,s,m,r,ss,MeanRR,BeatsNumber,Fs);

for ii = 1: BeatsNumber,

BeatStartTime = SumOfRR;

SumOfRR = SumOfRR + MeanRR;

for jj = 1: Events_Number,

Lambda(jj) = -4*log(epsval/abs(Events_ Amplitude(jj)))*(1/…

(Events_Duration(jj)*MeanRR)^2);

end,

KthBeatEventsLocation = BeatStartTime + Events_ Location*MeanRR;

for jj = 1: Events_Number,

SyntheticABP = SyntheticABP + …

Events_Amplitude(jj)*exp((-(Lambda(jj)))*…

(Allx-0.45 - KthBeatEventsLocation(jj)).^2);

end,

end,

SyntheticABP = T’+ SyntheticABP;

%========================================================%

% Generation of Synthetic ABP Signal with Alternans

%========================================================%

Events_Location = [0 14 37 39 42 60 65 67 77 84 90 98 ]/100;

Events_Duration = [1, 1,5,1, 0.9,1, 1,1,1,1,1,1];

Events_Amplitude = 10*[-5.5 7 -1 -4 3.5 -.5 1.2 -0.1 -2.5 .2 0.5 1.2];

c =92;tt=0.1;s=.4;r=0.45;m=-100;ss=-0.2450;

%========================================================%

Events_Number = length(Events_Amplitude);

SumOfRR = 0;

TotalTime=BeatsNumber*MeanRR;

Allx = (0: (1/Fs): TotalTime);

SizeOfAllx = length(Allx);

SyntheticABP0 = zeros(SizeOfAllx,1);

[T,Allx] = Twave(c1,tt,s,m,r,ss,MeanRR,BeatsNumber,Fs);

Lambda=zeros(1,length(Events_Amplitude));

for ii = 1: BeatsNumber,

BeatStartTime = SumOfRR;

x=[0 rand*(-1)^(2*round(rand))*3 rand*(-1)^ (2*round(rand))*12 rand*(-1)^(2*round(rand))*10 …

(-1)^(2*round(rand))*rand*10 0 0 0 0 0 0 0]; %Random Values adding to the events amplitude

y=[rand*(-1)^(round(2*rand))*0.02 rand*(-1)^ (round(2*rand))*0.02 rand*(-1)^(2*round(rand))*0.02 …

rand*(-1)^(2*round(rand))*0.02 (-1)^(2*round(rand))*rand*0.02 0 0 0 0 0 0 0]; %Random Values adding to the events location

Events_Amplitudeii = Events_Amplitude + x;

Events_Locationii = Events_Location + y;

SumOfRR = SumOfRR + MeanRR;

for jj = 1: Events_Number,

Lambda(jj) =  -4*log(epsval/abs(Events_Amplitudeii(jj)))*(1/…

(Events_Duration(jj)*MeanRR)^2);

end,

KthBeatEventsLocation = BeatStartTime + Events_ Location*MeanRR;

for jj = 1: Events_Number,

SyntheticABP0 = SyntheticABP0 + …

Events_Amplitudeii(jj)*exp((-(Lambda(jj)))*…

(Allx-0.45 - KthBeatEventsLocation(jj)).^2);

end,

end,

SyntheticABP0 = T’+ SyntheticABP0;

%=========================================================

%Finding Each Event’s Amplitude and Incidence Location

%=========================================================

N=75;

M=120;

t1=0;t2=0;t3=0;t4=0;

DBP=zeros(1,BeatsNumber); %Diastolic Blood Pressure

SBP=zeros(1,BeatsNumber); %Systolic Blood Pressure

DN=zeros(1,BeatsNumber); %Dicrotic Notch

DP=zeros(1,BeatsNumber); %Dicrotic Peak

DBPL=zeros(1,BeatsNumber); %Diastolic Blood Pressure Location

SBPL=zeros(1,BeatsNumber); %Systolic Blood Pressure Location

DNL=zeros(1,BeatsNumber); %Dicrotic Notch Location

DPL=zeros(1,BeatsNumber); %Dicrotic Peak Location

for i=1:1:length(Allx)-2

if SyntheticABP0(i+1)<N &&SyntheticABP0(i+1)<SyntheticABP0(i+2) && SyntheticABP0(i+1)<SyntheticABP0(i)

t1=t1+1;

DBP(t1)=SyntheticABP0(i+1);

DBPL(t1)=(i+1);

elseif SyntheticABP0(i+1)>M &&

SyntheticABP0(i+1)>SyntheticABP0(i+2) &&

SyntheticABP0(i+1)>SyntheticABP0(i)

t2=t2+1;

SBP(t2)=SyntheticABP0(i+1);

SBPL(t2)=(i+1);

elseif SyntheticABP0(i+1)>N &&

SyntheticABP0(i+1)<SyntheticABP0(i+2) &&

SyntheticABP0(i+1)<SyntheticABP0(i)

t3=t3+1;

DN(t3)=SyntheticABP0(i+1);

DNL(t3)=(i+1);

elseif SyntheticABP0(i+1)<M &&

SyntheticABP0(i+1)>SyntheticABP0(i+2) &&

SyntheticABP0(i+1)>SyntheticABP0(i)

t4=t4+1;

DP(t4)=SyntheticABP0(i+1);

DPL(t4)=(i+1);

end

end

DcBP=DP-DN; %Dicrotic Blood Pressure

MAP=(SBP+2*DBP)/3;

BN=1:BeatsNumber; %Beats Number

Events_AMP=[DBP’ SBP’ DcBP’ MAP’];

DPL(BeatsNumber)=BeatsNumber*Fs*MeanRR;

Events_Loc=[DBPL’ SBPL’ DNL’ DPL’];

%=========================================================

% function [T,x] = Twave(c1,t,s,ss,m,r,d,n,Fs)

% i = 0;

% de = 1/Fs;

% c2 = c1;

% c3 = m*(r - t) + c2;

%Tsize=d*n*Fs;

%T=zeros(1,Tsize);

%x=zeros(1,Tsize);

% for y = 0: de: d*n;

% i=i + 1;

% if y <= t+floor((i-1)*de/d)*d;

% T(i)=c1*exp(-(y-t-floor((i-1)*de/d)*d).^2./(2.*s.^4));

% end

% if t+floor((i-1)*de/d)*d < y;

% T(i)=m*(y-floor((i-1)*de/d)*d) + c2 - m*t;

% end

%

% if r+floor((i-1)*de/d)*d <= y;

% T(i)=c3*exp(-(y-r-floor((i-1)*de/d)*d).^2./(2.*ss.^3));

% end

%

% end

% for i = 1: (d*n/de) + 1;

% x(i) = (i-1) * de;

% end

% x = x’;

%========================================================%

%Generation of Synthetic Heart Sound Signal(PCG)%%

%========================================================%

% Primary Input for generating Synthetic PCG

FirstLocation_S4=0.161; %First Location of S4 is in 0.161s

SampleFrequency_S4 = 85;

SampleTotalTime_S4=0.1; %The Duration of each S4 signal is 0.1s

FirstLocation_S1=0.310; %First Location of S1 is in 0.310s

SampleFrequency_S1=95;

SampleTotalTime_S1=0.1; %The Duration of each S1 signal is 0.1s

FirstLocation_S2=0.625; %First Location of S2 is in 0.625s

SampleFrequency_S2=95;

SampleTotalTime_S2=0.1; %The Duration of each S2 signal is 0.1s

FirstLocation_S3=0.750; %First Location of S3 is in 0.750s

SampleFrequency_S3=85;

SampleTotalTime_S3=0.1; %The Duration of each S3 signal is 0.1s

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++%

X=0: 1/Fs: (BeatsNumber + 1)*MeanRR;

XLength = length(X);

SyntheticPCG = zeros(XLength, 1);

SyntheticPCG_0=zeros(XLength, 1);

%Note: “Synthetic PCG” stands for Synthetic Heart Sound Signal “with Alternation”

%Note: “SyntheticPCG_0” and All Variables With This Format Stands for Synthetic Heart Sound Signal “without Alternation”

%++++++++++++++++++++++++++++++++++++++++++++++++++++++++%

%========================================================%

%Random Numbers Generation for Using in Alternation

%========================================================%

Alternation_Amplitude_TotalTime = 0.02;

Alternation_Amplitude_SampleFrequency = 10;

n = BeatsNumber;

a = zeros(n, 1);

for i = 1: 8

a(((n*i)-(n-1)): (n*i)) = rand(1, n);

b(((n*i)-(n-1)): (n*i)) = rand(1, n);

for j=((n*i)-(n-1)): (n*i)

if b(j) <= 0.5

a(j) = -a(j);

end

end

end

RandNumForTotalTime_S4 = Alternation_Amplitude_TotalTime * a(0*n+1: 1*n);

RandNumForTotalTime_S1 = Alternation_Amplitude_TotalTime * a(1*n+1: 2*n);

RandNumForTotalTime_S2 = Alternation_Amplitude_TotalTime * a(2*n+1: 3*n);

RandNumForTotalTime_S3 = Alternation_Amplitude_TotalTime * a(3*n+1: 4*n);

RandNumForFrequency_S4 = Alternation_Amplitude_SampleFrequency * a(4*n+1: 5*n);

RandNumForFrequency_S1 = Alternation_Amplitude_SampleFrequency * a(5*n+1: 6*n);

RandNumForFrequency_S2 = Alternation_Amplitude_SampleFrequency * a(6*n+1: 7*n);

RandNumForFrequency_S3 = Alternation_Amplitude_SampleFrequency * a(7*n+1: 8*n);

%========================================================%

%S4 Generation & Computing Its Parameters

%========================================================%

Events_Location_S4 = [0 15 19 30 68 78 80]/100;

Events_Duration_S4 = [1,1,1,1.1,1,1,1];

Events_Amplitude_S4 = 0.1*[0.9 -2.4 -0.2 6 6.2 -3 0.1];

Events_Number_S4 = length(Events_Amplitude_S4);

SumOfRR_S4 = 0;

TotalTime_S4 = SampleTotalTime_S4 + RandNumForTotalTime_S4;

TotalTime_S4_0 = SampleTotalTime_S4;

BeatStartTime_S4 = SumOfRR_S4;

SumOfRR_S4 = SumOfRR_S4 + TotalTime_S4;

SumOfRR_S4_0 = SumOfRR_S4 + TotalTime_S4_0;

SilenceDuration_S4 = MeanRR; %Time distance between each S4 to the next one

FirstLocation_S4 = FirstLocation_S4 * Fs;

JumpS4 = FirstLocation_S4; %JumpS4 determine the location of S4 in each loop’s RUN

SDS4 = SilenceDuration_S4 * Fs;

Amplitude_S4 = zeros(BeatsNumber, 1); %The amplitude of each S4 signal will be recorded in this matrix

Duration_S4 = zeros(BeatsNumber, 1); %The duration of each S4 signal will be recorded in this matrix

First_Point_S4 = zeros(BeatsNumber, 1); %First Point of each S4 signal will be recorded in this matrix for signal detections

End_Point_S4 = zeros(BeatsNumber, 1); %End Point of each S4 signal will be recorded in this matrix for signal detections

ii = 0;

for i = 1: BeatsNumber

Allx_S4 = (0: 1/Fs: TotalTime_S4(i))’;

Allx_S4_0 = (0: 1/Fs: TotalTime_S4_0)’;

SizeOfAllx_S4 = length(Allx_S4);

SizeOfAllx_S4_0 = length(Allx_S4_0);

Push_S4 = zeros(SizeOfAllx_S4, 1);

Push_S4_0 = zeros(SizeOfAllx_S4_0, 1);

Lambda4 = zeros(Events_Number_S4, 1);

Lambda4_0 = Lambda4;

for jj = 1: Events_Number_S4,

Lambda4(jj) = - 4*log(epsval/abs(Events_Amplitude_S4(jj)))*(1/…

(Events_Duration_S4(jj)*TotalTime_S4(i))^2);

Lambda4_0(jj) = - 4*log(epsval/abs(Events_Amplitude_S4(jj)))*(1/…

(Events_Duration_S4(jj)*TotalTime_S4_0)^2);

end,

KthBeatEventsLocation4 = BeatStartTime_S4 + Events_Location_S4*TotalTime_S4(i);

KthBeatEventsLocation4_0 = BeatStartTime_S4 + Events_Location_S4*TotalTime_S4_0;

for jj = 1: Events_Number_S4,

Push_S4 = Push_S4 + …

Events_Amplitude_S4(jj)*exp((-(Lambda4(jj)))*…

(Allx_S4 - KthBeatEventsLocation4(jj)).^2);

Push_S4_0 = Push_S4_0 + …

Events_Amplitude_S4(jj)*exp((-(Lambda4_0(jj)))*…

(Allx_S4_0 - KthBeatEventsLocation4_0(jj)).^2);

end,

Frequency_S4 = (SampleFrequency_S4 + RandNumForFrequency_S4(i)) * exp(-Allx_S4).^3;

Frequency_S4_0 = SampleFrequency_S4 * exp(-Allx_S4_0).^3;

y4 = sin(2.*Frequency_S4* pi.* Allx_S4);

y4_0 = sin(2.*Frequency_S4_0* pi.* Allx_S4_0);

S4 = y4.*Push_S4;

S4_0 = y4_0.*Push_S4_0;

Length_S4 = length(S4);

Length_S4_0 = length(S4_0);

JumpS4 = JumpS4 + SDS4;

ii = ii + 1;

for jj = (JumpS4-SDS4): (JumpS4+Length_S4-SDS4-1)

SyntheticPCG(jj) = SyntheticPCG(jj) + S4(jj-(FirstLocation_S4-1) - SDS4*(ii-1));

end

for jj = (JumpS4-SDS4): (JumpS4+Length_S4_0-SDS4-1)

SyntheticPCG_0(jj) = SyntheticPCG_0(jj) + S4_0(jj-(FirstLocation_S4-1) - SDS4*(ii-1));

end

First_Point_S4(i) = (JumpS4-SDS4);

End_Point_S4(i) = (JumpS4+Length_S4-SDS4-1);

k = 1;

maxs_S4 = zeros(Length_S4 - 2, 1);

mins_S4 = zeros(Length_S4 - 2, 1);

for jj = 2: (Length_S4-1)

if S4(jj) > S4(jj-1) && S4(jj) > S4(jj+1)

maxs_S4(k) = S4(jj+1);

k = k + 1;

elseif S4(jj) < S4(jj-1) && S4(jj) < S4(jj+1)

mins_S4(k) = S4(jj+1);

k = k + 1;

end

end

MAX_S4 = max(maxs_S4);

MIN_S4 = min(mins_S4);

Amplitude_S4(i)= MAX_S4 - MIN_S4;

Duration_S4(i) = TotalTime_S4(i);

end

%========================================================%

%S1 Generation & Computing its Parameters

%========================================================%

Events_Location_S1 = [3 10 20 36 52 60 70 75 85]/100;

Events_Duration_S1 = [0.5,0.5,0.5,1,1,1,1,1,1];

Events_Amplitude_S1 = 0.12*[2 -4 5.5 3.8 1.6 1.8 0.1 0.2 0.1];

Events_Number_S1 = length(Events_Amplitude_S1);

SumOfRR_S1 = 0;

TotalTime_S1 = SampleTotalTime_S1 + RandNumForTotalTime_S1;

TotalTime_S1_0 = SampleTotalTime_S1;

BeatStartTime_S1 = SumOfRR_S1;

SumOfRR_S1 = SumOfRR_S1 + TotalTime_S1;

SumOfRR_S1_0 = SumOfRR_S1 + TotalTime_S1_0;

SilenceDuration_S1 = MeanRR; %Time Distance Between each S1 to the next one

FirstLocation_S1 = FirstLocation_S1 * Fs;

JumpS1 = FirstLocation_S1; %JumpS1 determine the location of S1 in each loop’s RUN

SDS1 = SilenceDuration_S1 * Fs;

Amplitude_S1 = zeros(BeatsNumber, 1); %The amplitude of each S1 signal will be recorded in this matrix

Duration_S1 = zeros(BeatsNumber, 1); %The duration of each S1 signal will be recorded in this matrix

First_Point_S1 = zeros(BeatsNumber, 1); %First Point of each S1 signal will be recorded in this matrix for signal detections

End_Point_S1 = zeros(BeatsNumber, 1); %End Point of each S1 signal will be recorded in this matrix for signal detections

S1_S1_Distance = zeros(BeatsNumber, 1);

ii = 0;

s = 2;

for i = 1: BeatsNumber

Allx_S1 = (0: 1/Fs: TotalTime_S1(i))’;

Allx_S1_0 = (0: 1/Fs: TotalTime_S1_0)’;

SizeOfAllx_S1 = length(Allx_S1);

SizeOfAllx_S1_0 = length(Allx_S1_0);

Push_S1 = zeros(SizeOfAllx_S1, 1);

Push_S1_0 = zeros(SizeOfAllx_S1_0, 1);

Lambda1 = zeros(Events_Number_S1, 1);

Lambda1_0 = Lambda1;

for jj = 1: Events_Number_S1,

Lambda1(jj) = - 4*log(epsval/abs(Events_Amplitude_S1(jj)))*(1/…

(Events_Duration_S1(jj)*TotalTime_S1(i))^2);

Lambda1_0(jj) = - 4*log(epsval/abs(Events_Amplitude_S1(jj)))*(1/…

(Events_Duration_S1(jj)*TotalTime_S1_0)^2);

end,

KthBeatEventsLocation1 = BeatStartTime_S1 + Events_Location_S1*TotalTime_S1(i);

KthBeatEventsLocation1_0 = BeatStartTime_S1 + Events_Location_S1*TotalTime_S1_0;

for jj = 1: Events_Number_S1,

Push_S1 = Push_S1 +…

Events_Amplitude_S1(jj)*exp((-(Lambda1(jj)))*…

(Allx_S1 - KthBeatEventsLocation1(jj)).^2);

Push_S1_0 = Push_S1_0 +…

Events_Amplitude_S1(jj)*exp((-(Lambda1_0(jj)))*…

(Allx_S1_0 - KthBeatEventsLocation1_0(jj)).^2);

end,

Frequency_S1 = (SampleFrequency_S1 + RandNumForFrequency_S1(i)) * exp(-Allx_S1).^3;

Frequency_S1_0 = SampleFrequency_S1 * exp(-Allx_S1_0).^3;

y1 = sin(2.*Frequency_S1* pi.* Allx_S1);

y1_0 = sin(2.*Frequency_S1_0* pi.* Allx_S1_0);

S1 = y1.*Push_S1;

S1_0 = y1_0.*Push_S1_0;

Length_S1 = length(S1);

Length_S1_0 = length(S1_0);

JumpS1 = JumpS1 + SDS1;

ii = ii + 1;

for jj = (JumpS1-SDS1): (JumpS1+Length_S1-SDS1-1)

SyntheticPCG(jj) = SyntheticPCG(jj) + S1(jj-(FirstLocation_S1-1) - SDS1*(ii-1));

end

for jj = (JumpS1-SDS1): (JumpS1+Length_S1_0-SDS1-1)

SyntheticPCG_0(jj) = SyntheticPCG_0(jj) + S1_0(jj-(FirstLocation_S1-1) - SDS1*(ii-1));

end

First_Point_S1(i) = (JumpS1-SDS1);

End_Point_S1(i) = (JumpS1+Length_S1-SDS1-1);

if i> = 2

S1_S1_Distance(i) = (1/Fs) * (End_Point_S1(s) - First_Point_S1(s-1));

s = s + 1;

end

k = 1;

maxs_S1 = zeros(Length_S1 - 2, 1);

mins_S1 = zeros(Length_S1 - 2, 1);

for jj = 2: (Length_S1-1)

if S1(jj) > S1(jj-1) && S1(jj) > S1(jj+1)

maxs_S1(k) = S1(jj+1);

k = k + 1;

elseif S1(jj) < S1(jj-1) && S1(jj) < S1(jj+1)

mins_S1(k) = S1(jj+1);

k = k + 1;

end

end

MAX_S1 = max(maxs_S1);

MIN_S1 = min(mins_S1);

Amplitude_S1(i)= MAX_S1 - MIN_S1;

Duration_S1(i) = TotalTime_S1(i);

end

%========================================================%

%S2 Generation & Computing Its Parameters

%========================================================%

Events_Location_S2 = [3 10 20 36 52 60 70 75 85]/100;

Events_Duration_S2 = [0.5,0.5,0.5,1,1,1,1,1,1];

Events_Amplitude_S2 = 0.12*[2 -4 5.5 3.8 1.6 1.8 0.1 0.2 0.1];

Events_Number_S2 = length(Events_Amplitude_S2);

SumOfRR_S2 = 0;

TotalTime_S2 = SampleTotalTime_S2 + RandNumForTotalTime_S2;

TotalTime_S2_0 = SampleTotalTime_S2;

BeatStartTime_S2 = SumOfRR_S2;

SumOfRR_S2 = SumOfRR_S2 + TotalTime_S2;

SumOfRR_S2_0 = SumOfRR_S2 + TotalTime_S2_0;

SilenceDuration_S2 = MeanRR; %Time Distance Between each S2 to the next one

FirstLocation_S2 = FirstLocation_S2 * Fs;

JumpS2 = FirstLocation_S2; %JumpS2 determine the location of S2 in each loop’s RUN

SDS2 = SilenceDuration_S2 * Fs;

Amplitude_S2 = zeros(BeatsNumber, 1); %The amplitude of each S2 signal will be recorded in this matrix

Duration_S2 = zeros(BeatsNumber, 1); %The duration of each S2 signal will be recorded in this matrix

First_Point_S2 = zeros(BeatsNumber, 1); %First Point of each S2 signal will be recorded in this matrix for signal detections

End_Point_S2 = zeros(BeatsNumber, 1); %End Point of each S2 signal will be recorded in this matrix for signal detections

S2_S2_Distance = zeros(BeatsNumber, 1);

ii = 0;

r = 2;

for i = 1: BeatsNumber

Allx_S2 = (0: 1/Fs: TotalTime_S2(i))’;

Allx_S2_0 = (0: 1/Fs: TotalTime_S2_0)’;

SizeOfAllx_S2 = length(Allx_S2);

SizeOfAllx_S2_0 = length(Allx_S2_0);

Push_S2 = zeros(SizeOfAllx_S2, 1);

Push_S2_0 = zeros(SizeOfAllx_S2_0, 1);

Lambda2 = zeros(Events_Number_S2, 1);

Lambda2_0 = Lambda2;

for jj = 1: Events_Number_S2,

Lambda2(jj) =  - 4*log(epsval/abs(Events_Amplitude_S2(jj)))*(1/…

(Events_Duration_S2(jj)*TotalTime_S2(i))^2);

Lambda2_0(jj) =  - 4*log(epsval/abs(Events_Amplitude_S2(jj)))*(1/…

(Events_Duration_S2(jj)*TotalTime_S2_0)^2);

end,

KthBeatEventsLocation2 = BeatStartTime_S2 + Events_Location_S2*TotalTime_S2(i);

KthBeatEventsLocation2_0 = BeatStartTime_S2 + Events_Location_S2*TotalTime_S2_0;

for jj = 1: Events_Number_S2,

Push_S2 = Push_S2 + …

Events_Amplitude_S2(jj)*exp((-(Lambda2(jj)))*…

(Allx_S2 - KthBeatEventsLocation2(jj)).^2);

Push_S2_0 = Push_S2_0 + …

Events_Amplitude_S2(jj)*exp((-(Lambda2_0(jj)))*…

(Allx_S2_0 - KthBeatEventsLocation2_0(jj)).^2);

end,

Frequency_S2 = (SampleFrequency_S2 + RandNumForFrequency_S2(i)) * exp(-Allx_S2).^3;

Frequency_S2_0 = SampleFrequency_S2 * exp(-Allx_S2_0).^3;

y2 = sin(2.*Frequency_S2* pi.* Allx_S2);

y2_0 = sin(2.*Frequency_S2_0* pi.* Allx_S2_0);

S2 = y2.*Push_S2;

S2_0 = y2_0.*Push_S2_0;

Length_S2 = length(S2);

Length_S2_0 = length(S2_0);

JumpS2 = JumpS2 + SDS2;

ii = ii + 1;

for jj = (JumpS2-SDS2): (JumpS2+Length_S2-SDS2-1)

SyntheticPCG(jj) = SyntheticPCG(jj) + S2(jj-(FirstLocation_S2-1) - SDS2*(ii-1));

end

for jj = (JumpS2-SDS2): (JumpS2+Length_S2_0-SDS2-1)

SyntheticPCG_0(jj) = SyntheticPCG_0(jj) + S2_0(jj-(FirstLocation_S2-1) - SDS2*(ii-1));

end

First_Point_S2(i) = (JumpS2-SDS2);

End_Point_S2(i) = (JumpS2+Length_S2-SDS2-1);

if i >= 2

S2_S2_Distance(i) = (1/Fs) * (End_Point_S2(r) - First_Point_S2(r-1));

r = r + 1;

end

k = 1;

maxs_S2 = zeros(Length_S2 - 2, 1);

mins_S2 = zeros(Length_S2 - 2, 1);

for jj = 2: (Length_S2-1)

if S2(jj) > S2(jj-1) && S2(jj) > S2(jj+1)

maxs_S2(k) = S2(jj+1);

k = k + 1;

elseif S2(jj) < S2(jj-1) && S2(jj) < S2(jj+1)

mins_S2(k) = S2(jj+1);

k = k + 1;

end

end

MAX_S2 = max(maxs_S2);

MIN_S2 = min(mins_S2);

Amplitude_S2(i)= MAX_S2 - MIN_S2;

Duration_S2(i) = TotalTime_S2(i);

end

%========================================================%

%S3 Generation & Computing Its Parameters

%========================================================%

Events_Location_S3 = [0 15 19 30 68 78 80]/100;

Events_Duration_S3 = [1,1,1,1.1,1,1,1];

Events_Amplitude_S3 = 0.1*[0.9 -2.4 -0.2 6 6.2 -3 0.1];

Events_Number_S3 = length(Events_Amplitude_S3);

SumOfRR_S3 = 0;

TotalTime_S3 = SampleTotalTime_S3 + RandNumForTotalTime_S3;

TotalTime_S3_0 = SampleTotalTime_S3;

BeatStartTime_S3 = SumOfRR_S3;

SumOfRR_S3 = SumOfRR_S3 + TotalTime_S3;

SumOfRR_S3_0 = SumOfRR_S3 + TotalTime_S3_0;

SilenceDuration_S3 = MeanRR; %Time Distance Between each S3 to the next one

FirstLocation_S3 = FirstLocation_S3 * Fs;

JumpS3 = FirstLocation_S3; %JumpS3 determine the location of S3 in each loop’s RUN

SDS3 = SilenceDuration_S3 * Fs;

Amplitude_S3 = zeros(BeatsNumber, 1); %The amplitude of each S2 signal will be recorded in this matrix

Duration_S3 = zeros(BeatsNumber, 1); %The duration of each S2 signal will be recorded in this matrix

First_Point_S3 = zeros(BeatsNumber, 1); %First Point of each S2 signal will be recorded in this matrix for signal detections

End_Point_S3 = zeros(BeatsNumber, 1); %End Point of each S2 signal will be recorded in this matrix for signal detections

ii = 0;

for i = 1: BeatsNumber

Allx_S3 = (0: 1/Fs: TotalTime_S3(i))’;

Allx_S3_0 = (0: 1/Fs: TotalTime_S3_0)’;

SizeOfAllx_S3 = length(Allx_S3);

SizeOfAllx_S3_0 = length(Allx_S3_0);

Push_S3 = zeros(SizeOfAllx_S3, 1);

Push_S3_0 = zeros(SizeOfAllx_S3_0, 1);

Lambda3 = zeros(Events_Number_S3, 1);

Lambda3_0 = Lambda3;

for jj = 1: Events_Number_S3,

Lambda3(jj) = - 4*log(epsval/abs(Events_Amplitude_S3(jj)))*(1/…

(Events_Duration_S3(jj)*TotalTime_S3(i))^2);

Lambda3_0(jj) = - 4*log(epsval/abs(Events_Amplitude_S3(jj)))*(1/…

(Events_Duration_S3(jj)*TotalTime_S3_0)^2);

end,

KthBeatEventsLocation3 = BeatStartTime_S3 + Events_Location_S3*TotalTime_S3(i);

KthBeatEventsLocation3_0 = BeatStartTime_S3 + Events_Location_S3*TotalTime_S3_0;

for jj = 1: Events_Number_S3,

Push_S3 = Push_S3 +…

Events_Amplitude_S3(jj)*exp((-(Lambda3(jj)))*…

(Allx_S3 - KthBeatEventsLocation3(jj)).^2);

Push_S3_0 = Push_S3_0 +…

Events_Amplitude_S3(jj)*exp((-(Lambda3_0(jj)))*…

(Allx_S3_0 - KthBeatEventsLocation3_0(jj)).^2);

end,

Frequency_S3 = (SampleFrequency_S3 + RandNumForFrequency_S3(i)) * exp(-Allx_S3).^3;

Frequency_S3_0 = SampleFrequency_S3 * exp(-Allx_S3_0).^3;

y3 = sin(2.*Frequency_S3* pi.* Allx_S3);

y3_0 = sin(2.*Frequency_S3_0* pi.* Allx_S3_0);

S3 = y3.*Push_S3;

S3_0 = y3_0.*Push_S3_0;

Length_S3 = length(S3);

Length_S3_0 = length(S3_0);

JumpS3 = JumpS3 + SDS3;

ii = ii + 1;

for jj = (JumpS3-SDS3): (JumpS3+Length_S3-SDS3-1)

SyntheticPCG(jj) = SyntheticPCG(jj) + S3(jj-(FirstLocation_S3-1) - SDS3*(ii-1));

end

for jj = (JumpS3-SDS3): (JumpS3+Length_S3_0-SDS3-1)

SyntheticPCG_0(jj) = SyntheticPCG_0(jj) + S3_0(jj-(FirstLocation_S3-1) - SDS3*(ii-1));

end

First_Point_S3(i) = (JumpS3-SDS3);

End_Point_S3(i) = (JumpS3+Length_S3-SDS3-1);

k = 1;

maxs_S3 = zeros(Length_S3 - 2, 1);

mins_S3 = zeros(Length_S3 - 2, 1);

for jj = 2: (Length_S3-1)

if S3(jj) > S3(jj-1) && S3(jj) > S3(jj+1)

maxs_S3(k) = S3(jj+1);

k = k + 1;

elseif S3(jj) < S3(jj-1) && S3(jj) < S3(jj+1)

mins_S3(k) = S3(jj+1);

k = k + 1;

end

end

MAX_S3 = max(maxs_S3);

MIN_S3 = min(mins_S3);

Amplitude_S3(i)= MAX_S3 - MIN_S3;

Duration_S3(i) = TotalTime_S3(i);

end

%======================================================%

S1_S2_Distance = (End_Point_S2 - End_Point_S1)/Fs;

S1_S1_Distance = S1_S1_Distance’;

S2_S2_Distance = S2_S2_Distance’;

S1_S2_Distance = S1_S2_Distance’;

%========================================================%

% Contaminating the Synthetic ECG with Noise and Baseline Wander

SyntheticECG_Noise = awgn(SyntheticECG,SNR_ECG);

%========================================================%

% Contaminating the Synthetic ABP with Noise and Baseline Wander

SyntheticABP_Noise = awgn(SyntheticABP,SNR_ABP);

%========================================================%

% Contaminating the Synthetic PCG with Noise and Baseline Wander

SyntheticPCG_Noise = awgn(SyntheticPCG,SNR_PCG);

SyntheticPCG_Noise_0 = awgn(SyntheticPCG_0,SNR_PCG);

%========================================================%

figure (1)

subplot(3,1,1);plot(Allx,SyntheticECG_Noise);axis([4.85 7.5 -.500 2.500]);xlabel(‘Time(sec)’);ylabel(‘SyntheticECG’)

subplot(3,1,2);plot(Allx,SyntheticABP_Noise);axis([4.85 7.5 50 145]);xlabel(‘Time(sec)’);ylabel(‘SyntheticABP’)

subplot(3,1,3);plot(X, SyntheticPCG_Noise_0);axis([4.85 7.5 -1 1]);xlabel(‘Time(sec)’);ylabel(‘SyntheticPCG’)

figure (2)

plot(Allx,SyntheticABP0);xlabel(‘Time(sec)’);ylabel(‘SyntheticABP with Alternance’);hold on

plot(Events_Loc(:,1)/1000,SyntheticABP0(Events_Loc(:,1)),’rp’);hold on

plot(Events_Loc(:,2)/1000,SyntheticABP0(Events_Loc(:,2)),’go’);hold on

plot(Events_Loc(:,3)/1000,SyntheticABP0(Events_Loc(:,3)),’m^’);hold on

plot(Events_Loc(:,4)/1000,SyntheticABP0(Events_Loc(:,4)),’cs’);hold on

figure (3)

plot(BN,DBP);axis([0 BeatsNumber -10 160]);xlabel(‘Time(sec)’);ylabel(‘SyntheticABP Events Amplitude’);hold on

plot(BN,SBP);axis([0 BeatsNumber -10 160]);hold on

plot(BN,DcBP);axis([0 BeatsNumber -10 160]);hold on

plot(BN,MAP);axis([0 BeatsNumber -10 160]);hold on

figure(4)

plot(X, SyntheticPCG), title(‘Synthetic PCG with Alternation’), xlabel(‘time(s)’);hold on

plot(First_Point_S4 / Fs, SyntheticPCG(First_Point_S4), ‘rs’, ‘MarkerFaceColor’,’r’);hold on

plot(End_Point_S4 / Fs, SyntheticPCG(End_Point_S4), ‘gs’, ‘MarkerFaceColor’,’g’);hold on

plot(First_Point_S1 / Fs, SyntheticPCG(First_Point_S1), ‘rs’, ‘MarkerFaceColor’,’r’);hold on

plot(End_Point_S1 / Fs, SyntheticPCG(End_Point_S1), ‘gs’, ‘MarkerFaceColor’,’g’);hold on

plot(First_Point_S2 / Fs, SyntheticPCG(First_Point_S2), ‘rs’, ‘MarkerFaceColor’,’r’);hold on

plot(End_Point_S2 / Fs, SyntheticPCG(End_Point_S2), ‘gs’, ‘MarkerFaceColor’,’g’);hold on

plot(First_Point_S3 / Fs, SyntheticPCG(First_Point_S3), ‘rs’, ‘MarkerFaceColor’,’r’);hold on

plot(End_Point_S3 / Fs, SyntheticPCG(End_Point_S3), ‘gs’, ‘MarkerFaceColor’,’g’)

figure(5)

subplot(2,1,1), plot(1: BeatsNumber, Amplitude_S4), title(‘S4′), xlabel(‘Beat Number’), ylabel(‘Signal Amplitude’)

subplot(2,1,2), plot(1: BeatsNumber, Duration_S4), xlabel(‘Beat Number’), ylabel(‘Signal Duration’)

figure(6)

subplot(2,1,1), plot(1: BeatsNumber, Amplitude_S1), title(‘S1’), xlabel(‘Beat Number’), ylabel(‘Signal Amplitude’)

subplot(2,1,2), plot(1: BeatsNumber, Duration_S1), xlabel(‘Beat Number’), ylabel(‘Signal Duration’)

figure(7)

subplot(2,1,1), plot(1: BeatsNumber, Amplitude_S2), title(‘S2’), xlabel(‘Beat Number’), ylabel(‘Signal Amplitude’)

subplot(2,1,2), plot(1: BeatsNumber, Duration_S2), xlabel(‘Beat Number’), ylabel(‘Signal Duration’)

figure(8)

subplot(2,1,1), plot(1: BeatsNumber, Amplitude_S3), title(‘S3’), xlabel(‘Beat Number’), ylabel(‘Signal Amplitude’)

subplot(2,1,2), plot(1: BeatsNumber, Duration_S3), xlabel(‘Beat Number’), ylabel(‘Signal Duration’)

figure(9)

subplot(3,1,1), plot(1: BeatsNumber, S1_S1_Distance), xlabel(‘Beat Number’), ylabel(‘S1-S1 Distance’), axis([0 BeatsNumber 0.85 0.95])

subplot(3,1,2), plot(1: BeatsNumber, S2_S2_Distance), xlabel(‘Beat Number’), ylabel(‘S2-S2 Distance’), axis([0 BeatsNumber 0.85 0.95])

subplot(3,1,3), plot(1: BeatsNumber, S1_S2_Distance), xlabel(‘Beat Number’), ylabel(‘S1-S2 Distance’), axis([0 BeatsNumber 0.25 0.4]