How to find Peak to Peak delay?

How to find peak to peak delay of below shown ECG signals?

The code to plot above ECG signal is

load 16786m; figure(); plot(val(1,:));

I also want to find peak to peak delay of another ECG signal shown below:

The code to plot above signal is

load 100m; figure(); plot(val(1,:));

Source Codes are available here to download.

Will the same code detect peak distances for both ECGs?

[dpb]

What Toolboxen you got at hand?

[Explorer]

I have following tool boxes installed in my MATLAB.

MATLAB Version 8.1 (R2013a)

Simulink Version 8.1 (R2013a)

Aerospace Blockset Version 3.11 (R2013a)

Aerospace Toolbox Version 2.11 (R2013a)

Bioinformatics Toolbox Version 4.3 (R2013a)

Communications System Toolbox Version 5.4 (R2013a)

Computer Vision System Toolbox Version 5.2 (R2013a)

Control System Toolbox Version 9.5 (R2013a)

Curve Fitting Toolbox Version 3.3.1 (R2013a)

DO Qualification Kit Version 2.1 (R2013a)

DSP System Toolbox Version 8.4 (R2013a)

Data Acquisition Toolbox Version 3.3 (R2013a)

Database Toolbox Version 4.1 (R2013a)

Datafeed Toolbox Version 4.5 (R2013a)

Econometrics Toolbox Version 2.3 (R2013a)

Embedded Coder Version 6.4 (R2013a)

Filter Design HDL Coder Version 2.9.3 (R2013a)

Financial Instruments Toolbox Version 1.1 (R2013a)

Financial Toolbox Version 5.1 (R2013a)

Fixed-Point Designer Version 4.0 (R2013a)

Fuzzy Logic Toolbox Version 2.2.17 (R2013a)

Global Optimization Toolbox Version 3.2.3 (R2013a)

HDL Coder Version 3.2 (R2013a)

HDL Verifier Version 4.2 (R2013a)

IEC Certification Kit Version 3.1 (R2013a)

Image Acquisition Toolbox Version 4.5 (R2013a)

Image Processing Toolbox Version 8.2 (R2013a)

Instrument Control Toolbox Version 3.3 (R2013a)

MATLAB Builder EX Version 2.3.1 (R2013a)

MATLAB Builder JA Version 2.2.6 (R2013a)

MATLAB Builder NE Version 4.1.3 (R2013a)

MATLAB Coder Version 2.4 (R2013a)

MATLAB Compiler Version 4.18.1 (R2013a)

MATLAB Distributed Computing Server Version 6.2 (R2013a)

MATLAB Report Generator Version 3.14 (R2013a)

Mapping Toolbox Version 3.7 (R2013a)

Model Predictive Control Toolbox Version 4.1.2 (R2013a)

Model-Based Calibration Toolbox Version 4.6 (R2013a)

Neural Network Toolbox Version 8.0.1 (R2013a)

OPC Toolbox Version 3.2 (R2013a)

Optimization Toolbox Version 6.3 (R2013a)

Parallel Computing Toolbox Version 6.2 (R2013a)

Partial Differential Equation Toolbox Version 1.2 (R2013a)

Phased Array System Toolbox Version 2.0 (R2013a)

RF Toolbox Version 2.12 (R2013a)

Real-Time Windows Target Version 4.2 (R2013a)

Robust Control Toolbox Version 4.3 (R2013a)

Signal Processing Toolbox Version 6.19 (R2013a)

SimBiology Version 4.3 (R2013a)

SimDriveline Version 2.4 (R2013a)

SimElectronics Version 2.3 (R2013a)

SimEvents Version 4.3 (R2013a)

SimHydraulics Version 1.12 (R2013a)

SimMechanics Version 4.2 (R2013a)

SimPowerSystems Version 5.8 (R2013a)

SimRF Version 4.0 (R2013a)

Simscape Version 3.9 (R2013a)

Simulink 3D Animation Version 6.3 (R2013a)

Simulink Code Inspector Version 1.3 (R2013a)

Simulink Coder Version 8.4 (R2013a)

Simulink Control Design Version 3.7 (R2013a)

Simulink Design Optimization Version 2.3 (R2013a)

Simulink Design Verifier Version 2.4 (R2013a)

Simulink PLC Coder Version 1.5 (R2013a)

Simulink Report Generator Version 3.14 (R2013a)

Simulink Verification and Validation Version 3.5 (R2013a)

Spreadsheet Link EX Version 3.1.7 (R2013a)

Stateflow Version 8.1 (R2013a)

Statistics Toolbox Version 8.2 (R2013a)

Symbolic Math Toolbox Version 5.10 (R2013a)

System Identification Toolbox Version 8.2 (R2013a)

SystemTest Version 2.6.5 (R2013a)

Trading Toolbox Version 1.0 (R2013a)

Vehicle Network Toolbox Version 2.0 (R2013a)

Wavelet Toolbox Version 4.11 (R2013a)

xPC Target Version 5.4 (R2013a)

xPC Target Embedded Option Version 5.4 (R2013a)

Answer 1

[Star Strider]

Using my own code:

D = load('16786m.mat');                                 % Read EKG File
EKG = D.val/200;                                        % Gain From the ‘info’ File
Fs = 128;                                               % Fs From the ‘info’ File
tv = linspace(0,size(EKG,2),size(EKG,2))/Fs;            % Time Vector

