how to calculate area integral using cumsum under a known matrix of peaks and it +-50ms surrounding

Hi guys,

I'm working on quite big medical data project and I got stuck on a problem...

Here it is:

I've got matrix of the peaks values on the y-axis, name 'amp' 1x81 dimension and I've got matrix of the peaks location on the x-axis, name 'pos' also 1x81 dimension.

I would like to find area, under every peak using function cumsum. Ok, it's easy for me, I use Y = cumsum(amp,1); -> ready but there's another problem:

I would like to have +-50ms frame from every peak position and then calculate the integral area.

So basically the question is how to obtain M matrix of +-50ms surroundings of the peaks including peaks.

Guys, anyone know how to do it quickly?

I've got some idea to subtrack and add 50 ms to indeks like: a = indeks(1:end)-50 && b = indeks(1:end)+50 and then merged it somehow into a frame and then calculate integral area.

Please give me some ideas.

Below is picture demonstration how it should look like.. Thanks in advance for any help!

Answer 1

[Star Strider]

You have your sampling frequency ‘Fs’ in Hz (and sampling interval ‘Ts’ in seconds), so I would use the Signal Processing Toolbox findpeaks function to find the Lead II R-waves, then calculate the window as ±0.05/Ts to get the ‘window’ in index units you want on each side of the R-wave peak. Use findpeaks on your signal without specifying the sampling frequency so you get the peak locations in index values.

The normal QRS complex is <120 ms in width, so you could increase the window from 50 ms to 55 or 60 ms. The Lead II EKG is normally isoelectric in the P-R and S-T intervals, so including a larger window should have no effect on the integral in a normal EKG.

Then use the trapz function to do the integration.

I am not aware of a particular reason to do the integral of any part of the EKG waveform, so I am interested in what you’re studying.

[Natalia Opacka]

Yes, I'm extracting ecg derived respiration(EDR) from ECG and one of the method is to integral the area under the R peak.

But the thing is that I already found R peaks along with their index.

So I have peak locations in index values but I would like to enlarge it to +-50ms before and after each peak and then use for example: Y = cumsum(result,1);

Moreover, the result should have been the same size matrix as my ' Rpeaks' and ' index' so I can further process it with function spline.

I do attach example data from what I have now: please use 'load mydata.mat'

Please not that Fs = 200Hz, I would like to make: yy = spline(x,Y,xx);

I've got: Y=cumsum(ecg_from_database(index(1:end)-50:index(1:end):index(1:end)+50)) but it give me one value instead of 1x81 dimension matrix

[Natalia Opacka]

Ok, I think I found solution by using for loop:

for i=1:1:length(index) Y(i)=cumsum(ecg_from_database(index(i)-50:index(i):index(i)+50));

end

Maybe it will be useful for someone

Anyway, thanks for your help :)

[Star Strider]

I would not use a spline approximation or cumsum. Use trapz or cumtrapz to get an accurate integral value. The cumtrapz vector will be the same length as the original vector.

This indexing is wrong:

index(i)-50:index(i):index(i)+50)

With Fs=200Hz,Ts=0.005, the window will be 0.05/Ts=0.05*Fs or ±10 index range. So the indexing should be:

index(i)-10:index(i)+10)

‘Respiratory sinus arrythmia’ seems to have been well-researched and well-documented. See for example ‘Influences of breathing patterns on respiratory sinus arrhythmia in humans during exercise’ American Journal of Physiology - Heart and Circulatory Physiology 1 February 2005 Vol. 288 no. 2, H887-H895.

I apologise for being brutally honest, but it will not be useful to anyone because your entire approach is wrong. A more appropriate experimental design would be to compare the EKG and ear-lobe photoplythysmography if you want to show the relationship of blood flow and the respiratory cycle. Changes in the R-wave amplitude or integrated area will not change in the respiratory cycle, because the EKG component amplitudes are determined by the myocardial mass (and health), not minute-to-minute heart rate or stroke volume.

Oh, well...

Answer 2

[Thorsten]

Your solution is wrong, because you use index(i) as the increment; you probably want

ecg_from_database(index(i)-50:index(i)+50)

and you need the final value of the cumsum

cs = cumsum(ecg_from_database(index(i)-50:index(i)+50));
Y(i) = cs(end);

and consider using trapz to approximate the integration:

Y(i) = trapz(ecg_from_database(index(i)-50:index(i)+50));;

[Natalia Opacka]

Well thank you for that answer, but I'm still quite not sure which one is correct because they both gave me finally the same result. ;)

I will consult it with my professor and update.

[Thorsten]

"if the function you are integrating is well-approximated by a cubic polynomial on every subinterval" the mid-point rule (i.e., cumsum) gives a better approximation than the trapezoidal rule; for your data it is probably better to use cumsum.

 http://math.stackexchange.com/questions/603830/why-does-trapezoidal-rule-have-potential-error-greater-than-midpoint