%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [d, scaleVec] = load_template_corr(d,templateFile) Subtracts a provided template from the data (either filtered or classic mode) in order to correct for cold load oscillations. I/O: - 'd' is the input/output data structure. - 'alphaType' (0/1/2) is the type of alpha correction applied. (FILTERED/CLASSIC/POLONLY) - 'scaleType' (0/1/2) is the type of amplitude scaling. Only '0' is currently implemented. MAS -- 29-March-2012 Updates: 01-May-2012 (MAS): Changed to work with latest ASCII template database. Also, removed optional scaling parameter, as gain scaling didn't work. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 function [d, scaleVec] = load_template_corr(d,alphaType,scaleType) 0002 0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0004 % 0005 % function [d, scaleVec] = load_template_corr(d,templateFile) 0006 % 0007 % Subtracts a provided template from the data (either filtered or classic 0008 % mode) in order to correct for cold load oscillations. 0009 % 0010 % I/O: 0011 % - 'd' is the input/output data structure. 0012 % - 'alphaType' (0/1/2) is the type of alpha correction applied. 0013 % (FILTERED/CLASSIC/POLONLY) 0014 % - 'scaleType' (0/1/2) is the type of amplitude scaling. Only '0' 0015 % is currently implemented. 0016 % 0017 % MAS -- 29-March-2012 0018 % 0019 % Updates: 0020 % 01-May-2012 (MAS): Changed to work with latest ASCII template 0021 % database. Also, removed optional scaling 0022 % parameter, as gain scaling didn't work. 0023 % 0024 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0025 0026 if nargin < 3 0027 scaleType = 0; 0028 end 0029 0030 0031 % Read in the template database 0032 [mjdList templates] = readTemplateDatabase(); 0033 0034 % Choose the template MJD. 0035 mjdList(1) = -Inf; 0036 mjdList(end+1) = Inf; 0037 %[~, mjdIndex] = max(histc(d.antenna0.receiver.utc,mjdList)); 0038 [aaa, mjdIndex] = max(histc(d.antenna0.receiver.utc,mjdList)); % me fiddling to get it to work 0039 0040 % Which templates to use... 0041 switch alphaType 0042 case 1 0043 alphaIndex = [17 18 10:15 19 20]; % CLASSIC 0044 case 2 0045 alphaIndex = 9:16; % POLONLY. 0046 otherwise 0047 alphaIndex = 1:8; % FILTERED 0048 end 0049 0050 0051 nPhase = size(templates,1); 0052 nChannels = length(alphaIndex); 0053 0054 % Grab the actual templates: 0055 thetaVec = ((1:nPhase)-0.5) / nPhase * 2 * pi; 0056 templateArray = squeeze(templates(:,mjdIndex,alphaIndex)); 0057 0058 % figure; 0059 % plot(templateArray); 0060 0061 % Now grab the relative utc (subtracting the minimum receiver.utc, which is 0062 % necessary for calculating theta from B and C). 0063 utcVec = d.antenna0.receiver.utc - min(d.antenna0.receiver.utc); 0064 0065 % How much data do we have? 0066 dataLen = length(d.antenna0.receiver.data); 0067 0068 0069 % This is a major part. Pay attention. 0070 % We calculate A, B, and C from the thermal.ccTemperatureLoad vector and 0071 % interpolate and smooth to match the 100 Hz sampling. 0072 [A B C] = load_measure(d, 0.5, 5); 0073 0074 % figure; 0075 % plot(detrend(d.antenna0.thermal.coldLoad)); 0076 % hold on; 0077 % plot(A.*sin(2*pi*B .* utcVec * 86400 + C),'r'); 0078 % input('adfas'); 0079 0080 0081 % Take the B and C vectors and calculate the theta of the load for the 0082 % 100 Hz sampling rate. 0083 thetaLoad = mod(2*pi*B .* utcVec * 86400 + C, 2*pi); 0084 0085 0086 % Now create the full-length templates with proper phase and sampling for 0087 % appropriate number of channels. 0088 template = zeros(dataLen,nChannels); 0089 for k=1:nChannels 0090 template(:,k) = A .* ... 0091 interp1(thetaVec,templateArray(:,k),thetaLoad,'nearest','extrap'); 0092 end 0093 0094 % plot(d.antenna0.receiver.data(:,1) - median(d.antenna0.receiver.data(:,1))); 0095 % hold on; 0096 % plot(template(:,1),'r'); 0097 0098 0099 % Calculate and apply the scaling factor for the templates. 0100 scaleVec = load_scale(d, template, scaleType); 0101 for k=1:nChannels 0102 template(:,k) = template(:,k) * scaleVec(k); 0103 end 0104 0105 % disp(scaleVec); 0106 % plot(template(:,1),'g'); 0107 % input('adfasdf'); 0108 0109 % for k=1:10 0110 % figure; 0111 % plot(detrend(d.antenna0.receiver.data(:,k))); 0112 % hold on; 0113 % plot(detrend(d.antenna0.receiver.data(:,k)) - template(:,k),'r'); 0114 % plot(template(:,k),'g'); 0115 % input('adsfasf'); 0116 % end 0117 0118 % Now do the actual subtraction of the full-length templates from the data. 0119 d.antenna0.receiver.data = d.antenna0.receiver.data - template; 0120 0121 0122 % And we're done! 0123 end