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/070001 function [output] = h2asq_fd(IFO,gpsSec,input,fs,sigType) 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/DARM_ERR counts. 0005 % 0006 % output = h2asq_fd(IFO,gpsSec,input,fs[,sigType]) 0007 % 0008 % IFO - interferometer 0009 % gpsSec - GPS time 0010 % input Input timeseries data in strain. 0011 % fs sampling rate of input timeseries. 0012 % sigType [OPTIONAL] signal type (AS_Q,DARM_ERR) 0013 % <default = AS_Q> 0014 % 0015 % output Output timeseries data in gw channel counts. 0016 % 0017 % Notes: 0018 % 1. Linear interpolation is used to find the response function at 0019 % frequencies desired. At some frequencies (particularly low) 0020 % the magnitude of the interpolated reponse is too low. This could 0021 % be fixed using a higher-order interpolation (spline). However, 0022 % since the reference calibrations are normally supplied at sub-1Hz 0023 % resolutions, I don't see the need to do anything here. 0024 % 2. Any DC data in the reference calibrations is ignored (contains 0025 % NaNs). 0026 % 3. Acausal filtering. Get ringing before start time of an impulsive 0027 % input. The peak time is after the impulse time, though. 0028 % 0029 % $Id: h2asq_fd.m,v 1.8 2005/04/07 01:01:55 thorne Exp $ 0030 % 2003/11/07 0031 0032 %----- Derived quantities. 0033 %---------- Number of data points 0034 N = length(input); 0035 %---------- Frequency resolution 0036 dF = fs/N; 0037 %---------- Vector holding frequencies in usual screwy FFT order: 0038 % vector element: [ 1 2 ... N/2-1 N/2 N/2+1 N/2+2 ... N-1 N ] 0039 % frequency (dF): [ 0 1 ... N/2-2 N/2-1 (N/2 or -N/2) -N/2+1 ... -2 -1 ] 0040 F = [ 0:N/2 , -N/2+1:-1 ]'*dF; 0041 %---------- Fourier transform 0042 input_f = fft(input); 0043 0044 % GET alpha,beta values for this time 0045 % AVERAGE these over the interval 0046 % - since alpha-beta values are updated every 60 sec, 0047 % there is likely only one time used 0048 % GET open loop gain as a function of frequency 0049 % GET sensing function as a function of frequency 0050 % IF frequency steps do not match 0051 % CREATE error message 0052 % ENDIF 0053 output = []; 0054 DEF_SIG = 'AS_Q'; 0055 if(nargin > 4) 0056 if( (strcmp(sigType,'AS_Q')==false) && ... 0057 (strcmp(sigType,'DARM_ERR') == false)) 0058 msgId = 'h2asq_fd:badSignal'; 0059 warning(msgId,'%s: signal type %s is not valid - assume AS_Q',msgId,sigType); 0060 sigType = DEF_SIG; 0061 end 0062 else 0063 sigType = DEF_SIG; 0064 end 0065 [gpsList,alpha_t,beta_t] = getalphabeta(IFO,gpsSec,1); 0066 alpha = mean(alpha_t); 0067 beta = mean(beta_t); 0068 [fH,H] = getopenloopgain(IFO,gpsSec); 0069 [fC,C] = getsensefunction(IFO,gpsSec,sigType); 0070 if (~isequal(fH,fC)) 0071 error('The open-loop gain and sensing function must be sampled at the same frequencies'); 0072 end 0073 0074 %----- Often reference cal files have an f=0 line, which will contain NaNs. 0075 if (fH(1)==0) 0076 fH(1) = []; 0077 H(1) = []; 0078 fC(1) = []; 0079 C(1) = []; 0080 end 0081 0082 % Choose response function from signal type 0083 switch sigType 0084 case 'AS_Q' 0085 %----- Response function (AS_Q/strain) 0086 Rtemp = alpha*C./(1+alpha*beta*H); 0087 case 'DARM_ERR' 0088 %----- Response function (DARM_ERR/strain) 0089 Rtemp = alpha*beta*C./(1+alpha*beta*H); 0090 otherwise 0091 %----- Response function (AS_Q/strain) 0092 Rtemp = alpha*C./(1+alpha*beta*H); 0093 end 0094 0095 %----- fH typically doesn't cover full range of frequencies. Select subset 0096 % of frequencies covered by the reference calibrations and interpolate 0097 % to them. 0098 k = find(F>=fH(1) & F<=fH(end)); 0099 Fcal = F(k); 0100 Rcal = interp1(fH,Rtemp,Fcal); 0101 %Rcal = interp1(fH,Rtemp,Fcal,'spline'); 0102 %figure; loglog(fH,abs(Rtemp),Fcal,abs(Rcal),'+'); 0103 %----- Copy interpolated response function into new holder of same size as 0104 % data. Leave response function as zero at frequencies for 0105 % which we have no calibration data. 0106 R = zeros(size(F)); 0107 R(k) = Rcal; 0108 %----- Sharp drop in response function produces ringing in 0109 % time domain at f ~ f_drop (around 4kHz). Found can eliminate by 0110 % making response function go smoothly to zero at DC and Nyquist 0111 % frequencies. (Smoothing at DC seems to have no effect, but do 0112 % anyway.) 0113 %---------- DC side. 0114 slope = R(k(1))/(k(1)-1); 0115 for j=2:k(1)-1 0116 R(j) = slope*(j-1); 0117 end 0118 %---------- Nyquist side. 0119 % ** NOTE: we only need to make this correction if our frequency range 0120 % is beyond the Nyquist of the calibration data. 0121 % Otherwise, we get a nasty divide-by-zero 0122 if( k(end) < ((N/2)+1)) 0123 slope = -R(k(end))/(N/2+1-k(end)); 0124 for j=k(end)+1:N/2+1 0125 R(j) = R(k(end))+slope*(j-k(end)); 0126 end 0127 end 0128 %----- Set R(-f) = R*(f) (real filter). 0129 R(N/2+2:N) = conj(R(N/2:-1:2)); 0130 %----- Plot original and inrepolated data. 0131 %figure; 0132 %subplot(2,1,1) 0133 %plot(fH,abs(Rtemp),'+',F,abs(R),'markersize',4); 0134 %grid 0135 %subplot(2,1,2) 0136 %plot(fH,angle(Rtemp),'+',F,angle(R),'markersize',4); 0137 %plot(fH,unwrap(angle(Rtemp)),'+',F,unwrap(angle(R)),'markersize',4); 0138 %grid 0139 0140 %----- Calibrate input data in frequency domain. 0141 output_f = input_f.*R; 0142 output = ifft(output_f); 0143 %---------- Cast off any imaginary part (introduced through roundoff errors). 0144 output = real(output); 0145 0146 %----- Plot powers. 0147 %[Pinput,w] = pwelch(input); 0148 %[Poutput,w] = pwelch(output); 0149 %R_chk = interp1(fH,Rtemp,w*fs/2/pi); 0150 %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)); 0151 %grid 0152 0153 return