function [rawData, highPassedData, whitenedData, calibration] = ...
    qscancondition(rawData, tiling, responseFunctions);
% QSCANCONDITION High pass filter and whiten time series data for Q scans
%
% QSCANCONDITION is a modified version of QCONDITION that also returned
% intermediate results required by QSCAN.
%
% QSCANCONDITION high pass filters and whitens time series data prior to
% analysis by the Q transform.  The data is first zero-phase high pass filtered
% at the minimum analysis frequency specified in the tiling structure.  The
% resulting high pass filtered data is then whitened by zero-phase linear
% prediction at a frequency resolution equal to the minimum analysis bandwidth
% requested in the tiling structure.  Note that the resulting whitened data is
% returned as a frequency series, not a time series.  In addition, the returned
% frequency series extend from zero frequency to the Nyquist frequency.  As a
% result, they are of length N / 2 + 1, where N is the length of the individual
% input time series.
%
% QSCANCONDITION also returns information on the estimated length of filter
% transients.  Note that, due to use of zero-phase filtering, these transients
% are present at both the beginning and end of the data stream.
%
% usage:
%
% [rawData, highPassedData, whitenedData, calibration] = ...
%     qscancondition(rawData, tiling, responseFunctions);
%
%   rawData               cell array of input time series
%   tiling                discrete Q transform tiling structure from QTILE
%   responseFunctions     cell array of calibration response functions
%
%   rawData               cell array of unconditioned frequency series
%   highPassedData        cell array of high passed filtered frequency series
%   whitenedData          cell array of whitened frequency series data
%   calibration           cell array of Q transform calibration structures
%
% The resulting calibration structures are parallel to the structure returned
% by QTILE and contain the following supplemental fields for each frequency
% row.
%
%   calibrationFactor     row specific amplitude correction factors
%
% If provided, responseFunctions should consist of a cell array of complex
% valued frequency response vectors, with one cell per channel.  These vectors
% should extend from zero frequency to the Nyquist frequency at a frequency
% resolution equal to the reciprocal of the data duration.  As a result they
% should be of length N / 2 + 1, where N is the length of the individual input
% time series.  By default, QSCANCONDITION assumes that the raw data is already
% calibrated and that responseFunctions is unity.
%
% See also QTILE, QCONDITION, QTRANSFORM, QTHRESHOLD, QSELECT, and QSCAN.

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

% $Id: qscancondition.m,v 1.2 2007/08/10 19:02:30 shourov Exp $

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

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

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

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

% provide default response functions
if isempty(responseFunctions),
  responseFunctions = cell(size(rawData));
end

% force cell array of response functions
if ~iscell(responseFunctions),
  responseFunctions = mat2cell(responseFunctions, ...
                               size(responseFunctions, 1), ...
                               size(responseFunctions, 2));
end

% force one dimensional cell array
rawData = rawData(:);
responseFunctions = responseFunctions(:);

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

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

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

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

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

% validate calibration response channels
if length(responseFunctions) ~= numberOfChannels,
  error('number of channels not consistent with calibration response');
end

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

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

% create empty cell array of Q transform calibration structures
calibration = cell(size(rawData));

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

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

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

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

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

  % end loop over Q planes
  end

% end loop over channels
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             hardcoded parameters                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% high pass filter order
hpfOrder = 12;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              derived parameters                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% nyquist frequency
nyquistFrequency = tiling.sampleFrequency / 2;

% linear predictor error filter order
lpefOrder = ceil(tiling.sampleFrequency * tiling.whiteningDuration);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                 window data                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% design window
tukeyWindow = ...
    tukeywin(dataLength, 2 * ...
             max(tiling.transientDuration - tiling.whiteningDuration, 0) / ...
             tiling.duration).';

% apply window
for channelNumber = 1 : numberOfChannels,
  rawData{channelNumber} = rawData{channelNumber} .* tukeyWindow;
end

% release window memory
clear tukeyWindow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               high pass filter                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% design high pass filter
[hpfZeros, hpfPoles, hpfGain] = butter(hpfOrder, tiling.highPassCutoff / ...
                                       nyquistFrequency, 'high');
hpfSOS = zp2sos(hpfZeros, hpfPoles, hpfGain);

% apply high pass filter
for channelNumber = 1 : numberOfChannels,
  highPassedData{channelNumber} = sosfiltfilt(hpfSOS, rawData{channelNumber});
end

