Matlab Program - Calculating Beat Rate of Heart Using Discrete Wavelet Transform

%Calculating Beat Rate of heart Using Discrete Wavelet Transform
% ECGDATA is the input ecg data file
clear;close all;clc;
load ECGDATA

%Eliminate Baseline Drift
s1=ECG_1;
s2=smooth(s1,150);
ecgsmooth=s1-s2;

%Apply Wavelet Transform
[C,L]=wavedec(ecgsmooth,8,'db4');
[d1,d2,d3,d4,d5,d6,d7,d8]=detcoef(C,L,[1,2,3,4,5,6,7,8]);

%Denoise
[thr,sorh,keepapp]=ddencmp('den','wv',ecgsmooth);
cleanecg=wdencmp('gbl',C,L,'db4',8,thr,sorh,keepapp);

%Thresholding
max_value=max(cleanecg);
mean_value=mean(cleanecg);
threshold=(max_value-mean_value)/2;

%R detection
a5=appcoef(C,L,'db4',5);
C1=[a5;d5;d4;d3];
L1=[length(a5);length(d5);length(d4);length(d3);length(cleanecg)];
R_detect_signal=waverec(C1,L1,'db4');
R_detect_squared=R_detect_signal.^2;

%Beat Rate Extraction
for a=1:length(R_detect_squared)
  if R_detect_squared(a)>threshold
  R_detect_new(a)=R_detect_squared(a);
  else
  R_detect_new(a)=0;
  end
end
mean_R_detect=5*mean(R_detect_new);
for q=1:length( R_detect_new)-1
  if R_detect_new(q)< mean_R_detect
  R_detect_new(q)=0;
  end

end

d=0;
for b=1:length( R_detect_new)-1
  if ( R_detect_new(b)==0) & ( R_detect_new(b+1)~=0)
  d=d+1;
  indext(d)= b+1;
  end
end
fs_R_detect=length(R_detect_new)/20;
time=indext.*1/fs_R_detect;
ind=0;
for z=1:length(time)-1
  ind=ind+1;
  time_diff(ind)=time(z+1)-time(z);
end
av_time=mean(time_diff);
Square_Number=av_time/.2;
beat_Rate=300/Square_Number;
high=max(R_detect_new);

% Plot the Orginal Signal and Eliminating Baseline Drift signal
% Plot Denoised signal & R Detected Signal
subplot(4,1,1);plot(s1);title('ORIGINAL SIGNAL');
subplot(4,1,2);plot(s1-s2);title('BASELINE DRIFT ELIMINATION');
subplot(4,1,3);plot(cleanecg);title('MAIN SIGNAL(DENOISED ECG)');
subplot(4,1,4);plot(R_detect_new);title('R DETECTED SIGNAL');
text(length(R_detect_new)/2,high,['BEAT RATE = ',num2str(fix(beat_Rate))],'EdgeColor','Blue');