% TESTRECOVERY Test recovery of signal to noise ratio
%
% TESTRECOVERY tests the accuracy with which the Q transform recovers
% the signal to noise ratio of sinusoidal Gaussian signals.
%
% TESTRECOVERY repeatedly applies the Q transform to sinusoidal
% Gaussian signals injected with random center time, center frequency,
% Q, phase, and signal to noise ratio into stationary white noise.
% For every injection, TESTRECOVERY records the measured signal to
% noise ratio in the template closest to the injected signal
% parameters.
%
% As it runs, TESTRECOVERY writes the signal center time, center
% frequency, Q, injected SNR, and measured SNR a log file
% testrecovery/testrecovery.log.  TESTRECOVERY can also produce
% intermediate scatter plots of measured vs. injected SNR after each
% injection if the hardcoded watchFlag is set to a non-zero value.
% Upon completion, TESTRECOVERY also writes results to a MAT file
% testrecovery/testrecovery.mat.
%
% Upon completion, TESTRECOVERY also produces the following plot.
%
%   testrecovery/testrecovery_1.png
%     Scatter plot of measured vs. injected signal to noise ratio.
%
%   testrecovery/testrecovery_2.png
%     Histogram of measured vs. injected signal to noise ratio.
%
% See also QTILE and QTRANSFORM.

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

% $Id: testrecovery.m,v 1.4 2007/04/11 14:36:38 shourov Exp $

% path to qtransform code base
path('..', path);

% create output directory
unix('mkdir -p testrecovery');

% open log file
logFileName = 'testrecovery/testrecovery.log';
logFileID = fopen(logFileName, 'w');

% number of trials
numberOfTrials = 1000;

% boolean flag to produce intermediate plots
watchFlag = 0;

% signal space parameters
timeRange = 64;
qRange = [sqrt(11) 100];
frequencyRange = [0 Inf];
sampleFrequency = 4096;

% desired tiling accuracy
maximumMismatch = 0.2;

% tile the signal space
tiling = qtile(timeRange, qRange, frequencyRange, sampleFrequency, ...
               maximumMismatch);

% time vector for signal construction
time = 0 : 1 / sampleFrequency : timeRange - 1 / sampleFrequency;

% initialize random number generator
rand('state', sum(100 * clock));
randn('state', sum(100 * clock));

% initialize results vectors
centerTimes = [];
centerFrequencies = [];
qs = [];
phases = [];
injectedSNRs = [];
measuredSNRs = [];

% begin loop over trials
for trialNumber = 1 : numberOfTrials,

  % select random center time
  centerTime = unifrnd(0 + tiling.transientDuration, ...
                       tiling.duration - tiling.transientDuration, ...
                       1, 1);;

  % select random quality factor
  q = 2.^unifrnd(log2(tiling.minimumQ), ...
                 log2(tiling.maximumQ));

  % select random center frequency
  minimumFrequency = 50 * q / (2 * pi * tiling.duration);
  maximumFrequency = tiling.sampleFrequency / (2 * (1 + sqrt(11) / q));
  centerFrequency = 2.^unifrnd(log2(max(minimumFrequency, ...
                                        tiling.highPassCutoff)), ...
                               log2(min(maximumFrequency, ...
                                        tiling.lowPassCutoff)));

  % select random phase
  phase = unifrnd(0, 2 * pi, 1, 1);

  % select random SNR
  snr = 10.^unifrnd(log10(sqrt(2)), log10(100));

  % construct noise
  noise = normrnd(0, sqrt(tiling.sampleFrequency / 2), 1, length(time));

  % construct signal
  signal = exp(-(time - centerTime).^2 * ...
               4 * pi^2 * centerFrequency^2 / q^2) .* ...
           sin(2 * pi * centerFrequency * (time - centerTime) + phase);

  % scale signal to desired snr
  signal = sqrt(snr^2 / 2 - 1) * ...
           signal / sqrt(sum(signal.^2) / sampleFrequency);

  % combine signal and noise
  data = signal + noise;

  % transform to one-sided frequency domain signal
  data = fft(data);
  data = data(1 : length(data) / 2 + 1);

  % apply q transform
  calibration = [];
  outlierFactor = 2.0;
  transform = qtransform(data, tiling, calibration, outlierFactor);

  % determine measured snr
  [ignore, plane] = min(abs(log(tiling.qs / q)));
  [ignore, row] = min(abs(log(tiling.planes{plane}.frequencies / ...
                              centerFrequency)));
  [ignore, tile] = min(abs(tiling.planes{plane}.rows{row}.times - centerTime));
  measuredSNR = ...
      sqrt(2 * transform{1}.planes{plane}.rows{row}.normalizedEnergies(tile));

  % report results
  fprintf(logFileID, '%#09.6f %9.3e %9.3e %9.3e %9.3e\n', ...
          centerTime, centerFrequency, q, snr, measuredSNR);

  % append result to results vectors
  centerTimes = [centerTimes centerTime];
  centerFrequencies = [centerFrequencies centerFrequency];
  qs = [qs q];
  injectedSNRs = [injectedSNRs snr];
  measuredSNRs = [measuredSNRs measuredSNR];

  % display scatter plot of measured vs. injected snr
  if watchFlag ~= 0,
    clf;
    set(gca, 'FontSize', 18);
    handle = loglog(injectedSNRs, max(measuredSNRs, 1), 'b.');
    axis([1 100 1 100]);
    grid on;
    set(handle, 'MarkerSize', 9);
    xlabel('Injected signal to noise ratio');
    ylabel('Measured signal to noise ratio');
    title('Measured vs. injected signal to noise ratio');
    drawnow;
  end

% end loop over trials
end

% save data
save testrecovery/testrecovery.mat ...
     tiling numberOfTrials centerTimes centerFrequencies qs ...
     injectedSNRs measuredSNRs

% display scatter plot of measured vs. injected snr
clf;
set(gca, 'FontSize', 18);
handle = loglog(injectedSNRs, max(measuredSNRs, 1), 'b.');
axis([1 100 1 100]);
grid on;
set(handle, 'MarkerSize', 9);
xlabel('Injected signal to noise ratio');
ylabel('Measured signal to noise ratio');
title('Measured vs. injected signal to noise ratio');
drawnow;
print -dpng testrecovery/testrecovery_1.png
unix(['convert -resize 500x500 testrecovery/testrecovery_1.png ' ...
      'testrecovery/testrecovery_1_thumbnail.png']);

% display histogram of measured vs. injected snr
clf;
set(gca, 'FontSize', 18);
hist2d(log10(injectedSNRs), log10(measuredSNRs), 'image', ...
       'catX', linspace(0, 2, 100), 'catY', linspace(0, 2, 100));
set(gca, 'XTick', log10([1 2 5 10 20 50 100]));
set(gca, 'YTick', log10([1 2 5 10 20 50 100]));
set(gca, 'XTickLabel', [1 2 5 10 20 50 100]);
set(gca, 'YTickLabel', [1 2 5 10 20 50 100]);
set(gca, 'TickDir', 'out');
set(gca, 'TickLength', [0.02 0.00]);
colormap jet;
handle = colorbar;
set(handle, 'FontSize', 18);
set(get(handle, 'YLabel'), 'String', 'Number of injections');
set(get(handle, 'YLabel'), 'FontSize', 18);
xlabel('Injected signal to noise ratio');
ylabel('Measured signal to noise ratio');
title('Measured vs. injected signal to noise ratio');
drawnow;
print -dpng testrecovery/testrecovery_2.png
unix(['convert -resize 500x500 testrecovery/testrecovery_2.png ' ...
      'testrecovery/testrecovery_2_thumbnail.png']);