function transforms = qtransform(data, tiling, calibration, outlierFactor);
% QTRANSFORM Discrete Q transform
%
% QTRANSFORM applies the discrete Q transform described by the predetermined
% tiling structure to frequency domain data.  The tiling structure contains
% the transform parameters and is generated by the QTILE function.  The input
% data should be the fourier transform of the time series data to be analyzed.
% However, the input frequency series should only extend from zero frequency
% to the Nyquist frequency.  As a result, the input frequency series should be
% of length N / 2 + 1, where N is the length of the original input time series.
%
% usage: transforms = qtransform(data, tiling, calibration, outlierFactor);
%
%   data                 cell array of input frequency series data
%   tiling               discrete Q tranform tiling structure from QTILE
%   calibration          cell array of calibration structures from QCONDITION
%   outlierFactor        Tukey whisker multiplier for outlier rejection
%
%   transforms           cell array of discrete Q transform structures
%
% The resulting discrete Q transform structures are parallel to the structure
% returned by QTILE and contain the following supplemental fields for each
% frequency row.
%
%   coefficients          vector of complex valued tile coefficients
%   energies              vector of tile energies
%   meanEnergy            mean of tile energies
%   sigmaEnergy           standard deviation of tile energies
%   normalizedEnergies    vector of normalized tile energies
%   rayleighStatistic     ratio of standard deviation to mean tile energy
%
% In addition, the mean energy and Rayleigh statistic from each frequency row
% are collected for convenience and returned as the following two supplemental
% fields in each plane.
%
%   powerSpectrum         vector of mean energies from each frequency row
%   rayleighStatistics    ratio of standard deviation to mean tile energies
%
% No calibration correction factors are applied if the calibration argument
% is the empty matrix [].
%
% See also QTILE, QCONDITION, QTHRESHOLD, QSELECT, and QPIPELINE.

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

% $Id: qtransform.m,v 1.2 2006/10/16 03:27:20 shourov Exp $

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        process command line arguments                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% verify correct number of input arguments
error(nargchk(4, 4, nargin));

% force cell array
if ~iscell(data),
  data = mat2cell(data, size(data, 1), size(data, 2));
end

% force one dimensional cell array
data = data(:);

% determine number of channels
numberOfChannels = length(data);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       validate command line arguments                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% validate tiling structure
if ~strcmp(tiling.id, 'Discrete Q-transform tile structure'),
  error('input argument is not a discrete Q transform tiling structure');
end

% validate calibration structures
if ~isempty(calibration),
  if length(calibration) ~= numberOfChannels,
    error('number of channels not consistent with calibration structure');
  end
  for channelNumber = 1 : numberOfChannels,
    if ~strcmp(calibration{channelNumber}.id, ...
               'Discrete Q-transform calibration structure'),
      error('input argument is not a discrete Q transform calibration structure');
    end
  end
end

% determine required data lengths
dataLength = tiling.sampleFrequency * tiling.duration;
halfDataLength = dataLength / 2 + 1;

% validate data lengthand force row vectors
for channelNumber = 1 : numberOfChannels,
  data{channelNumber} = data{channelNumber}(:).';
  if length(data{channelNumber}) ~= halfDataLength,
    error('data length not consistent with tiling');
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                      initialize Q transform structures                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% create empty cell array of Q transform structures
transforms = cell(size(data));

% begin loop over channels
for channelNumber = 1 : numberOfChannels,

  % insert structure identification string
  transforms{channelNumber}.id = 'Discrete Q-transform transform structure';

  % create empty cell array of Q plane structures
  transforms{channelNumber}.planes = cell(size(tiling.planes));

  % begin loop over Q planes
  for plane = 1 : tiling.numberOfPlanes,

    % create empty cell array of frequency row structures
    transforms{channelNumber}.planes{plane}.rows = ...
        cell(size(tiling.planes{plane}.numberOfRows));

  % end loop over Q planes
  end

