(***********************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which typically
will have the same name as this file except ending in ".nb" instead of
".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent to
using the Evaluate Initialiation Cells menu command in the front end.

DO NOT EDIT THIS FILE.  This entire file is regenerated automatically 
each time the parent Notebook file is saved in the Mathematica front end.
Any changes you make to this file will be overwritten.
***********************************************************************)

(* :Title: QRSdetection.m -- EKG QRS detection algorithm *)

(* :Context: Electrocardiography`QRSdetection`*)

(* :Author: Taras V. Pogorelov *)

(* :Summary:
  This package implements a QRS complex detection algorithm for the electrocardiograms,
  described in [1]. The package template ([2]) was used.
 *)

(* :Copyright: \[Copyright] 1998 by Wolfram Research Inc. *)

(* :Package Version: 1.0 *)

(* :Mathematica Version: 3.0 *)

(* :History: originally written by Taras V. Pogorelov, 1998 *)

(* :Keywords: EKG, QRS complex, detection, algorithm *)

(* :Sources:
   1. Willis J. Tompkins (Ed.). Biomedical Digital Signal 
   Processing, Prentice-Hall, NJ, 1993. (ISBN 0-13-067216-5)
   2. Roman E. Maeder. Programming in Mathematica, 3rd ed. 
   Addison-Wesley, 1996.(ISBN 0-201-85449-X)
*)

(* :Warnings: all functions of the package designed for EKG signals sampled with the standard sampling frequency of 200Hz.

*)

(* :Discussion:
  The package contains all the functions necessary for non-real time implementation of the QRS complex detection algorithm for an EKG signal sampled with 200Hz sampling frequency. This includes the following: lowpass (cutoff frequency ~11Hz) and highpass (cutoff frequency ~5Hz) digital IIR integer filters; digital differentiator; squaring function; digital moving window integrator (the window width is 30 samples, which corresponds to 150 msec); the complete detection function. For more details see [1].
*)

(* set up the package context, including public imports *)

BeginPackage["QRSdetection`"]

(* usage messages for the exported functions and the context itself *)

QRSdetection::usage = "QRSdetection.m is a package that contains all functions necessary for implementation of the QRS complex detection algorithm for electrocardiograms."

lowpassfilter::usage = "lowpassfilter[ekg] applies a digital lowpass filter to the sampled signal."

highpassfilter::usage = "highpassfiter[ekg] applies a digital highpass filter to the sampled signal."

derivative::usage = "derivative[ekg] differentiates the sampled signal."

squarelist::usage = "squarelist[ekg] squares elements of the list."

windowintegral::usage = "windointegral[ekg] applies a moving window integrator to the sampled signal."

QRSdetect::usage = "QRSdetect[ekg, opts] applies a QRS detection algorithm to the sampled signal, returns and displays the results."


(* error messages for the exported objects *)

lowpassfilter::list = "List of real numbers expected at position 1 in `1`"

highpassfilter::list = "List of real numbers expected at position 1 in `1`"

derivative::list = "List of real numbers expected at position 1 in `1`"

squarelist::list = "List of real numbers expected at position 1 in `1`"

windowintegral::list = "List of real numbers expected at position 1 in `1`"

QRSdetect::list = "List of real numbers expected at position 1 in `1`"


Begin["`Private`"]    (* begin the private context (implementation part) *)

