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

Question

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

Solution 1

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]);

Solution 2

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.