% QEXAMPLE Example application of Q transform related functions
%
% QEXAMPLE is a Matlab script which provides a simple example of how to use
% functions in the Q transform toolbox.  QEXAMPLE is useful as both a
% simple test of toolbox functions as well as a simple template for
% developing Q transform based applications.
%
% QEXAMPLE preforms the following steps.
%
%    1. tile the search space
%    2. read in example data
%    3. display raw time series data
%    4. condition the data
%    5. display conditioned time series data
%    6. apply the Q transform
%    7. display spectrogram of maximum Q plane
%    8. threshold on tile significance
%    9. display eventgram of significant tiles
%   10. exclude less significant overlapping tiles
%   11. display eventgram of non-overlapping tiles
%   12. write triggers to the standard output
%
% QEXAMPLE is provided with sample data that contains an optimally oriented 1.4,
% 1.4 solar mass inspiral waveform at a distance of 2.0 Mpc.  The peak signal
% energy occurs at a GPS time of approximately 822191381.5625.
%
% Execute "type qexample" from the Matlab plot to view the actual script.

% Shourov K. Chatterji
% shourov@ligo.caltech.edu

% $Id: qexample.m,v 1.7 2007/06/27 23:18:09 shourov Exp $

% tile signal space
disp('tiling signal space...');
reSampleFrequency = 4096;
timeRange = 32;
frequencyRange = [64 1024];
qRange = [4 64];
maximumMismatch = 0.2;
tiling = qtile(timeRange, qRange, frequencyRange, ...
               reSampleFrequency, maximumMismatch);

% load detector data
disp('loading detector data...');
frameCacheFile = '../framecaches/test.txt';
frameCache = loadframecache(frameCacheFile);
frameCache.directories = strcat('../', frameCache.directories);
frameType = 'TEST';
channelName = 'H1:LSC-DARM_ERR';
centerTime = 822191381.5625;
dataDuration = timeRange;
startTime = centerTime - dataDuration / 2;
stopTime = centerTime + dataDuration / 2;
timeShift = 0;
debugLevel = 1;
[data, sampleFrequency] = ...
    qreaddata(frameCache, channelName, frameType, ...
              startTime, stopTime, timeShift, debugLevel);

% resample data
disp('resampling data...');
data = qresample(data, sampleFrequency, reSampleFrequency);

% plot resampled time series
disp('plotting resampled data...');
plotReferenceTime = centerTime;
plotTimeRange = [-2 +2];
qtimeseries(data, tiling, startTime, plotReferenceTime, ...
            plotTimeRange, channelName);
drawnow;

% condition data
disp('conditioning data...');
calibration = [];
[data, calibration] = qcondition(data, tiling, calibration);

% plot conditioned time series
disp('plotting conditioned data...');
qtimeseries(qifft(data), tiling, startTime, plotReferenceTime, ...
            plotTimeRange, channelName);
drawnow;

% apply q transform
disp('transforming data...');
outlierFactor = 2.0;
transform = qtransform(data, tiling, calibration, outlierFactor);

% plot spectrogram
disp('plotting spectrogram...');
plotFrequencyRange = [];
plotQRange = max(qRange);
plotNormalizedEnergyRange = [];
qspectrogram(transform, tiling, startTime, plotReferenceTime, ...
             plotTimeRange, plotFrequencyRange, plotQRange, ...
             plotNormalizedEnergyRange, channelName);
drawnow;

% threshold on significance
disp('thresholding on significance...');
falseRate = 1e-2;
thresholdReferenceTime = [];
thresholdTimeRange = [];
thresholdFrequencyRange = [];
thresholdQRange = [];
maximumSignificants = 1e5;
analysisMode = [];
falseVetoRate =[];
uncertaintyFactor = [];
channelNames = [];
applyVeto = [];
debugLevel = 1;
significants = qthreshold(transform, tiling, startTime, falseRate, ...
                          thresholdReferenceTime, thresholdTimeRange, ...
                          thresholdFrequencyRange, thresholdQRange, ...
                          maximumSignificants, analysisMode, falseVetoRate, ...
                          uncertaintyFactor, applyVeto, channelNames, debugLevel);

% plot significants eventgram
disp('plotting significants...');
plotDurationInflation = 1.0;
plotBandwidthInflation = 1.0;
qeventgram(significants, tiling, startTime, plotReferenceTime, ...
           plotTimeRange, plotFrequencyRange, ...
           plotDurationInflation, plotBandwidthInflation, ...
           plotNormalizedEnergyRange, channelName);
drawnow;

% select mosaic events
disp('selecting non-overlapping tiles...');
durationInflation = 1.0;
bandwidthInflation = 1.0;
maximumMosaics = 1e3;
debugLevel = 1;
mosaics = qselect(significants, durationInflation, bandwidthInflation, ...
                  maximumMosaics, channelNames, debugLevel);

% plot mosaic eventgram
disp('plotting mosaics...');
qeventgram(mosaics, tiling, startTime, plotReferenceTime, ...
           plotTimeRange, plotFrequencyRange, ...
           plotDurationInflation, plotBandwidthInflation, ...
           plotNormalizedEnergyRange, channelName);
drawnow;

% write results
disp('writing results...');
fileName = 'example.txt';
triggerFields = {'time', 'frequency', 'duration', 'bandwidth', 'normalizedEnergy'};
qwritetriggers(mosaics, fileName, triggerFields, 'txt');

% done
disp('done.');

% return to calling function
return;