QRS detection(peaks) of a raw ecg signal in matlab

2
votes

I want to find the peaks of the raw ecg signal so that I can calculate the beats per minute(bpm). I Have written a code in matlab which I have attached below.In the code below I am unable to find threshold point correctly which will help me in finding the peaks and hence the bpm.

%input the signal into matlab

[x,fs]=wavread('heartbeat.wav');
subplot(2,1,1)
plot(x(1:10000),'r-')
grid on



%lowpass filter the input signal with cutoff at 100hz
h=fir1(30,0.3126);     %normalized cutoff freq=0.3126
y=filter(h,1,x);

subplot(2,1,2)
plot(y(1:10000),'b-')
grid on



% peaks are seen as pulses(heart beats)
beat_count=0;
for p=2:length(y)-1
    th(p)=abs(max(y(p)));
    if(y(p) >y(p-1) && y(p) >y(p+1) && y(p)>th(p))
        beat_count=beat_count+1;
    end
end

N = length(y);
duration_seconds=N/fs;
duration_minutes=duration_seconds/60;
BPM=beat_count/duration_minutes;
bpm=ceil(BPM);

Please help me as I am new to matlab

3
matlabsignal-processing
Hi, I suggest that you explain the logic behind your peak detection algorithm rather than expecting people to understand it just from the code. Also abs(max(y(p))) doesn't make a lot of sense as y(p) is a scalar so max(y(p)) is just y(p), perhaps you meant max(y)? – Dan
@Dan yup...i meant by max(y)....By finding the threshold ,I can compare a value of y with its previous value,its next value and the threshold value.if the value is greater than all the 3 conditions then it is considered as a peak point. – R J
100 Hz is a crazy high cutoff for filtering QRSs from ECG data. You probably want more like about 15-30 Hz. You also want a high-pass filter with a cut-off of about 5 Hz to get rid of baseline drift etc. Further to that, you should probably be using filtfilt so there is no phase distortion - after all, in this case you ARE interested in what the time history looks like. You might also consider using the Pan & Tompkins algorithm, which works very well. – wakjah
i am also applying same thing in ios i hope you got the solution but there are two ways either you find the slope or you get the distance from the old coordinates to the new , i am also glad if i get help in ios – Waleed Ahmed

3 Answers

2
votes

I suggest changing this section of your code

beat_count=0;
for p=2:length(y)-1
    th(p)=abs(max(y(p)));
    if(y(p) >y(p-1) && y(p) >y(p+1) && y(p)>th(p))
        beat_count=beat_count+1;
    end
end

This is definitely flawed. I'm not sure of your logic here but what about this. We are looking for peaks, but only the high peaks, so first lets set a threshold value (you'll have to tweak this to a sensible number) and cull everything below that value to get rid of the smaller peaks:

th = max(y) * 0.9; %So here I'm considering anything less than 90% of the max as not a real peak... this bit really depends on your logic of finding peaks though which you haven't explained
Yth = zeros(length(y), 1);
Yth(y > th) = y(y > th);

OK so I suggest you now plot y and Yth to see what that code did. Now to find the peaks my logic is we are looking for local maxima i.e. points at which the first derivative of the function change from being positive to being negative. So I'm going to find a very simple numerical approximation to the first derivative by finding the difference between each consecutive point on the signal:

Ydiff = diff(Yth);

No I want to find where the signal goes from being positive to being negative. So I'm going to make all the positive values equal zero, and all the negative values equal one:

Ydiff_logical = Ydiff < 0;

finally I want to find where this signal changes from a zero to a one (but not the other way around)

Ypeaks = diff(Ydiff_logical) == 1;

Now count the peaks:

sum(Ypeaks)

note that for plotting purpouse because of the use of diff we should pad a false to either side of Ypeaks so 

Ypeaks = [false; Ypeaks; false];

OK so there is quite a lot of matlab there, I suggest you run each line, one by one and inspect the variable by both plotting the result of each line and also by double clicking the variable in the matlab workspace to understand what is happening at each step.

Example: (signal PeakSig taken from http://www.mathworks.com/help/signal/ref/findpeaks.html) and plotting with:

plot(x(Ypeaks),PeakSig(Ypeaks),'k^','markerfacecolor',[1 0 0]);

enter image description here

2
votes

What do you think about the built-in 

findpeaks(data,'Name',value)

function? You can choose among different logics for peak detection:

'MINPEAKHEIGHT'
'MINPEAKDISTANCE'
'THRESHOLD'
'NPEAKS'
'SORTSTR'

I hope this helps.

0
votes

You know, the QRS complex does not always have the maximum amplitude, for pathologic ECG it can be present as several minor oscillations instead of one high-amplitude peak.

Thus, you can try one good algothythm, tested by me: the detection criterion is assumed to be high absolute rate of change in the signal, averaged within the given interval.

Algorithm: - 50/60 Hz filter (e.g. for 50 Hz sliding window of 20 msec will be fine) - adaptive hipass filter (for baseline drift) - find signal's first derivate x' - fing squared derivate (x')^2 - apply sliding average window with the width of QRS complex - approx 100-150 msec (you will get some signal with 'rectangles', which have width of QRS) - use simple threshold (e.g. 1/3 of maximum of the first 3 seconds) to determine approximate positions or R - in the source ECG find local maximum within +-100 msec of that R position.

However, you still have to eliminate artifacts and outliers (e.g. surges, when the electrod connection fails).

Also, you can find a lot of helpful information from this book: "R.M. Rangayyan - Biomedical Signal Analysis"