scripts/epochClasses/ec/Beats.java
/////////////////////////////////////////////////////////////////
// Define epochs time-locked to heart beats
/////////////////////////////////////////////////////////////////
package ec;
import java.io.*;
import java.text.*;
import java.util.*;
import epochClasses.*;
import generalClasses.*;
import recordingClasses.Recording;
import seriesClasses.*;
import static epochClasses.EpochSelector.*;
/////////////////////////////////////////////////////////////////
/** Define epochs time-locked to heart beats
* This script generates a list of epoch times; and for each epoch
* generates a standard set of attributes. The attribute table is
* a flexible way for attaching all manner of information to each epoch.
* This information may be used later when subaveraging epochs.
*
* <p>This script also defines the start and duration of epochs through
* <ul><li><tt>defaultEpochStart</tt>: the offset of epoch beginnings from
* the times listed in the the (obligatory) attribute <tt>Time</tt>.
* A value of 0.0 is common for EEGs.</li>
* <li><tt>defaultEpochDur</tt>: the duration of epochs, in seconds.
* A value of 2.0 is common for ERPs.</li>
* <li><tt>dt</tt>: the increment between successsive epochs, in sec</li>
* </ul>
*
* <p>Attribute values, and epoch start and duration, are reported
* to the calling program via the methods defined in the superclass
* <tt>EpochScript</tt>
*/
public class Beats extends EpochScript
{
private static float defaultEpochStart = -1.0f;
private static float defaultEpochDur = 2.0f; // cf loHzFit if fitting
//private static float dt = 5.0f; // step between successive epochs, sec
private float[] times = null; // t=0 point of epochs, in lieu of evs
private String[][] th = null; // column header information
private Object[][] td = null; // main results array
private int iteratorEventN = 0; // used only by iterator, next()
private int nAttrib; // number of attribute columns
private ArrayList<Event> runtimeEvents = null; // obtained from Recording
////////////////////////////////////////////////////////////////////
/** This script generates a standard set of attributes:
* <ul><li><tt>Time</tt>: points within the span of the recording. In
* the case of ERP recordings, these times usually coincide with
* significant stimuli. However they can be offset for stimuli, or
* else identify EDA onsets, K-complexes, hypoxic events, etc. <b>This
* attribute is obligatory.</b></li>
* <li><tt>Stim</tt>: a string label that is commonly used for
* identifying epochs. It is here that an event (which contains
* information like Tone, 1000 Hz) is assigned the label "Bg".</li>
* <li><tt>AlphaPow, AlphaPha, AlphaQ</tt>: values that quantify
* the magnitude and phase of alpha.</li>
* <li><tt>ScalpSD, ScalpPP</tt>: maximal SD and peak-to-peak range
* within each epoch. Only scalp channels are considered.</li>
* <li>And so on …</li>
* </ul>
* (Check code for a definitive list.) It also sets the epoch start
* offset and epoch duration.
* @param rec Reference to the data. This is used to retrieve both time
* series and event information.
*/
public Beats(Recording rec) {
super(rec, defaultEpochStart, defaultEpochDur);
float allTimes[] = null;
ArrayList ss = getSeries();
for(int chanN=0; chanN<ss.size(); chanN++) {
if(ss.get(chanN) instanceof SeriesAnalog) {
SeriesAnalog s = (SeriesAnalog)(ss.get(chanN));
if(s.getMode()==DataMode.ECG) {
allTimes = QRSpt.getTimes(s,null);
break;
}
}
}
if(allTimes==null) {
CommonLog.get().severe("In ec/Standard.java: no ECG channel found");
System.exit(2);
}
float t0 = 0;
float t1 = 0;
try {
t0 = rec.getX0();
t1 = rec.getX1();
} catch(recordingClasses.RecordingException e) {
CommonLog.get().severe("In ec/Standard.java: "+e.toString());
System.exit(2);
}
// Count the number of QRS events for which the associated epoch
// is wholely contained within the duration of the recording.
// (It is likely that the first and last QRS will need to be
// dropped, as the epoch would extend beyond the available range.)
int nEvents = 0;
for(int i=0; i<allTimes.length; i++) {
float thisT = allTimes[i];
if(thisT+defaultEpochStart>t0 &&
thisT+defaultEpochStart+defaultEpochDur<t1) nEvents++;
}
// Build new list of times
times = new float[nEvents];
int j = 0;
for(int i=0; i<allTimes.length; i++) {
float thisT = allTimes[i];
if(thisT+defaultEpochStart>t0 &&
thisT+defaultEpochStart+defaultEpochDur<t1) times[j++] = thisT;
}
// Generate meta information about attribute table
String[] names = {
"Time", "Stim",
//"AlphaPow", "AlphaPha", "AlphaQ",
"AlphaPow", "AlphaBig",
"ScalpSD", "ScalpPP"};
String[] types = {
"FLOAT", "VARCHAR6",
//"FLOAT", "FLOAT", "VARCHAR3",
"FLOAT", "VARCHAR14",
"FLOAT", "FLOAT"};
assert names.length==types.length: "Standard.java: length mismatch";
nAttrib = names.length;
th = new String[2][nAttrib];
for(int i=0; i<nAttrib; i++) {
th[0][i] = names[i];
th[1][i] = types[i];
}
// Allocate attribute table
assert (nEvents>=1): "Logic error in ec/Standard: nEvents="+nEvents;
td = new Object[nEvents][nAttrib];
// Initialize attribute table, one (or a few) columns at a time
updateTime(0);
updateStim(1);
//updateSpectralPowerAtStim(2,3,4);
updateSpectralPowerInEpoch(2, 3);
updateArtif(4,5);
} // Standard
////////////////////////////////////////////////////////////////////
/** Dump summary of this object, to be used for diagnostics
* @return String representation of this object
*/
@Override
public String toString() {
String s = super.toString(); // get superclass description
return s;
} // toString
////////////////////////////////////////////////////////////////////
/** This returns the name and datatype of any number of attributes.
* Each name and datatype characterizes one column of the table
* of attributes.
* @return The name and datatype for each of <tt>nCols</tt> columns
* of the attribute table. The array is dimensioned [2][nCol].
*/
@Override
public String[][] getAttribNameType() {
return th;
} // getAttribNameType
////////////////////////////////////////////////////////////////////
/** The attribute table is returned one row at a time. This checks
* if there are more rows available.
*/
@Override
public boolean hasNext() {
return td!=null && iteratorEventN < td.length;
} // hasNext
////////////////////////////////////////////////////////////////////
/** Returns the next row of attribute table. If there are no more
* rows remaining, then <tt>null</tt> is returned. Note that is
* possible that each element in the returned array is itself
* <tt>null</tt>, so calling code should be prepared to deal with
* such cases too.
*/
@Override
public Object[] next() {
return (td!=null && iteratorEventN<td.length)?
td[iteratorEventN++]: null;
} // next
////////////////////////////////////////////////////////////////////
// private methods //
////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////
/** Initializes td[*][colTime] with the event times, in seconds.
* NOTE: this is obligatory, unlike all other array columns.
*/
private void updateTime(int colTime) {
for(int row=0; row<td.length; row++) {
td[row][colTime] = new Float(times[row]);
}
} // updateTime
////////////////////////////////////////////////////////////////////
/** Initializes td[*][colStim].
*
* This version is specific to EC/EO paradigms where the 'stimuli'
* are pseudo.
*/
private void updateStim(int colStim) {
// Compute average R-R time
float avRR = 0;
for(int i=0; i<times.length-1; i++) avRR += times[i+1]-times[i];
avRR /= (times.length-1);
for(int row=0; row<td.length; row++) {
float rr = 0;
if(row<td.length-1) rr = times[row+1]-times[row];
else rr =times[td.length-1]-times[td.length-2];
td[row][colStim] = (rr>avRR)? "SlowRR": "FastRR";
}
} // updateStim
////////////////////////////////////////////////////////////////////
/** Initializes td[*][colPow] td[*][colPha] and td[*][colPhaseBand],
* which refer to spectral band power, phase, and a categorical version
* of phase. Alpha is quantified for some brief (e.g. 1-sec) period
* following the times listed in <tt>float[] times</tt>.
*
* This version is specific to 'Pz' and to the alpha band.
*/
private void updateSpectralPowerAtStim(int colPow, int colPha,
int colPhaseBand) {
// Which site to look at
String site = "Pz";
// Definition of spectral band: binWidth * binN = 10.0 +/- binWidth/2
float binWidth = 1.0f;
int binN = 10;
// Calculate alpha power and phase at Pz
ArrayList ss = getSeries();
int chanN = 0;
for(chanN=0; chanN<ss.size(); chanN++) {
Series s = (Series)(ss.get(chanN));
if(s.getPrimaryLabel().equalsIgnoreCase(site)) break;
}
Float[][] alpha = new Float[td.length][2];
if(chanN<ss.size())
getEegAtChan((Series)(ss.get(chanN)), binWidth, binN, times, alpha);
// Update power, phase and categorized phase
for(int row=0; row<td.length; row++) {
td[row][colPow] = alpha[row][0]; // power
td[row][colPha] = alpha[row][1]; // phase
if(alpha[row][1] != null) {
int cat = Math.round(alpha[row][1]/180); if(cat<0) cat+=2;
td[row][colPhaseBand] = (cat==1)? "Out": "In";
} else
td[row][colPhaseBand] = null;
}
} // updateSpectralPowerAtStim
////////////////////////////////////////////////////////////////////
/** Initializes td[*][colPow] and td[*][colPowBig],
* which refer to spectral band power, and a categorical version
* of power. Alpha is quantified for the period of duration
* <tt>getEpochDur()</tt> following the times listed in
* <tt>float[] times</tt>.
*
* This version is specific to 'Pz' and to the alpha band.
*/
private void updateSpectralPowerInEpoch(int colPow, int colPowBig) {
// Which site to look at
String site = "Pz";
// Definition of spectral band: bandMin ... bandMax
float bandMin = 8.0f; // Hz, in range 0...Nyquist
float bandMax = 13.0f; // Hz, in range 0...Nyquist
float fDelta = 1.0f; // Something < bandMax-bandMin
// Threshold between 'small' and 'big' alpha powers
float thresh; // uV^2 (computed)
// Calculate alpha power at Pz
ArrayList ss = getSeries();
int chanN = 0;
for(chanN=0; chanN<ss.size(); chanN++) {
Series s = (Series)(ss.get(chanN));
if(s.getPrimaryLabel().equalsIgnoreCase(site)) break;
}
Float[] pow = new Float[td.length];
float[] pow2 = new float[td.length];
if(chanN<ss.size()) {
SeriesAnalog s = (SeriesAnalog)(ss.get(chanN));
for(int row=0; row<td.length; row++) {
// No need to check for exceptions from segment()...
float tStart = times[row]+getEpochStart();
SeriesAnalog seg = s.segment(tStart,getEpochDur(),
getEpochStart());
pow[row] = seg.fourierPower(fDelta).integral(bandMin,bandMax);
pow2[row] = pow[row];
}
// Find the median power
Arrays.sort(pow2);
if(td.length % 2 == 0)
thresh = (pow2[td.length/2-1]+pow2[td.length/2])/2;
else
thresh = pow2[(td.length-1)/2];
String th = ""+((int)(thresh))+" uV^2";
// Update
for(int row=0; row<td.length; row++) {
td[row][colPow] = pow[row];
td[row][colPowBig] = (pow[row].floatValue()<thresh)?
"<"+th: ">"+th;
}
}
} // updateSpectralPowerInEpoch
////////////////////////////////////////////////////////////////////
/** Initializes td[*][colScalpSd] and td[*][colScalpPP], which refer
* to the maximal standard deviation and peak-to-peak values, considering
* all scalp signals. These values are suitable for threshold-based
* epoch rejection.
*/
private void updateArtif(int colScalpSd, int colScalpPP) {
ArrayList ss = getSeries();
for(int row=0; row<td.length; row++) {
// Calculate SD and P-P measures by finding max in EEG chans only
float maxSd = 0;
float maxPP = 0;
float tStart = times[row]+getEpochStart();
for(int chanN=0; chanN<ss.size(); chanN++) {
if(ss.get(chanN) instanceof SeriesAnalog) {
SeriesAnalog s = (SeriesAnalog)(ss.get(chanN));
if(s.getMode()==DataMode.EEG &&
s.indexOf(tStart) >= 0 &&
s.indexOf(tStart+getEpochDur()) <= s.getNIndexes()) {
BasicStats bs = s.segment(tStart,getEpochDur(),
getEpochStart())
.getStats();
if(bs.sd > maxSd) maxSd = bs.sd;
if(bs.max-bs.min > maxPP) maxPP = bs.max-bs.min;
}
}
}
// Update SD and P-P measures
td[row][colScalpSd] = maxSd;
td[row][colScalpPP] = maxPP;
}
} // updateArtif
}