Needs["Utilities`FilterOptions`"](*Import the FilterOptions` package*)


(* definition of the exported functions *)

realQ=(NumericQ[#]&&Im[#]\[Equal]0)&;

(*realQ[expr] gives True if expr is a numeric quantity and not a complex number, and False otherwise*)

lowpassfilter[ekg:{__?realQ}] := 
  Module[{padded, lowpass}, 
  		padded = Join[Table[0, {12}], ekg, Table[0, {12}]]; 
    	lowpass = Table[0, {Length[padded]}]; 
    	Do[
    		lowpass[[i]] = 
      		2*lowpass[[i - 1]] - lowpass[[i - 2]] + 
       		padded[[i]] - 2*padded[[i - 6]] + padded[[i - 12]], 
    	{i, 13, Length[padded]}
    	]; 
  N[1/36]*Drop[lowpass, 12]
  ]
  
lowpassfilter[]:=(Message[lowpassfilter::list, lowpassfilter];Return[$Failed])

lowpassfilter[_]:=(Message[lowpassfilter::list, lowpassfilter];Return[$Failed])

highpassfilter[ekg:{__?realQ}]:=
	Module[{padded,lowpass2,highpass},
		padded=Join[Table[0,{32}],ekg,Table[0,{32}]];
    	lowpass2=Table[0,{Length[padded]}];
		highpass=Table[0,{Length[padded]}];
		Do[ 			
			lowpass2[[i]]=lowpass2[[i-1]]+ padded[[i]]- padded[[i-32]];
			highpass[[i]]=padded[[i-16]]-N[1/32]*lowpass2[[i]],
	 	{i,33,Length[padded]}];
	Drop[highpass,32]
	]
	
highpassfilter[]:=(Message[highpassfilter::list, highpassfilter];Return[$Failed])
 
highpassfilter[_]:=(Message[highpassfilter::list, highpassfilter];Return[$Failed])
 
derivative[ekg:{__?realQ}]:=
	Module[{padded,deriv},
		padded=Join[Table[0,{4}],ekg,Table[0,{4}]];
		deriv=Table[0,{Length[padded]}];
		Do[ 
	deriv[[i]]=N[1/8]*(2*padded[[i]]+padded[[i-1]]-padded[[i-3]]-2*padded[[i-4]]),			{i,5,Length[padded]}
		];
	Drop[deriv,4]
	]
	
derivative[]:=(Message[derivative::list, derivative];Return[$Failed])

derivative[_]:=(Message[derivative::list, derivative];Return[$Failed])

squarelist[ekg:{__?realQ}]:=Map[#^2&,ekg]

squarelist[]:=(Message[squarelist::list, squarelist];Return[$Failed])

squarelist[_]:=(Message[squarelist::list, squarelist];Return[$Failed])

windowintegral[ekg:{__?realQ}]:=
	Module[{padded,integral,sum,j},
		padded=Join[Table[0,{29}],ekg,Table[0,{29}]];
    	integral=Table[0,{Length[padded]}];
		Do[ 
			integral[[i]]=N[1/29]*Sum[padded[[j]],{j,i-29,i}],
		{i,30,Length[padded]}
		];
	Drop[integral,29]
	]
	
windowintegral[]:=(Message[windowintegral::list, windowintegral];Return[$Failed])

windowintegral[_]:=(Message[windowintegral::list, windowintegral];Return[$Failed])

QRSdetect[ekg:{__?realQ}, opts___?OptionQ]:=
	Module[{integralekg,options},
	integralekg=windowintegral[squarelist[derivative[highpassfilter[ lowpassfilter[ekg]]]]];
	
	options=(PlotRange-> (PlotRange/.{FilterOptions[ListPlot,opts]}/.(PlotRange|Automatic|{})->{-35,110}));
	
	Show[GraphicsArray[ {
		ListPlot[ekg, options, FilterOptions[ListPlot,opts],
			PlotJoined\[Rule]True,GridLines->Automatic, Frame->True, 
			DisplayFunction-> Identity, 
			FrameLabel\[Rule]{"Sample number","Amplitude", FontForm["EKG Data 					Points",{"Times-Bold",10}],""}],
		ListPlot[integralekg, options, FilterOptions[ListPlot,opts],
			PlotJoined\[Rule]True, GridLines->Automatic, Frame->True, 
			DisplayFunction-> Identity,
			FrameLabel\[Rule]{"Sample number","Amplitude", FontForm["QRSdetect 					Result",{"Times-Bold",10}],""}] },
		ImageSize->{72 5.5, 72 2}
	]];
	
	integralekg]

(*options assigns the requested value to the PlotRange or keeps the default value of PlotRange->{-35,110}*)

QRSdetect[]:=(Message[QRSdetect::list, QRSdetect];Return[$Failed])
	
QRSdetect[_,x___]:=(Message[QRSdetect::list, QRSdetect];Return[$Failed])


End[ ]         (* end the private context *)

Protect[lowpassfilter, highpassfilter, derivative, squarelist, windowintegral, QRSdetect, realQ]    
(* protect exported symbols *)

EndPackage[ ]  (* end the package context *)