scripts/channelClasses/standard/HRwqrs.java
/////////////////////////////////////////////////////////////////
// Processes ECG data by the WQRS algorithm, to generate an HR timeseries.
//
// This code is compiled at runtime, but compilation can be tested by:
// javac -cp build scripts/channelClasses/standard/HRwqrs.java
// rm scripts/channelClasses/standard/HRwqrs.class
//
// Try
// java -ea -cp build:lib/derby.jar:jri/JRI.jar -Dserver.cache.enable=false -Djava.library.path=c:jri frontendClasses/CLI -fn data/81237443-1.EC.NS5 -scriptChannel AvRef:HRwqrs -paradigm ec -reviewSeries
//
// and see, from the additional series, that the mean HR is about 59.8 bpm.
// Then try
//
// java -ea -cp build:lib/derby.jar:jri/JRI.jar -Dserver.cache.enable=false -Djava.library.path=c:jri frontendClasses/CLI -fn data/81237443-1.EC.NS5 -scriptChannel AvRef:HRwqrs -paradigm ec -scriptEpoch OneLong -reviewSeries -chanSelection "not Mode='EEG'" -transform power -binZ Site -score "TradSpect()" -display "TiledStack(offsetV=0.2, offsetH=0.25)"
//
// and see if ECG peaks at 0.970 Hz, and respiratory sinus arrhythmia is
// apparent in ECG_Rate at 0.2 -- 0.3 Hz.
/////////////////////////////////////////////////////////////////
package standard;
import java.io.*;
import java.util.*;
import generalClasses.*;
import recordingClasses.Recording;
import seriesClasses.*;
import channelClasses.ChannelScript;
import static channelClasses.Channel.*;
/////////////////////////////////////////////////////////////////
/** Processes ECG data by the WQRS algorithm, to generate an HR timeseries.
* The algorithm involves computing the 'length transform' of the the
* ECG, and applying an (adaptive) threshold to the transformed signal.
* It is completely independent of the Pan-Tompkins algorithm, which is
* used by the {@link scorerClasses.HeartRate} scoring algorithm.
* [It would be a simple matter to move the core of this algorithm into
* the <tt>generalClasses</tt> package, so that it could both be called
* from this script and from the <tt>HeartRate</tt> scoring option.]
*
* <p>The primary reference for the WQRS algorithm is<br>
* Gritzali, F; Frangakis, G; Papakonstantinou, G.
* Detection of the P and T waves in an ECG.
* <i>Computers and Biomedical Research</i>. 1989;22:83-91<br>
* PMID: 2914427<br>
* The code here is derived from <tt>wqrs.c</tt>, which is part of
* the PhysioToolkit at <a href="http://www.physionet.org">Physionet</a>.
*
* <p>The length transform at any time is simply the length of a plotted
* line segment, sqrt(deltaT^2 + deltaY^2). This transform has the
* effect of emphasising edges; it also results in a purely positive
* function. For both reasons it is a useful operation prior to
* threshold-based peak detection.
*/
public class HRwqrs extends ChannelScript
{
/** Recording instance to be operated on */
Recording rec = null;
/** Approximate duration of QRS complex, in seconds */
private float eyeClosing_s;
private static float eyeClosing_def = 0.35f;
/** Maximum QRS width, in seconds */
private float maxQRSw_s;
private static float maxQRSw_def = 0.13f;
/** If no QRS is detected over this period, the threshold is automatically
* reduced to a minimum value (tm_); the threshold is restored upon
* detection */
private float expectPeriod_s;
private static float expectPeriod_def = 2.5f;
/** Power frequency, in Hz */
private float pwFreq_Hz;
private static float pwFreq_def = 50.0f;
/** Minimum threshold */
private float tm_;
private static float tm_def = 2.0f;
/** Log of warnings generated during update */
private ArrayList<String> warnings = null;
////////////////////////////////////////////////////////////////////
/** Initialize instance by setting its parameters to specified values.
*/
public HRwqrs(Recording rec, float eyeClosing, float maxQRSw,
float expectPeriod, float pwFreq, float tm) {
this.rec = rec;
this.eyeClosing_s = eyeClosing;
this.maxQRSw_s = maxQRSw;
this.expectPeriod_s = expectPeriod;
this.pwFreq_Hz = pwFreq;
this.tm_ = tm;
warnings = new ArrayList<String>();
} // HRwqrs
////////////////////////////////////////////////////////////////////
/** Initialize instance by setting its parameters to default values.
*/
public HRwqrs(Recording rec) {
this(rec, eyeClosing_def, maxQRSw_def, expectPeriod_def,
pwFreq_def, tm_def);
} // HRwqrs
////////////////////////////////////////////////////////////////////
/** Compute the 'length transform' of the ECG signal.
* The 'length transform' is basically an edge-enhanced version of
* a time series, which in this case is the derivative of the ECG.
* This results in a simple triangular pulse at each beat, which
* makes threshold detection easy. The quantity <tt>lfsc</tt>, used
* in calculating the length transform, is mysterious: I reckon that
* 'dy' should be normalized and lfsc set to 1.0, and therefore have
* made it so.
* <p>The quantities <tt>v0,v1,v2</tt> are used to perform notch
* filtering at the mains frequency, while <tt>yn,yn1,yn2</tt>
* are involved in an autoregressive lowpass filter (with cutoff
* possibly at 32 Hz, according to one paper).
*
* @param series An ECG time series
* @return lbuf[i] contains the length-transformed signal, summed
* over the interval from i-ltWindow+1 to i. The values are
* dimensionless, and range from 1 to about 10.
*/
private float[] computeLT(SeriesAnalog series) {
float[] ecg = series.getY(); // get reference to time series
int nPoints = series.getNIndexes();
float sps = 1/series.getXDelta(); // sampling rate
//float gain = 0.01f; // uV/ADC_unit (approx)
//float lfsc = 1.25f*gain*gain/sps; // combined voltage and time resol
float lfsc = 1.0f; // alternative
float norm = 100000.0f; // characteristic slope of QRS waveform
// ebuf[] contains the LP derivative of the input, length transformed
float[] ebuf = new float[nPoints];
// lbuf[] contains a boxcar average of ebuf[]
float[] lbuf = new float[nPoints];
int ltWindow = Math.round(maxQRSw_s*sps); // boxcar length
float aet = ltWindow*lfsc; // latest boxcar sum
// Offsets used to perform notch filtering of the input
int lp1n = Math.round(sps/pwFreq_Hz);
float yn, yn1, yn2; // current and recent filtered input values
yn = yn1 = yn2 = 0;
for(int i=0; i<nPoints; i++) {
yn2 = yn1;
yn1 = yn;
if(i>=lp1n && i<nPoints-lp1n) {
float v0 = ecg[i-lp1n];
float v1 = ecg[i];
float v2 = ecg[i+lp1n];
yn = v0 - 2*v1 + v2 + 2*yn1 - yn2; // notch and LP filter
} else
yn = 0;
// Compute the normalized derivative of the filtered ECG signal
float dy = (yn-yn1)/series.getXDelta(); // differentiate => uV/s
dy /= norm; // normalize using some characteristic slope
// Compute length transform, and smooth
float et = (float)Math.sqrt(lfsc+dy*dy);
ebuf[i] = et;
aet += et - ((i>=ltWindow)? ebuf[i-ltWindow]: lfsc);
lbuf[i] = aet/ltWindow;
}
return lbuf;
} // computeLT
////////////////////////////////////////////////////////////////////
/** Update recording data by performing channel-oriented operations.
*/
public void update() {
// Find one and only ECG
ArrayList<SeriesAnalog> subset = selectMode(rec, DataMode.ECG);
if(subset.size()==0) {
warnings.add("HRwqrs: unable to find any ECG channels");
return;
} else if(subset.size()>1) {
String list = "";
for(SeriesAnalog s: subset) list += s.getPrimaryLabel()+" ";
warnings.add("HRwqrs: "+list+"all ECG; only one is allowed");
return;
}
SeriesAnalog series = subset.get(0);
int nPoints = series.getNIndexes();
float sps = 1/series.getXDelta(); // sampling rate
int eyeClosing = Math.round(eyeClosing_s*sps);
int expectPeriod = Math.round(expectPeriod_s*sps);
// Compute the length transform
float[] lt = computeLT(series);
// Initialize learning dur (t1) and thresholds (thresh0 and threshA)
int t1 = Math.round(8.0f*sps);
if(t1>nPoints*9/10) t1 = nPoints/2; // t1 > dur*0.9, so t1=dur/2
if(t1<expectPeriod) {
warnings.add("HRwqrs: duration of available data is too short");
return;
}
float thresh0 = 0;
for(int i=0; i<t1; i++) thresh0 += lt[i]/t1;
float threshA = 3*thresh0;
// Main loop
boolean learning = true;
float thresh1 = 0;
int timer = 0; // interval since the last detected QRS
Map<Integer,Float> beats = new TreeMap<Integer,Float>();
for(int i=0; i<nPoints; i++) {
// Set the current threshold, thresh1, to be 1 or 2*thresh0
if(learning) {
if(i<=t1)
thresh1 = 2*thresh0;
else { // Terminate learning and start over
thresh1 = thresh0;
learning = false;
i = 0;
}
}
// Compare length transformed value to current threshold
if(lt[i]>thresh1) {
// Found a possible QRS wavelet
// Check extrema in the region of this point
float max = lt[i];
float min = max;
for(int ii=i+1; ii<i+eyeClosing/2 && ii<nPoints; ii++)
if(lt[ii]>max) max = lt[ii];
for(int ii=i-1; ii<i-eyeClosing/2 && ii>=0; ii--)
if(lt[ii]<min) min = lt[ii];
// Are extrema sufficiently separated? [I am not sure
// what consideration motivates this extra criterion, so
// I set 'requiredGap' to something small. If lt[] is
// properly normalized, then the value should be O(1).]
float requiredGap = 0.10f;
if(max>min+requiredGap) {
// QRS detection is confirmed
// Now look for the onset of the QRS (PQ junction)
float onset = max/100;
int tpq = i-5;
for(int ii=i; ii>i-eyeClosing/2; ii--) {
if(lt[ii] - lt[ii-1] < onset &&
lt[ii-1] - lt[ii-2] < onset &&
lt[ii-2] - lt[ii-3] < onset &&
lt[ii-3] - lt[ii-4] < onset) {
tpq = ii;
break;
}
}
if(!learning) {
if(tpq<0 || tpq>=nPoints) break;
// Record tpq somehow
beats.put(i,tpq/sps);
timer = 0;
}
// Adjust thresholds
threshA += (max-threshA)/10;
thresh1 = threshA/3;
// Prevent further detections during the eye-closing period
i += eyeClosing;
}
} else if(!learning) {
timer++;
// If no longer learning, and no QRS has been detected
// recently, then decrease threshold towards tm_
if(timer>expectPeriod && threshA>tm_) {
// Adjust thresholds
threshA -= 1.0f;
thresh1 = threshA/3;
}
}
}
// Create an array of heart rates
float[] hr = new float[nPoints];
int last = 0;
float beatTime = 0;
boolean first = true;
for(Map.Entry<Integer,Float> entry: beats.entrySet()) {
int i = entry.getKey();
float rate = 60/(entry.getValue()-beatTime);
beatTime = entry.getValue();
if(first) // The first entry is unreliable, so should be skipped
first = false;
else
for(; last<i; last++) hr[last] = rate;
//System.out.println("Beat at "+beatTime+", bpm="+rate);
}
if(last>0)
for(; last<nPoints; last++) hr[last] = hr[last-1];
else {
for(; last<nPoints; last++) hr[last] = 0;
warnings.add("HRwqrs: no QRS events found");
}
// Build a HR series with args:
// SiteSet sites, float x0, float xDelta, Units xUnits,
// float[] y, Units yUnits, DataMode mode)
SiteSet siteLT=new SiteSet(new Site(series.getPrimaryLabel()+"_LT"));
SiteSet siteHR =new SiteSet(new Site(series.getPrimaryLabel()+"_Rate"));
float x0 = series.getFirstX();
float xDelta = series.getXDelta();
Units xUnits = series.getXUnits();
Units yUnits = new Units(Unit.bpm);
DataMode mode = DataMode.CED; // This avoids creating multiple ECG chans
SeriesAnalog ltSeries = new SeriesAnalog(siteLT,x0,xDelta,xUnits,lt,
new Units(Unit.none), mode);
SeriesAnalog hrSeries = new SeriesAnalog(siteHR,x0,xDelta,xUnits,hr,
yUnits, mode);
// Update Recording object
//appendToCurrentSeries(rec, ltSeries);
appendToCurrentSeries(rec, hrSeries);
} // update
////////////////////////////////////////////////////////////////////
/** Dump summary of this class or object
* @return String representation of this object
*/
public String toString() {
String s = "<<<"+this.getClass().toString()+">>>\n";
s += "WQRS algorithm assumes:\n";
s += "Approx QRS duration (s) = "+ eyeClosing_s+"\n";
s += "Max QRS duration (s) = "+ maxQRSw_s+"\n";
s += "Min threshold () = "+tm_+"\n";
s += "Mains freq (Hz) = "+ pwFreq_Hz+"\n";
for(String w: warnings)
s += w+"\n";
return s;
} // toString
}