function responseFunctions = qreadcalibration(calibrationFiles, tiling, ...
                                              debugLevel)
% QREADCALBIRATION Read calbiration response functions for multiple detectors
%
% QREADCALIBRATION reads the frequency domain calibration response functions
% from the specified calibration text files.  The calibraiton text files are
% assumed to be in the form provided by the calibration team.  They consist
% of three columns of ASCII text specifying the frequency, magnitude, and
% phase of the response function.
%
% usage:
%
%   responseFunctions = qreadcalibration(calibrationFiles, tiling, debugLevel);
%
%     calibrationFiles        cell array of calibration file names
%     tiling                  q transform tiling structure
%     debugLevel              verboseness of debug output
%
%     responseFunctions       cell array of response functions
%
% QREADCALIBRATION constructs a one sided complex valued response function
% from the magnitude and phase data.  If necessary, this response function
% is extended from zero frequency to the Nyquist frequency of the analysis
% and interpolated to the frequency resolution of the analysis.  The
% resulting response functions are then returned as a cell array of row
% vectors in the same order as the cell array of calibration names.
%
% See also QTILE, QREADDATA, QCONDITION, QTRANSFORM, QTHRESHOLD, QSELECT,
% and QEXAMPLE.

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

% $Id: qreadcalibration.m,v 1.2 2007/04/02 05:03:42 shourov Exp $

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

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

% apply default arguments
if (nargin < 3) || isempty(debugLevel),
  debugLevel = 1;
end

% force cell array of calibration names
if ~isempty(calibrationFiles) & ~iscell(calibrationFiles),
  calibrationFiles = mat2cell(calibrationFiles, size(calibrationFiles, 1), ...
                              size(calibrationFiles, 2));
end

% force one dimensional cell arrays
calibrationFiles = calibrationFiles(:);

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

% check for valid discrete Q-transform tile structure
if ~strcmp(tiling.id, 'Discrete Q-transform tile structure'),
  error('The first argument is not a valid Q-transform tiling.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                      response function frequency vector                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% vector of frequencies required by tiling
responseFrequencies = 0 : 1 / tiling.duration : tiling.sampleFrequency / 2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         initialize result structures                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% initialize cell array of response functions
responseFunctions = cell(numberOfChannels, 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                  read data                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

  % if calibration requested
  if ~strcmp(upper(calibrationFiles{channelNumber}), 'NONE'),

    % report status
    if debugLevel >=2,
      fprintf(1, 'reading %s\n', calibrationFiles{channelNumber});
    end

    % read calibration data
    [frequency, magnitude, phase] = ...
        dataread('file', calibrationFiles{channelNumber}, ...
                 '%n%n%n', 'commentstyle', 'matlab');

    % test for error reading calibration data
    if length(frequency) < 2,
      error('error reading %s\n', calibrationFiles{channelNumber});
    end

    % discard NaNs and Infs
    keepIndices = find(isfinite(magnitude) & isfinite(phase));
    frequency = frequency(keepIndices);
    magnitude = magnitude(keepIndices);
    phase = phase(keepIndices);

    % interpolate response function to correspond to tiling
    magnitude = interp1(frequency, magnitude, responseFrequencies, ...
                        'nearest', 'extrap');
    phase = interp1(frequency, phase, responseFrequencies, ...
                    'nearest', 'extrap');

    % construct complex valued response function
    responseFunctions{channelNumber} = magnitude .* exp(sqrt(-1) * phase);

    % force row data
    responseFunctions{channelNumber} = responseFunctions{channelNumber}(:);

  % end test for requested calibration
  end

% end loop over channels
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          return to calling function                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% return to calling function
return;