genTemplateEstimates

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 function [numVec templateArray] = genTemplateEstimates(utcStart,utcEnd,newEstimatesDir,nTheta)

   Reads in the 2-hour cold-load template estimates from the archive.   


   I/O:

   MAS -- 20-April-2012

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [numVec templateArray] = genTemplateEstimates(utcStart,utcEnd,newEstimatesDir,nTheta)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % function [numVec templateArray] = genTemplateEstimates(utcStart,utcEnd,newEstimatesDir,nTheta)
0006 %
0007 %   Reads in the 2-hour cold-load template estimates from the archive.
0008 %
0009 %
0010 %   I/O:
0011 %
0012 %   MAS -- 20-April-2012
0013 %
0014 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0015 
0016 
0017 % Parameters for rejecting RFI affected / bright sky data.
0018 tScale = 10;    % Seconds
0019 nClip = 5;      % Number of iterations of data rejection.
0020 clipThresh = 3; % Sigma threshold for data rejection.
0021 tN = tScale * 100; % Length of the std deviation rejection in 10ms samples.
0022 
0023 nChannels = 20;
0024 
0025 % Read in the filter...
0026 [home,installeddir] = where_am_i();
0027 f = open([home '/' installeddir '/reduc/load/filters/load_filter_hp.mat']);
0028 
0029 
0030 % Get the data.
0031 disp('Reading in data...');
0032 % I know I'm not supposed to do this, but read_arc can crash when
0033 % there's no data in the archive...
0034 try
0035     d = read_arc(utcStart,utcEnd);
0036 catch
0037     disp('Encountered an error, moving on...');
0038     numVec = -1;
0039     templateArray = -1;
0040     return;
0041 end
0042 % Make sure that we have enough data points...
0043 if length(d.antenna0.receiver.data) < ...
0044         86400 * 100 * (tstr2mjd(utcEnd) - tstr2mjd(utcStart)) / 2
0045     disp('Data appear to have a problem, moving on...');
0046     numVec = -1;
0047     templateArray = -1;
0048     return;
0049 end
0050 
0051 
0052 disp('Applying alpha corrections.')
0053 d_filt = assembleAlphaStreams(d,'FILTERED');
0054 d_filt = applyAlpha(d_filt,'FILTERED');
0055 d_polonly = assembleAlphaStreams(d,'POLONLY');
0056 d_polonly = applyAlpha(d_polonly,'POLONLY');
0057 d_classic = assembleAlphaStreams(d,'CLASSIC');
0058 d_classic = applyAlpha(d_classic,'CLASSIC');
0059 
0060 
0061 
0062 
0063 % Bring noise diode events in line with rest of data so transitions
0064 % don't interfere with later calculations.
0065 disp('Removing noise diode events...');
0066 d_filt = noiseRemove(d_filt);
0067 d_classic = noiseRemove(d_classic);
0068 d_polonly = noiseRemove(d_polonly);
0069 
0070 % Now assemble our channels:
0071 channelData = [d_filt.antenna0.receiver.data, ...
0072     d_polonly.antenna0.receiver.data, ...
0073     d_classic.antenna0.receiver.data(:,[1 2 9 10])];
0074 
0075 clear d_filt;
0076 clear d_classic;
0077 clear d_polonly;
0078 
0079 
0080 % High-pass filter to remove the 1/f noise and most of the sky signal.
0081 disp('Filtering the data below 1.1Hz.');
0082 filtData = filter(f.Num, 1, channelData);
0083 
0084 
0085 filtLen = length(f.Num);
0086 filtData = filtData(filtLen:end,:);
0087 utcVec = ...
0088     d.antenna0.receiver.utc(ceil(filtLen/2):end-floor(filtLen/2)) - ...
0089     min(d.antenna0.receiver.utc);
0090 
0091 
0092 
0093 % Now we ignore any remaining data with anomalously high standard
0094 % deviations.  Causes of this would include RFI or very bright sky
0095 % signal.
0096 disp('Rejecting data based on standard deviations.');
0097 dataLen = length(utcVec);
0098 
0099 nSub = floor(dataLen / tN);
0100 subStart = floor((dataLen - nSub * tN) / 2) + 1;
0101 subEnd = dataLen - ceil((dataLen - nSub * tN) / 2);
0102 
0103 stdData = squeeze(std( ...
0104     reshape(filtData(subStart:subEnd,:), tN, nSub, nChannels)));
0105 stdUtc = mean(reshape(utcVec(subStart:subEnd), tN, nSub))';
0106 
0107 
0108 % Initialize the clip vector.
0109 clipVec = zeros(nSub,1);
0110 
0111 % Do the clipping...
0112 for l=1:nClip
0113     stdDataStd = std(stdData(clipVec == 0,:));
0114     stdDataMed = median(stdData(clipVec == 0,:));
0115 
0116     clipVec = max(abs(stdData - repmat(stdDataMed,nSub,1)) > ...
0117         clipThresh * repmat(stdDataStd,nSub,1),[],2);
0118 end
0119 
0120 % Lengthen the flagging vector.  Also, rejection via proximity.
0121 flagVec = interp1(stdUtc, double(smooth(clipVec,2) > 0), ...
0122     utcVec, 'nearest', 'extrap');
0123 
0124 
0125 % Calculate the A, B, and C vectors from the ccColdPlate vector,
0126 % allowing for calculation of theta(t).
0127 disp('Calculating the load parameters.');
0128 [A B C] = load_measure(d, 0.5, 5);
0129 
0130 % Clip down to match the lengths of the filtered data.
0131 A = A(ceil(filtLen/2):end-floor(filtLen/2));
0132 B = B(ceil(filtLen/2):end-floor(filtLen/2));
0133 C = C(ceil(filtLen/2):end-floor(filtLen/2));
0134 
0135 % Calculate theta(t)
0136 thetaLoad = mod(2*pi*B .* utcVec * 86400 + C, 2*pi);
0137 
0138 
0139 % Remove the mean from the data.
0140 filtData = filtData - repmat(mean(filtData(flagVec == 0,:)),dataLen,1);
0141 % Normalize the data by the A vector.
0142 normData = filtData ./ repmat(A,1,nChannels);
0143 % Also normalize by standard deviation.
0144 normData = normData ./ repmat(std(normData(flagVec == 0,:)),dataLen,1);
0145 
0146 
0147 % Prepare the output arrays.
0148 %thetaVec = 2*pi * ((1:nTheta)'-0.5)/nTheta;
0149 templateArray = zeros(nTheta,nChannels);
0150 numVec = zeros(nTheta,1);
0151 
0152 % This is the binning stage.
0153 disp('Binning the data to create templates.');
0154 % Note NaN values:
0155 notnan = ~isnan(sum(normData,2));
0156 for l=1:nTheta
0157     templateArray(l,:) = nanmean( ...
0158         normData(thetaLoad >= (l-1)*2*pi/nTheta & ...
0159         thetaLoad < l*2*pi/nTheta & flagVec == 0 & ...
0160         notnan,:));
0161 
0162     numVec(l) = sum(thetaLoad >= (l-1)*2*pi/nTheta & ...
0163         thetaLoad < l*2*pi/nTheta & flagVec == 0 & notnan);
0164 end
0165 
0166 
0167 [year, month, day, hour, minute, second] = mjd2date(tstr2mjd(utcStart));
0168 
0169 if (minute > 30)
0170     [year, month, day, hour, minute, second] = mjd2date(tstr2mjd(utcStart)+0.1/12);
0171 end
0172 
0173 
0174 fileName = sprintf('%04d-%02d-%02dT%02d.bin', ...
0175     year, month, day, hour);
0176 
0177 fid = fopen([newEstimatesDir '/' fileName],'wb');
0178 fwrite(fid, [double(numVec), templateArray],'float');
0179 fclose(fid);
0180 
0181 end