function [transforms, channelNames] = qcollocate(transforms, tiling, ...
                                                channelNames)
% QCOLLOCATE Coherent Q transform sum and nulls for collocated detectors
%
% QCOLLOCATE computes the coherent Q transform sum and nulls of two or
% more collocated detectors.  The coherent sum is computed to maximize
% detectable signal to noise ratio, while the coherent nulls are
% computed to cancel out any signal.
%
% usage:
%
%  [transforms, channelNames] = qcollocate(transforms, tiling, channelNames);
%
%   transforms           cell array of detector Q transform structures
%   tiling               discrete Q tranform tiling structure from QTILE
%   channelNames         cell array of single detector channel names
%
%   transforms           cell array of coherent Q transform structures
%   channelNames         cell array of coherent transform channel names
%
% The results are returned as a cell array with the first cell containing
% the coherent sum transform and the remaining cells containing coherent
% null transforms.  The number of null transforms is one less than the
% number of detectors.
%
% The resulting 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
%   incoherentEnergies    vector of tile incoherent energies
%
% QCOLLOCATE also returns a modified list of coherent transform channel
% names based on original single detector channel names.
%
% See also QTRANSFORM, QTILE, QCONDITION, QTHRESHOLD, QSELECT, and QPIPELINE.

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

% Leo C. Stein
% lstein@ligo.mit.edu

% $Id: qcollocate.m,v 1.4 2007/01/19 23:56:49 shourov Exp $

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

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

% apply default arguments
if (nargin < 3) || isempty(channelNames),
  channelNames = [];
end

% force cell arrays
if ~iscell(transforms),
  transforms = mat2cell(transforms, size(transforms, 1), ...
                        size(transforms, 2));
end
if ~isempty(channelNames) & ~iscell(channelNames),
  channelNames = mat2cell(channelNames, size(channelNames, 1), ...
                          size(channelNames, 2));
end

% force one dimensional cell arrays
transforms = transforms(:);
channelNames = channelNames(:);

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

% provide default channel names
if isempty(channelNames),
  channelNames = cell(numberOfChannels, 1);
  for channelNumber = 1 : numberOfChannels,
    channelNames{channelNumber} = ['X' int2str(channelNumber) ':STRAIN'];
  end
end

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

% determine number of null transforms
numberOfNulls = numberOfChannels - 1;

% validate number of channels
if numberOfNulls < 1,
  error('at least two collocated detectors are required');
end

% 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 transform structures
for channelNumber = 1 : numberOfChannels,
  if ~strcmp(transforms{channelNumber}.id, ...
             'Discrete Q-transform transform structure'),
    error('input argument is not a discrete Q transform structure');
  end
end

% validate channel names
if ~isempty(channelNames) & (length(channelNames) ~= numberOfChannels),
  error(['channel names is inconsistent with number of transform channels']);
end

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

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

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

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

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                  determine coherent coefficients                         %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % collect vector of mean q transform energies for each channel
    meanEnergies = zeros(1, numberOfChannels);
    for channelNumber = 1 : numberOfChannels,
      meanEnergies(channelNumber) = ...
        transforms{channelNumber}.planes{plane}.rows{row}.meanEnergy;
    end

    % singular value decomposition of whitened response matrix
    [U, S, V] = svd(1 ./ sqrt(meanEnergies));

    % coherent sum and null stream weighting factors
    coherentWeightingFactors = V.' ./ ...
        (ones(numberOfChannels, 1) * sqrt(meanEnergies));
    incoherentWeightingFactors = (V.').^2;

    % initialize matrices of q transform coefficients and normalized energies
    transformCoefficients = ...
        zeros(numberOfChannels, tiling.planes{plane}.rows{row}.numberOfTiles);
    transformNormalizedEnergies = ...
        zeros(numberOfChannels, tiling.planes{plane}.rows{row}.numberOfTiles);

    % collect matrices of q transform coefficients and normalized energies
    for channelNumber = 1 : numberOfChannels,
      transformCoefficients(channelNumber, :) = ...
        transforms{channelNumber}.planes{plane}.rows{row}.coefficients;
      transformNormalizedEnergies(channelNumber, :) = ...
        transforms{channelNumber}.planes{plane}.rows{row}.normalizedEnergies;
    end

    % provide coherent sum in units of calibrated strain
    normalizationFactor = ones(numberOfChannels, 1);
    normalizationFactor(1) = 1 / S(1);

    % matrix of coherent transform coefficients
    coherentCoefficients = ...
        coherentWeightingFactors * transformCoefficients;

    % matrix of incoherent transform energies
    incoherentCoefficients = ...
        incoherentWeightingFactors * transformNormalizedEnergies;

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

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

      % insert coherent transform coefficients
      transforms{channelNumber}.planes{plane}.rows{row}.coefficients = ...
          coherentCoefficients(channelNumber, :) * ...
          normalizationFactor(channelNumber);

      % energy of coherent transform coefficients
      transforms{channelNumber}.planes{plane}.rows{row}.energies = ...
          (real(coherentCoefficients(channelNumber, :)).^2 + ...
           imag(coherentCoefficients(channelNumber, :)).^2) * ...
          normalizationFactor(channelNumber).^2;

      % mean energy of coherent transform coefficients
      transforms{channelNumber}.planes{plane}.rows{row}.meanEnergy = ...
          normalizationFactor(channelNumber).^2;

      % sigma energy of coherent transform coefficients
      transforms{channelNumber}.planes{plane}.rows{row}.sigmaEnergy = ...
          normalizationFactor(channelNumber).^2;

      % rayleigh statistic of coherent q transform
      transforms{channelNumber}.planes{plane}.rows{row}.rayleighStatistic = ...
          transforms{channelNumber}.planes{plane}.rows{row}.sigmaEnergy / ...
          transforms{channelNumber}.planes{plane}.rows{row}.meanEnergy;

      % normalized energies of coherent q transform
      transforms{channelNumber}.planes{plane}.rows{row}.normalizedEnergies = ...
          transforms{channelNumber}.planes{plane}.rows{row}.energies ./ ...
          normalizationFactor(channelNumber).^2;

      % incoherent energies of coherent q transform
      transforms{channelNumber}.planes{plane}.rows{row}.incoherentEnergies = ...
          incoherentCoefficients(channelNumber, :) * ...
          normalizationFactor(channelNumber).^2;

    % end loop over coherent channels
    end

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

  % end loop over frequency rows
  end

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

% end loop over Q planes
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             modify channel names                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% basename of colloated channels
channelBaseName = [channelNames{1}(1) ':COHERENT'];

% name of coherent sum channel
channelNames{1} = [channelBaseName '-SUM'];

% name of null stream channels
for channelNumber = 2 : numberOfChannels,
  channelNames{channelNumber} = [channelBaseName '-NULL'];
  if numberOfChannels > 2,
    channelNames{channelNumber} = [channelNames{channelNumber} '_' ...
                                   int2str(channelNumber - 1)];
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   return collocated Q transform structures                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% return to calling function
return