function [output] = h2asq_fd(IFO,gpsSec,input,fs,sigType)
% H2ASQ_FD - Apply the inverse calibration to the specified timeseries data.
% This function reads an input timeseries in strain and using the specified 
% calibration converts it to AS_Q/DARM_ERR counts.
%
%   output = h2asq_fd(IFO,gpsSec,input,fs[,sigType])
%
%       IFO - interferometer
%       gpsSec - GPS time
%       input   Input timeseries data in strain.
%       fs      sampling rate of input timeseries.
%       sigType [OPTIONAL] signal type (AS_Q,DARM_ERR)
%               <default = AS_Q>
%
%   output  Output timeseries data in gw channel counts.
%
% Notes: 
% 1. Linear interpolation is used to find the response function at 
%    frequencies desired.  At some frequencies (particularly low) 
%    the magnitude of the interpolated reponse is too low.  This could
%    be fixed using a higher-order interpolation (spline).  However, 
%    since the reference calibrations are normally supplied at sub-1Hz 
%    resolutions, I don't see the need to do anything here.
% 2. Any DC data in the reference calibrations is ignored (contains 
%    NaNs).
% 3. Acausal filtering.  Get ringing before start time of an impulsive 
%    input.  The peak time is after the impulse time, though.
%
% $Id: h2asq_fd.m,v 1.8 2005/04/07 01:01:55 thorne Exp $
% 2003/11/07

%----- Derived quantities.
%---------- Number of data points
N = length(input);
%---------- Frequency resolution
dF = fs/N;  
%---------- Vector holding frequencies in usual screwy FFT order:
%           vector element:   [ 1  2  ...  N/2-1  N/2    N/2+1            N/2+2   ... N-1  N ]
%           frequency (dF):   [ 0  1  ...  N/2-2  N/2-1  (N/2 or -N/2)   -N/2+1  ... -2   -1 ]
F = [ 0:N/2 , -N/2+1:-1 ]'*dF;
%---------- Fourier transform
input_f = fft(input);

% GET alpha,beta values for this time
% AVERAGE these over the interval
%      - since alpha-beta values are updated every 60 sec, 
%       there is likely only one time used
% GET open loop gain as a function of frequency
% GET sensing function as a function of frequency
% IF frequency steps do not match
%   CREATE error message
% ENDIF
output = [];
DEF_SIG = 'AS_Q';
if(nargin > 4)
    if( (strcmp(sigType,'AS_Q')==false) && ...
        (strcmp(sigType,'DARM_ERR') == false))
        msgId = 'h2asq_fd:badSignal';
        warning(msgId,'%s: signal type %s is not valid - assume AS_Q',msgId,sigType);
        sigType = DEF_SIG;
    end
else
    sigType = DEF_SIG;
end
[gpsList,alpha_t,beta_t] = getalphabeta(IFO,gpsSec,1);
alpha = mean(alpha_t);
beta = mean(beta_t);
[fH,H] = getopenloopgain(IFO,gpsSec);
[fC,C] = getsensefunction(IFO,gpsSec,sigType);
if (~isequal(fH,fC))
    error('The open-loop gain and sensing function must be sampled at the same frequencies');
end

%----- Often reference cal files have an f=0 line, which will contain NaNs.
if (fH(1)==0)
    fH(1) = [];
    H(1) = [];
    fC(1) = [];
    C(1) = [];
end

% Choose response function from signal type
switch sigType
    case 'AS_Q'
%----- Response function (AS_Q/strain)
        Rtemp = alpha*C./(1+alpha*beta*H);
    case 'DARM_ERR'
%----- Response function (DARM_ERR/strain)
        Rtemp = alpha*beta*C./(1+alpha*beta*H);
    otherwise
%----- Response function (AS_Q/strain)
        Rtemp = alpha*C./(1+alpha*beta*H);
end

%----- fH typically doesn't cover full range of frequencies.  Select subset  
%      of frequencies covered by the reference calibrations and interpolate
%      to them.  
k = find(F>=fH(1) & F<=fH(end));
Fcal = F(k);
Rcal = interp1(fH,Rtemp,Fcal);
%Rcal = interp1(fH,Rtemp,Fcal,'spline');
%figure; loglog(fH,abs(Rtemp),Fcal,abs(Rcal),'+');
%----- Copy interpolated response function into new holder of same size as  
%      data.  Leave response function as zero at frequencies for 
%      which we have no calibration data.
R = zeros(size(F));
R(k) = Rcal;
%----- Sharp drop in response function produces ringing in 
%      time domain at f ~ f_drop (around 4kHz).  Found can eliminate by 
%      making response function go smoothly to zero at DC and Nyquist 
%      frequencies.  (Smoothing at DC seems to have no effect, but do 
%      anyway.)
%---------- DC side.
slope = R(k(1))/(k(1)-1);
for j=2:k(1)-1
    R(j) = slope*(j-1);
end
%---------- Nyquist side.
% ** NOTE: we only need to make this correction if our frequency range 
%   is beyond the Nyquist of the calibration data.
%  Otherwise, we get a nasty divide-by-zero
if( k(end) < ((N/2)+1))
    slope = -R(k(end))/(N/2+1-k(end));
    for j=k(end)+1:N/2+1
        R(j) = R(k(end))+slope*(j-k(end));
    end
end
%----- Set R(-f) = R*(f) (real filter).
R(N/2+2:N) = conj(R(N/2:-1:2));
%----- Plot original and inrepolated data.
%figure; 
%subplot(2,1,1)
%plot(fH,abs(Rtemp),'+',F,abs(R),'markersize',4);
%grid
%subplot(2,1,2)
%plot(fH,angle(Rtemp),'+',F,angle(R),'markersize',4);
%plot(fH,unwrap(angle(Rtemp)),'+',F,unwrap(angle(R)),'markersize',4);
%grid

%----- Calibrate input data in frequency domain.
output_f = input_f.*R;
output = ifft(output_f);
%---------- Cast off any imaginary part (introduced through roundoff errors).
output = real(output);

%----- Plot powers.
%[Pinput,w] = pwelch(input);
%[Poutput,w] = pwelch(output);
%R_chk = interp1(fH,Rtemp,w*fs/2/pi);
%figure; loglog(w*fs/2/pi,Pinput.^0.5,w*fs/2/pi,Poutput.^0.5,w*fs/2/pi,Pinput.^0.5.*abs(R_chk));
%grid

return
