scripts/channelClasses/mock_ec/DualAlpha.java

/////////////////////////////////////////////////////////////////
// This script generates (somewhat) realistic EC data.
//
// Run with (for example)
//    java -cp build:lib/derby.jar -Dserver.cache.enable=false -Djava.library.path=c frontendClasses/CLI -scriptChannel mock_ec/DualAlpha -reviewSeries -paradigm ec -binZ site -transform power -display "TiledStack()" -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.  The EEG 
 * is a blend of:
 * <ol><li>a frontal source with a slow alpha</li>
 *     <li>a posterior source with a fast alpha</li>
 * </ol>
 *
 * <p>The intended effect is:
 * <ul><li>Two alpha sources, nearly overlapping in time, frequency 
 *         and topography</li>
 * </ul>
 */
public class DualAlpha 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 DualAlpha(Recording rec) {
        this.rec = rec;
    } // DualAlpha

    ////////////////////////////////////////////////////////////////////
    /** 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.090f),     // alpha1
            getEegWithAlpha(template,0.110f),     // alpha2
            getEcg(template),                     // ECG
        };

        // How recording sites are derived from sources
        // recording = factors*source
        //                   alpha1 alpha2  ECG blink pendular
        float[][] factors = {{0.0f, 1.0f, 0.0f},   // Fz
                             {0.4f, 0.4f, 0.0f},   // C3
                             {0.5f, 0.5f, 0.0f},   // Cz
                             {0.4f, 0.4f, 0.0f},   // C4
                             {0.9f, 0.1f, 0.0f},   // Pz
                             {0.0f, 0.0f, 1.0f}    // ECG
        };
 
        // Channels resulting from mixing sources: row labels of mixing matrix
        String[] labels = {"Fz", "C3", "Cz", "C4", "Pz", "ECG"};

        // 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);

    } // 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 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
}