% supress high pass filter transients
for channelNumber = 1 : numberOfChannels,
  rawData{channelNumber}(1 : lpefOrder) = ...
      zeros(1, lpefOrder);
  rawData{channelNumber}(dataLength - lpefOrder + 1 : dataLength) = ...
      zeros(1, lpefOrder);
  highPassedData{channelNumber}(1 : lpefOrder) = ...
      zeros(1, lpefOrder);
  highPassedData{channelNumber}(dataLength - lpefOrder + 1 : dataLength) = ...
      zeros(1, lpefOrder);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            fast fourier transform                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% fourier transform
for channelNumber = 1 : numberOfChannels,
  rawData{channelNumber} = fft(rawData{channelNumber});
  highPassedData{channelNumber} = fft(highPassedData{channelNumber});
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        linear predictor error filter                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% auto-correlation estimates
for channelNumber = 1 : numberOfChannels,
  autoCorrelation{channelNumber} = real(highPassedData{channelNumber}).^2 + ...
                                   imag(highPassedData{channelNumber}).^2;
end
for channelNumber = 1 : numberOfChannels,
  autoCorrelation{channelNumber} = real(ifft(autoCorrelation{channelNumber}));
end

% solve for lpef coefficients
filterLength = lpefOrder + 1;
for channelNumber = 1 : numberOfChannels,
  lpefCoefficients{channelNumber} = ...
      levinson(autoCorrelation{channelNumber}(1 : filterLength), lpefOrder);
end

% release auto-correlation vector memory
clear autoCorrelation;

% create zero-phase lpef
for channelNumber = 1 : numberOfChannels,
  lpefCoefficients{channelNumber} = [lpefCoefficients{channelNumber} ...
                                     zeros(1, dataLength - filterLength)];
end
for channelNumber = 1 : numberOfChannels,
  lpefCoefficients{channelNumber} = fft(lpefCoefficients{channelNumber});
end
for channelNumber = 1 : numberOfChannels,
  lpefCoefficients{channelNumber} = abs(lpefCoefficients{channelNumber});
end

% extract zero to nyquist frequency data and lpef coefficients
for channelNumber = 1 : numberOfChannels,
  rawData{channelNumber} = ...
      rawData{channelNumber}(1 : halfDataLength);
  highPassedData{channelNumber} = ...
      highPassedData{channelNumber}(1 : halfDataLength);
  lpefCoefficients{channelNumber} = ...
      lpefCoefficients{channelNumber}(1 : halfDataLength);
end

% apply zero-phase lpef
for channelNumber = 1 : numberOfChannels,
  whitenedData{channelNumber} = lpefCoefficients{channelNumber} .* ...
      highPassedData{channelNumber};
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                             apply calibration                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% magnitude response of high pass filter
minimumFrequencyStep = 1 / tiling.duration;
frequencies = 0 : minimumFrequencyStep : nyquistFrequency;
hpfArgument = (frequencies / tiling.highPassCutoff).^(2 * hpfOrder);
hpfResponse = hpfArgument ./ (1 + hpfArgument);

% compute calibration correction factors
for channelNumber = 1 : numberOfChannels,
  combinedResponse = hpfResponse .* lpefCoefficients{channelNumber};
  if ~isempty(responseFunctions{channelNumber}),
    combinedResponse = combinedResponse .* ...
        abs(responseFunctions{channelNumber});
  end
  % smoothingBandwidth = 100 / tiling.duration;
  % smoothingLength = ceil(smoothingBandwidth / minimumFrequencyStep);
  % combinedResponse = smooth(combinedResponse, smoothingLength, 'moving').';
  for plane = 1 : tiling.numberOfPlanes,
    for row = 1 : tiling.planes{plane}.numberOfRows,
      calibrationIntegrand = (tiling.planes{plane}.rows{row}.window ./ ...
          combinedResponse(tiling.planes{plane}.rows{row}.dataIndices)).^2;
      calibration{channelNumber}.planes{plane}.rows{row}.calibrationFactor = ...
          sqrt(sum(calibrationIntegrand * minimumFrequencyStep) / 2) * ...
          dataLength / tiling.planes{plane}.rows{row}.numberOfTiles;
    end
  end
end

% apply calibration phase
for channelNumber = 1 : numberOfChannels,
  if ~isempty(responseFunctions{channelNumber}),
    rawData{channelNumber} = rawData{channelNumber} .* ...
        exp(-sqrt(-1) * angle(responseFunctions{channelNumber}));
    highPassedData{channelNumber} = highPassedData{channelNumber} .* ...
        exp(-sqrt(-1) * angle(responseFunctions{channelNumber}));
    whitenedData{channelNumber} = whitenedData{channelNumber} .* ...
        exp(-sqrt(-1) * angle(responseFunctions{channelNumber}));
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       release lpef coefficient memory                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% release lpef coefficient memory
clear lpefCoefficients;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                    return                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% return to calling function
return;