h2asq_fd

 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 counts.

   output = h2asq_fd(fname_olg,fname_sf,alpha,beta,input,fs)

 fname_olg  
         File containing frequency (Hz), modulus (dimensionless), and 
         phase (radians) of the open loop gain.
 fname_sf   
         File containing frequency (Hz), modulus (AS_Q counts/strain),  
         and phase (radians) of the sensing function.
 alpha   Calibration parameter alpha.
 beta    Calibration parameter beta.
 input   Input timeseries data in strain.
 fs      sampling rate of input timeseries.
 
 output  Output timeseries data in AS_Q 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.

0001 function [output] = h2asq_fd(fname_olg,fname_sf,alpha,beta,input,fs)
0002 % H2ASQ_FD - Apply the inverse calibration to the specified timeseries data.
0003 % This function reads an input timeseries in strain and using the specified
0004 % calibration converts it to AS_Q counts.
0005 %
0006 %   output = h2asq_fd(fname_olg,fname_sf,alpha,beta,input,fs)
0007 %
0008 % fname_olg
0009 %         File containing frequency (Hz), modulus (dimensionless), and
0010 %         phase (radians) of the open loop gain.
0011 % fname_sf
0012 %         File containing frequency (Hz), modulus (AS_Q counts/strain),
0013 %         and phase (radians) of the sensing function.
0014 % alpha   Calibration parameter alpha.
0015 % beta    Calibration parameter beta.
0016 % input   Input timeseries data in strain.
0017 % fs      sampling rate of input timeseries.
0018 %
0019 % output  Output timeseries data in AS_Q counts.
0020 %
0021 % Notes:
0022 % 1. Linear interpolation is used to find the response function at
0023 %    frequencies desired.  At some frequencies (particularly low)
0024 %    the magnitude of the interpolated reponse is too low.  This could
0025 %    be fixed using a higher-order interpolation (spline).  However,
0026 %    since the reference calibrations are normally supplied at sub-1Hz
0027 %    resolutions, I don't see the need to do anything here.
0028 % 2. Any DC data in the reference calibrations is ignored (contains
0029 %    NaNs).
0030 % 3. Acausal filtering.  Get ringing before start time of an impulsive
0031 %    input.  The peak time is after the impulse time, though.
0032 
0033 % $Id: h2asq_fd.m,v 1.2 2004/01/15 19:31:20 thorne Exp $
0034 % 2003/11/07
0035 %----- Derived quantities.
0036 %---------- Number of data points
0037 N = length(input);
0038 %---------- Frequency resolution
0039 dF = fs/N;  
0040 %---------- Vector holding frequencies in usual screwy FFT order:
0041 %           vector element:   [ 1  2  ...  N/2-1  N/2    N/2+1            N/2+2   ... N-1  N ]
0042 %           frequency (dF):   [ 0  1  ...  N/2-2  N/2-1  (N/2 or -N/2)   -N/2+1  ... -2   -1 ]
0043 F = [ 0:N/2 , -N/2+1:-1 ]'*dF;
0044 %---------- Fourier transform
0045 input_f = fft(input);
0046 
0047 %----- Load in reference calibrations.
0048 if(exist(fname_olg)>0)
0049     [fH, mH, pH]=ctextread(fname_olg,'%f %f %f');
0050 else
0051     fprintf(' Oops - OpenLoopGain file %s not found\n',fname_olg);
0052     exit
0053 end    
0054 if(exist(fname_sf)>0)
0055     [fC, mC, pC]=ctextread(fname_sf,'%f %f %f');
0056 else
0057     fprintf(' Oops - Sensing Function file %s not found\n',fname_sf);
0058     exit
0059 end        
0060 if (~isequal(fH,fC)) 
0061     error('The open-loop gain and sensing function must be sampled at the same frequencies') 
0062 end
0063 %----- Often reference cal files have an f=0 line, which will contain NaNs.
0064 if (fH(1)==0)
0065     fH(1) = [];
0066     mH(1) = [];
0067     pH(1) = [];
0068     fC(1) = [];
0069     mC(1) = [];
0070     pC(1) = [];
0071 end
0072 H = mH.*exp(i*pH);
0073 C = mC.*exp(i*pC);
0074 %----- Response function (AS_Q/strain)
0075 Rtemp = alpha*C./(1+alpha*beta*H);
0076 
0077 %----- fH typically doesn't cover full range of frequencies.  Select subset
0078 %      of frequencies covered by the reference calibrations and interpolate
0079 %      to them.
0080 k = find(F>=fH(1) & F<=fH(end));
0081 Fcal = F(k);
0082 Rcal = interp1(fH,Rtemp,Fcal);
0083 %Rcal = interp1(fH,Rtemp,Fcal,'spline');
0084 %figure; loglog(fH,abs(Rtemp),Fcal,abs(Rcal),'+');
0085 %----- Copy interpolated response function into new holder of same size as
0086 %      data.  Leave response function as zero at frequencies for
0087 %      which we have no calibration data.
0088 R = zeros(size(F));
0089 R(k) = Rcal;
0090 %----- Sharp drop in response function produces ringing in
0091 %      time domain at f ~ f_drop (around 4kHz).  Found can eliminate by
0092 %      making response function go smoothly to zero at DC and Nyquist
0093 %      frequencies.  (Smoothing at DC seems to have no effect, but do
0094 %      anyway.)
0095 %---------- DC side.
0096 slope = R(k(1))/(k(1)-1);
0097 for j=2:k(1)-1
0098     R(j) = slope*(j-1);
0099 end
0100 %---------- Nyquist side.
0101 slope = -R(k(end))/(N/2+1-k(end));
0102 for j=k(end)+1:N/2+1
0103     R(j) = R(k(end))+slope*(j-k(end));
0104 end
0105 %----- Set R(-f) = R*(f) (real filter).
0106 R(N/2+2:N) = conj(R(N/2:-1:2));
0107 %----- Plot original and inrepolated data.
0108 %figure;
0109 %subplot(2,1,1)
0110 %plot(fH,abs(Rtemp),'+',F,abs(R),'markersize',4);
0111 %grid
0112 %subplot(2,1,2)
0113 %plot(fH,angle(Rtemp),'+',F,angle(R),'markersize',4);
0114 %plot(fH,unwrap(angle(Rtemp)),'+',F,unwrap(angle(R)),'markersize',4);
0115 %grid
0116 
0117 %----- Calibrate input data in frequency domain.
0118 output_f = input_f.*R;
0119 output = ifft(output_f);
0120 %---------- Cast off any imaginary part (introduced through roundoff errors).
0121 output = real(output);
0122 
0123 %----- Plot powers.
0124 %[Pinput,w] = pwelch(input);
0125 %[Poutput,w] = pwelch(output);
0126 %R_chk = interp1(fH,Rtemp,w*fs/2/pi);
0127 %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));
0128 %grid
0129 
0130 return