scripts/channelClasses/mock_ec/Grat.java
/////////////////////////////////////////////////////////////////
// This script generates (somewhat) realistic EC data.
//
// Run with (for example)
// setenv CLASSPATH build:lib/nr.jar:lib/derby.jar
// java -cp build:lib/derby.jar -Dserver.cache.enable=false -Djava.library.path=c frontendClasses/CLI -scriptChannel mock_ec/Grat -reviewSeries -paradigm ec -binZ site -v
/////////////////////////////////////////////////////////////////
package mock_ec;
import java.io.*;
import java.util.*;
import epochClasses.*;
import generalClasses.*;
import recordingClasses.Recording;
import seriesClasses.*;
import seriesClasses.seriesGeneration.Erlang;
import channelClasses.ChannelScript;
import static channelClasses.Channel.*;
/////////////////////////////////////////////////////////////////
/** This script generates (somewhat) realistic EC data.
* This script tests the Gratton algorithm
* <ol><li>EEG channels, with realistic spectra</li>
* <li>the EOGv channel, in which there are several blinks</li>
* <li>the EOGh channel, which contains only noise with f<2 Hz</li>
* <li>ECG channel, which should be ignored by Gratton</li>
* </ol>
*
* <p>The intended effect is:
* <ol><li>vertical blink factors and vertical saccade factors should be about
* equal, and equal to the second last column of factors[][]</li>
* <li>horizontal blink factors and horizontal saccade factors should be
* about equal, and equal to the last column of factors[][]</li>
* </ol>
* Except it wasn't, at first. Horizontal blink factors were random numbers
* in the range (-1,+1), approximately. This is because the horizontal source
* contributes just 20 to 30% of the variance of any of the time series;
* and because there are only 0.5% of points clasified as belonging to
* blinks. The horizontal blink factors can be tidied up by increasing
* the size of the horizontal noise source (10x), or increasing the number of
* blinks. An example of the expected match is:<pre>
* Site ------ Blinks ------ ---- Saccacdes -----
* Vertical Horizontal Vertical Horizontal
* Fz 0.46 0.40 0.61 0.32
* Cz 0.39 0.37 0.51 0.32
* Pz 0.39 0.30 0.39 0.30
* EOGv 1.00 0.00 1.00 0.00
* EOGh 0.00 1.00 -0.00 1.00
* </pre>
* which also illustrates how the blink estimates are never as good as
* those for saccades.
*/
public class Grat extends ChannelScript
{
/** Recording instance to be operated on */
Recording rec = null;
/** Time series */
ArrayList<SeriesAnalog> list = new ArrayList<SeriesAnalog>();
/** Events. May be left empty */
ArrayList<Event> ev = new ArrayList<Event>();
////////////////////////////////////////////////////////////////////
/** Initialize instance by setting its parameters to default values.
*/
public Grat(Recording rec) {
this.rec = rec;
} // Grat
////////////////////////////////////////////////////////////////////
/** Update recording data by performing channel-oriented operations.
*/
public void update() {
// Template - used to encapsulate all sampling characteristics
float x0 = 10.0f; // in seconds: times start at x0
float xDelta = 0.004f; // in seconds: times increment by xDelta
float duration = 120.0f; // in seconds: times end at x0+duration
int nIndexes = Math.round(duration/xDelta);
SeriesAnalog template = new SeriesAnalog(new SiteSet(), // sites
x0, // x0
xDelta, // xDelta
new Units(Unit.s), // xUnits
nIndexes, // # samples
new Units(Unit.uV), // yUnits
DataMode.EEG); // DataMode
// Source signals, corresponding to columns of mixing matrix
SeriesAnalog[] source = {
getEegWithAlpha(template,0.085f), // alpha1
getEegWithAlpha(template,0.095f), // alpha2
getEcg(template), // ECG
getEogV(template), // blink
getEogH(template) // pendular
};
// How recording sites are derived from sources
// recording = factors*source
// alpha1 alpha2 ECG blink pendular
float[][] factors = {{0.2f, 1.0f, 0.0f, 0.6f, 0.3f}, // Fz
{0.5f, 0.8f, 0.0f, 0.5f, 0.3f}, // Cz
{0.9f, 0.2f, 0.0f, 0.4f, 0.3f}, // Pz
{0.0f, 0.0f, 1.0f, 0.0f, 0.0f}, // ECG
{0.0f, 0.0f, 0.0f, 1.0f, 0.0f}, // EOGv
{0.0f, 0.0f, 0.0f, 0.0f, 1.0f} // EOGh
};
// Channels resulting from mixing sources: row labels of mixing matrix
String[] labels = {"Fz", "Cz", "Pz", "ECG", "EOGv", "EOGh"};
// Generate time series from source[] and factors[][]
for(int site=0; site<labels.length; site++) {
// What modality?
DataMode mode = DataMode.EEG;
if(labels[site].matches("[eE][oO][gG].*")) mode=DataMode.EOG;
else if(labels[site].matches("[eE][cCkK][gG].*"))mode=DataMode.ECG;
// Create time series, with general character matching template
SeriesAnalog sum = template.clone();
// Form linear combination of sources
for(int s=0; s<factors[0].length; s++)
sum.add( source[s].clone().mul(factors[site][s]) );
sum.setSites(new SiteSet(new Site(labels[site])));
sum.setMode(mode);
// Append sum to result
list.add(sum);
}
// Add synthetic time series to the currently empty Recording
replaceAllSeries(rec, list);
replaceAllEvents(rec, ev);
doGratton(rec, 0);
} // update
////////////////////////////////////////////////////////////////////
/** Dump summary of this class or object
* @return String representation of this object
*/
public String toString() {
String s = "<<<"+this.getClass().toString()+">>>\n";
return s;
} // toString
////////////////////////////////////////////////////////////////////
/** Generate pseudo EEG containing alpha
*/
private SeriesAnalog getEegWithAlpha(SeriesAnalog template, float t0) {
float lnorm = 2.8f;
return SeriesAnalog.getModelledEegTimeseries(template,t0,lnorm);
} // getEegWithAlpha
////////////////////////////////////////////////////////////////////
/** Generate pseudo EOGv containing 'blinks' and low frequency noise
* If the blink duration is blinkDur=0.1 s and impulses are scaled by
* blinkArea=10 uV.s, then we should see a clear 100 uV peak in EOG, and
* proportionate peaks in EEG. The corresponding variances should be
* 2/9 amp^2.
*/
SeriesAnalog getEogV(SeriesAnalog template) {
float blinkDur = 0.1f; // width of blink, in sec
float blinkArea = 10.0f; // width*height of blink, in uV sec
float eogLP = 2.0f; // LP cutoff for EOG signals
float eogRms = 5.0f; // rms of ongoing EOG signals
float x0 =template.getFirstX(); // Initial time
float xDelta = template.getXDelta(); // sampling interval, in seconds
float white = eogRms/(float)Math.sqrt(eogLP*xDelta*2); // equiv noise
float duration = template.getXDelta()*template.getNIndexes();
// Blink times
float[] timesBl = {
x0+0.1003f*duration, x0+0.3015f*duration,x0+0.7027f*duration,
x0+0.75f*duration, x0+0.77f*duration, x0+0.79f*duration,
x0+0.81f*duration, x0+0.83f*duration, x0+0.85f*duration,
x0+0.87f*duration, x0+0.89f*duration, x0+0.91f*duration,
x0+0.93f*duration, x0+0.95f*duration, x0+0.97f*duration};
return SeriesAnalog.getImpulseTimeseries(template,timesBl)
.filterBox(blinkDur).mul(blinkArea)
.add(SeriesAnalog.getWhiteTimeseries(template,white).filterLP(eogLP));
} // getEogV
////////////////////////////////////////////////////////////////////
/** Generate pseudo EOGh containing low frequency noise only
*/
SeriesAnalog getEogH(SeriesAnalog template) {
float eogLP = 2.0f; // LP cutoff for EOG signals
float eogRms = 5.0f; // rms of ongoing EOG signals
float xDelta = template.getXDelta(); // sampling interval, in seconds
float white = eogRms/(float)Math.sqrt(eogLP*xDelta*2); // equiv noise
return SeriesAnalog.getWhiteTimeseries(template,white).filterLP(eogLP);
} // getEogH
////////////////////////////////////////////////////////////////////
/** Generate pseudo ECG.
* The wavelet is a standard ERP, but deliberately truncated to
* make it more spike-like.
* The ECG will have R-R intervals of 0.80+var seconds, where 'var'
* has a gamma distribution, with mean equal to 0.1 seconds.
*/
private SeriesAnalog getEcg(SeriesAnalog template) {
float ecgDuration = 0.2f;
float xDelta = template.getXDelta(); // sampling interval, in seconds
int nIndexes = Math.round(ecgDuration/xDelta);
if((nIndexes&1)==0) nIndexes++;
float x0 = -(nIndexes-1)*xDelta/2; // so wavelet is centred on 0 secs
SeriesAnalog wavelet = new SeriesAnalog(new SiteSet(), // sites
x0, // x0
xDelta, // xDelta
new Units(Unit.s), // xUnits
nIndexes, // # samples
new Units(), // yUnits
DataMode.ECG); // DataMode
float lnorm = -5.0f;
SeriesAnalog erp =SeriesAnalog.getModelledErpTimeseries(wavelet,lnorm);
Erlang ran = new Erlang(0.1/5, 5, 0); // scale = mean/5
ArrayList<Float> timesList = new ArrayList<Float>();
float t = template.getFirstX()+0.3f; // time of first ECG event
float duration = template.getXDelta()*template.getNIndexes();
while(t<duration) {
timesList.add(new Float(t));
t += (float)(0.8+ran.gen());
}
float[] times = new float[timesList.size()];
for(int i=0; i<timesList.size(); i++) times[i] = timesList.get(i);
return SeriesAnalog.getImpulseTimeseries(template,times)
.mul(xDelta).convolveAsym(erp);
} // getEcg
}