1. Template functions and Signal Processing

The first example demonstrates the use of template functions and some of the signal processing capabilities of the DMT version of root. I start by generating a Time series containing two sinusoids and displaying the spectrum using the following commands.

    root[0] Sine sin3(3.0);
    root[1] TSeries ts(Time(0), Interval(0.005), 8192, sin3);
    root[2] ts += TSeries(Time(0), Interval(0.005), 8192, Sine(60.0));
    root[3] Spectrum(ts);

The first command defines a 3.0 Hz sine function template, "sin3". The second command
creates a time series (TSeries class) object called "ts" and fills it with the data from sin3. The first three argument of the TSeries constructor are the absolute (GPS) start time, the sample interval (1/sample-rate) and the number of samples. Next, an equal length 60 Hz sinusoid is added to "ts". The right hand side (rhs) of this statement is equivalent to the first two statements in that it generates a temporary TSeries from a 60 Hz sinusoidal template. The "+=" operator is overloaded in the TSeries class definition to add the rhs TSeries to the lhs series on an element by element basis. Finally, the spectrum is displayed using the Spectrum(FSeries& fs) function. Note that "ts" is converted automatically to a frequency series (FSeries), including the implied Fourier transform, by an FSeries class constructor. The resulting spectrum is shown below:

Note that the two frequency components are spread out by the Fourier transform. This behavior is to be expected in an unwindowed transform of a signal that doesn't meet the periodic boundary conditions. This can be remedied by windowing the TSeries data with e.g. a Hanning window. The following statements demonstrate how this may be done.

    root[4] Hanning Wn;
    root[5] TSeries tw = Wn(ts);
    root[6] Spectrum(tw);

The first statement defines a Hanning window object. The windows are represented as classes derived from the Window base class to make the similarity between window classes apparent as well as to allow the sharing of the windowing code. In the second statement, I create a new windowed TSeries, "tw". This is again displayed using Spectrum(), giving the plot shown below.

The resulting spectrum still shows prominent signals at 3 Hz and 60 Hz. The signal diffusion due to the aperiodicity is reduced to the level of precision of the Fourier transform (this produces the jagged spectrum below the 10-15 level).

2. Filters

The DMT Sandbox contains a Filter class that allows the definition of filters which may then be used to filter signals in the time domain. At present (22-July-1999) the filter class supports only FIR filters, but additional (IIR) filter support should be possible in the future. The first step is to design a filter. As an example, a 10 Hz Low-Pass filter is constructed using the following commands:

    root[7] Filter f(25, 200.0);
    root[8] f.dFirW(25, "Hamming", "LowPass", 10.0);
    root[9] Bode(f.Xfer(0.1,100.0,0.1));

The first statement defines a filter object. The two arguments give the order (number of zeros) of the filter and the sampling frequency of the input stream. The second statement invokes an FIR design method to set the filter coefficients for a 10.0 Hz low-pass filter. The dFirW method uses a window-based design technique. In this case a Hamming window is used to design the low-pass filter. The following statement calculates the transfer function of the filter and display it as a Bode plot. The arguments to the Xfer method specify the frequency range (0.1-100 Hz) and the frequency increment (0.1 Hz), respectively. The resulting Bode plot is shown below.

This filter can then be applied to a time series (for example the time series "ts" creates in the example above) with the following commands.

    root[10] TSeries tf = f(ts);
    root[11] Spectrum(tf);

The spectrum of the resulting time-series is shown below


3. Spectral Power Strip Chart

The second example plots the spectral power in 6 different frequency bands of a single channel (in this case the 40m IFO_DMRO signal) as a function of time. This example uses
many of the DMT classes to simplify the coding. The C++ function used for this display
is listed below:

//>>>>>>>>>>      File:  $DMTHOME/root/macros/tscart.cc
//    Make a constantly updating spectrum plot
//    dTime   contains the time bin to be plotted
//    nbinF   contains the number of frequency bands
//    binf[]  contains the frequency band limits.
void tscart(void) {
  //---------------------------------  parameters.
    double dTime = 1.0;
    float binf[7]={1.00, 3.00, 10.00, 30.00, 100.00, 300.0, 1000.0};
    int   nbinF = 6;
    int   NStep = 1000;

  //---------------------------------  define constants.
    float pi = 3.14159265358;
    Interval FillTime(dTime);
    double SampleRate = 16384.0;
    double tsecs = dTime;

  //---------------------------------  Create a TSeries.
    int nbin = int(SampleRate*tsecs);
    TSeries* ts = new TSeries(Time(0), Interval(0.0), nbin);
    In.addChannel("IFO_DMRO", 0, &ts);    // Ask for data, say where it goes

  //---------------------------------  Create the strip chart.
    StripChart sc(nbinF, "Seismic bands");
    char line[64];
    for (int p=0 ; p<nbinF ; p++) {
        sprintf(line, "%.0f-%.0f Hz", binf[p], binf[p+1]);
        sc.setTitle(p, line);
    sc.setAxisType(1);    // x-axis interpreted as a date/time
    bool first = true;

  //---------------------------------  Loop over sampling intervals
    FSeries fs;
    while (1) {
        *ts += 0.0;            // subtract out DC offset
        *ts *= 1.0;            // Default (non-)calibration
        fs.setData(*ts);         // FFT data into FSeries

      //------------------------------  Get average values from FFT
        float df  = fs.getFStep();       //  Frequency step size
        unsigned long nwd = fs.getNStep(); //  Number of frequency bins
        float* sum = new float[nbinF];
        for (int i=0 ; i<nbinF ; i++) {
            float lf = binf[i];
            float hf = binf[i+1];
            sum[i] = fs.Power(lf, hf);  // get the power, divide by Hz
            sum[i] = sqrt(sum[i])*df/((hf-lf)*pi*(lf+hf));
            if (first) sc.setStripe(i, sum[i], 10.0); // center pens on the first sample
        Time t = ts.getStartTime(); // Use Unix date representation for
        double tbin = getUnix(t);   //     the x-axis
        if (first) {
            sc.setMaxX(tbin + NStep*dTime);
            first = false;
        sc.plot(tbin, sum);  // Plot the data to the strip chart
        delete[] sum;

The macro displayed above can be run under root as follows:

Note the reuse of the (automatically created) Dacc object "In".  The close() method
disconnects root from the online data partition, and the addFile() method reconnects it to the specified file(s). addFile() may be called multiple times to line up several frame files (or sets of files) for screening.

The picture below was generated by the root session listed above. The data used in
this test run were stored in frame files, and were recorded during the 1994 running of the 40m interferometer.