% end loop over channels
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           begin loop over Q planes                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% begin loop over Q planes
for plane = 1 : tiling.numberOfPlanes,

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                    initialize Q plane result vectors                       %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % begin loop over channels
  for channelNumber = 1 : numberOfChannels,

    % initialize power spectrum field in Q plane structure
    transforms{channelNumber}.planes{plane}.powerSpectrum = [];

    % initialize rayleigh statistics field in Q plane structure
    transforms{channelNumber}.planes{plane}.rayleighStatistics = [];

  % end loop over channels
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                      begin loop over frequency rows                        %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % begin loop over frequency rows
  for row = 1 : tiling.planes{plane}.numberOfRows,

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %               extract and window frequency domain data                   %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % number of zeros to pad at negative frequencies
    leftZeroPadLength = (tiling.planes{plane}.rows{row}.zeroPadLength - 1) / 2;

    % number of zeros to pad at positive frequencies
    rightZeroPadLength = (tiling.planes{plane}.rows{row}.zeroPadLength + 1) / 2;

    % begin loop over channels
    for channelNumber = 1 : numberOfChannels,

      % extract and window in-band data
      windowedData{channelNumber} = tiling.planes{plane}.rows{row}.window .* ...
          data{channelNumber}(tiling.planes{plane}.rows{row}.dataIndices);

      % zero pad windowed data
      windowedData{channelNumber} = [zeros(1, leftZeroPadLength) ...
                               windowedData{channelNumber} ...
                               zeros(1, rightZeroPadLength)];

      % reorder indices for fast fourier transform
      windowedData{channelNumber} = ...
          windowedData{channelNumber}([(end / 2 : end) (1 : end / 2 - 1)]);

    % end loop over channels
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %               inverse fourier transform windowed data                    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % begin loop over channels
    for channelNumber = 1 : numberOfChannels,

    % complex valued tile coefficients
    coefficients{channelNumber} = ifft(windowedData{channelNumber});

    % apply calibration correction factor to coefficients
    if ~isempty(calibration),
      coefficients{channelNumber} = coefficients{channelNumber} * ...
          calibration{channelNumber}.planes{plane}.rows{row}.calibrationFactor;
    end

    % tile energies
    energies{channelNumber} = real(coefficients{channelNumber}).^2 + ...
                              imag(coefficients{channelNumber}).^2;

    % end loop over channels
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %        exclude outliers and filter transients from statistics            %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % begin loop over channels
    for channelNumber = 1 : numberOfChannels,

      % indices of non-transient tiles
      validIndices{channelNumber} = ...
          find((tiling.planes{plane}.rows{row}.times > ...
                tiling.transientDuration) & ...
               (tiling.planes{plane}.rows{row}.times < ...
                tiling.duration - tiling.transientDuration));

      % identify lower and upper quartile energies
      sortedEnergies = ...
          sort(energies{channelNumber}(validIndices{channelNumber}));
      lowerQuartile{channelNumber} = ...
          sortedEnergies(round(0.25 * length(validIndices{channelNumber})));
      upperQuartile{channelNumber} = ...
          sortedEnergies(round(0.75 * length(validIndices{channelNumber})));

      % determine inter quartile range
      interQuartileRange{channelNumber} = upperQuartile{channelNumber} - ...
                                          lowerQuartile{channelNumber};

      % energy threshold of outliers
      outlierThreshold{channelNumber} = upperQuartile{channelNumber} + ...
          outlierFactor * interQuartileRange{channelNumber};

      % indices of non-outlier and non-transient tiles
      validIndices{channelNumber} = ...
          find((energies{channelNumber} < ...
                outlierThreshold{channelNumber}) & ...
               (tiling.planes{plane}.rows{row}.times > ...
                tiling.transientDuration) & ...
               (tiling.planes{plane}.rows{row}.times < ...
                tiling.duration - tiling.transientDuration));

    % end loop over channels
    end

    % for reasonable outlier factors,
    if outlierFactor < 100,

      % mean energy correction factor for outlier rejection bias
      meanCorrectionFactor = (4 * 3^outlierFactor - 1) / ...
                             ((4 * 3^outlierFactor - 1) - ...
                             (outlierFactor * log(3) + log(4)));

      % sigma energy correction factor for outlier rejection bias
      sigmaCorrectionFactor = (4 * 3^outlierFactor - 1)^2 / ...
                              ((4 * 3^outlierFactor - 1)^2 - ...
                              (outlierFactor * log(3) + log(4))^2 * ...
                              4 * 3^outlierFactor);

    % otherwise, for large outlier factors
    else

      % mean energy correction factor for outlier rejection bias
      meanCorrectionFactor = 1;


      % sigma energy correction factor for outlier rejection bias
      sigmaCorrectionFactor = 1;

    % continue
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %          determine tile statistics and normalized energies               %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % begin loop over channels
    for channelNumber = 1 : numberOfChannels,

      % mean of valid tile energies
      meanEnergy{channelNumber} = ...
          mean(energies{channelNumber}(validIndices{channelNumber}));

      % correct for bias due to outlier rejection
      meanEnergy{channelNumber} = meanEnergy{channelNumber} * ...
          meanCorrectionFactor;

      % standard deviation of valid tile energies
      sigmaEnergy{channelNumber} = ...
          std(energies{channelNumber}(validIndices{channelNumber}));

      % correct for bias due to outlier rejection
      sigmaEnergy{channelNumber} = sigmaEnergy{channelNumber} * ...
          sigmaCorrectionFactor;

      % rayleigh statistic of frequency bin
      rayleighStatistic{channelNumber} = sigmaEnergy{channelNumber} / ...
          meanEnergy{channelNumber};

      % normalized tile energies
      normalizedEnergies{channelNumber} = energies{channelNumber} / ...
          meanEnergy{channelNumber};

    % end loop over channels
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %               insert results into transform structure                    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % begin loop over channels
    for channelNumber = 1 : numberOfChannels,

    % insert tile coefficients into frequency row structure
    transforms{channelNumber}.planes{plane}.rows{row}.coefficients = ...
        coefficients{channelNumber};

    % insert tile energies into frequency row structure
    transforms{channelNumber}.planes{plane}.rows{row}.energies = ...
        energies{channelNumber};

    % insert mean tile energy into frequency row structure
    transforms{channelNumber}.planes{plane}.rows{row}.meanEnergy = ...
        meanEnergy{channelNumber};

    % insert tile energy standard deviation into frequency row structure
    transforms{channelNumber}.planes{plane}.rows{row}.sigmaEnergy = ...
        sigmaEnergy{channelNumber};

    % insert rayleigh statistic of row into frequency row structure
    transforms{channelNumber}.planes{plane}.rows{row}.rayleighStatistic = ...
        rayleighStatistic{channelNumber};

    % insert normalized tile energies into frequency row structure
    transforms{channelNumber}.planes{plane}.rows{row}.normalizedEnergies = ...
        normalizedEnergies{channelNumber};

    % append mean tile energy to power spectrum of Q plane
    transforms{channelNumber}.planes{plane}.powerSpectrum = ...
       [transforms{channelNumber}.planes{plane}.powerSpectrum; ...
        meanEnergy{channelNumber};];

    % append rayleigh statistic to rayleigh statistics field of Q plane
    transforms{channelNumber}.planes{plane}.rayleighStatistics = ...
        [transforms{channelNumber}.planes{plane}.rayleighStatistics; ...
         rayleighStatistic{channelNumber};];

    % end loop over channels
    end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                       end loop over frequency rows                         %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % end loop over frequency rows
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            end loop over Q planes                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% end loop over Q planes
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    return discrete Q transform structure                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% return to calling function
return