[pks,t] = findpeaks(EKG(1,:), Fs, 'MinPeakHeight',1);   % See ‘findpeaks’ Documentation
RR = diff(t);                                           % RR Intervals (sec)
BBRate = 1./RR;                                         % Beat-To-Beat Rate (Beats/sec)
Rate = mean(BBRate);                                    % Mean Rate (Beats/sec)

figure(1)
plot(tv,EKG(1,:))                                       % Plot EKG
hold on
plot(t, pks, '+r')                                      % Plot R Peaks
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude (mV)')
legend('EKG', 'R', 'Location','SE')

[Explorer]

Thanks for replying but I am getting error with your code.

[Star Strider]

The findpeaks funciton changed since R2013a. Replace the original findpeaks single-line call with these two lines to get the equivalent ‘R2013a’ behaviour:

[pks,idx] = findpeaks(EKG(1,:), 'MinPeakHeight',1);
t = tv(idx);

I ran this and it works (in R2015b), although you might also have to rename 'MinPeakHeight' with 'MINPEAKHEIGHT' for R2013a.

[Explorer]

Thanks! Code works fine on first ECG signal but it's not working on second one (100m.mat). R Peaks are not detected correctly.

Is it possible to make a single code of peak delay measurement for both the above mentioned ECG signals?

[Star Strider]

My pleasure!

There’s a ‘Gain’ scaling value in the .info file:

Source: record nsrdb/16786  Start: [11:48:00.000]
val has 2 rows (signals) and 1280 columns (samples/signal)
Duration:     0:10
Sampling frequency: 128 Hz  Sampling interval: 0.0078125 sec
Row  Signal  Gain  Base  Units
1  ECG1  200  0  mV
2  ECG2  200  0  mV

To convert from raw units to the physical units shown
above, subtract 'base' and divide by 'gain'.

For 16786, it’s 200, but it may change for others. One possibility for making making a ‘universal’ code might be to do something like this:

Gain = 200;                                             % Gain From the ‘info’ File
Base = 0;                                               % Baseline Value From ‘info’ File
EKG = (D.val-Base)/Gain;                                % Base & Gain Scaled EKG
Fs = 128;                                               % Fs From the ‘info’ File
EKG_1 = EKG(1,:);                                       % Choose Trace To Study
. . .
[pks,idx] = findpeaks(EKG_1, 'MinPeakHeight',max(EKG_1)/3);

See if the ‘1/Gain’ threshold works. (It would work in 16786.) There are several other options with findpeaks, but EKGs can vary considerably in morphology and RR intervals, so I hesitate to use other than the 'MinPeakHeioght' threshold values. I chose to divide the maximum by 3 in order to correct for any electrical alternans that may exist. Choose the value that works best for you. Any value between 2 and 3 should work. If there is considerable baseline drift, it will be necessary for you to design a Butterworth bandpass filter with a passband of about [5 45] Hz and a stopband of about [1 50] Hz. See How to design a lowpass filter for ocean wave data in Matlab? for my filter design and implementation procedure.

—————————————————————————————————————————————————

EDIT — GetEKGinfo

I’m certain this can be made more efficient, but it’s the best I can come up with just now. (At least some of the inefficiency is the result of my needing to set different variables in order to troubleshoot it, and some of it is just due to fatigue at this point.) It was quite definitely an educational experience for me!

It gets the information from the .info files and returns necessary information for plotting and analysing the EKG records. I only tested it on ‘16786m.info’, but if the others are all of similar format, it should work for them as well. The idea is to make my code more flexible, so that you can use it without having to look up and extract the .info data manually.

First, you need to save it in your MATLAB path as GetEKGinfo.m.

To use it in my code, the function call is:

[NrEKGs,EKG_RecLen,Gain,Base,Fs] = GetEKGinfo('16786m')

where ‘NrEKGs’ is the number of EKG records in the .mat file, ‘EKG_RecLen’ is the length of each EKG, and ‘Gain’, ‘Base’ and ‘Fs’ for each have already been defined. You can call the function in my code, and then assign the returned variables as necessary.

      function [recnrs,reclen,gain,base,Fs] = GetEKGinfo(info_file_prefix)
      fidi = fopen([info_file_prefix '.info']);
      k1 = 0;
      while ~feof(fidi)
          k1 = k1 + 1;
          Line{k1} = fgets(fidi);
          units1 = strfind(Line{k1}, 'val has');
          if ~isempty(units1)
              [recst1,recend1] = regexp(Line{k1}, '\d*');
              recnrs = str2num(Line{k1}(recst1(1):recend1(1)+1));
              reclen = str2num(Line{k1}(recst1(2):recend1(2)));
          end
          units2 = strfind(Line{k1}, 'Sampling frequency:');
          if ~isempty(units2)
              [recst2,recend2] = regexp(Line{k1}, '\d*');
              Fs = str2num(Line{k1}(recst2(1):recend2(1)+1));
          end
          units3 = strfind(Line{k1}, 'Row');
          if ~isempty(units3)
              for k2 = 1:recnrs
                  k1 = k1 + k2;
                  Line{k1} = fgets(fidi);
                  [recst3,recend3] = regexp(Line{k1}, '\d*');
                  gain(k2) = str2num(Line{k1}(recst3(3):recend3(3)));
                  base(k2) = str2num(Line{k1}(recst3(4):recend3(4)));
              end
          end
      end
      fclose(fidi);
      end