load_scale

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

 function scaleVec = load_scale(d, templateArray, mode)

   Calculates the scale factors for the data arrays given a set of
   templates.
   
   I/O:
    - 'd' is the data structure to which we wish to apply the template.
    - 'templateArray' is array of N templates.
    - 'scaleType' (0,1,2) is the scaling method 
    (fit to data / fit to switchData / alpha database).  Only type '0' is
    currently implemented.
    - 'scaleVec' is the N scale factors.

   MAS -- 29-March-2012

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

0001 function scaleVec = load_scale(d, templateArray, scaleType)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % function scaleVec = load_scale(d, templateArray, mode)
0006 %
0007 %   Calculates the scale factors for the data arrays given a set of
0008 %   templates.
0009 %
0010 %   I/O:
0011 %    - 'd' is the data structure to which we wish to apply the template.
0012 %    - 'templateArray' is array of N templates.
0013 %    - 'scaleType' (0,1,2) is the scaling method
0014 %    (fit to data / fit to switchData / alpha database).  Only type '0' is
0015 %    currently implemented.
0016 %    - 'scaleVec' is the N scale factors.
0017 %
0018 %   MAS -- 29-March-2012
0019 %
0020 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0021 
0022 % XXX This is on MAS's todo list.
0023 if scaleType ~= 0
0024     disp('Terribly sorry; only scaling by direct fitting is currently supported.');
0025 end
0026 
0027 
0028 
0029 
0030 % Parameters for ignoring data based on standard deviation.
0031 tScale = 10;    % Seconds
0032 nClip = 5;      % Number of iterations of data rejection.
0033 clipThresh = 3; % Sigma threshold for data rejection.
0034 
0035 
0036 % Length of the std deviation rejection in 10ms samples.
0037 tN = tScale * 100;
0038 
0039 
0040 % How many channels are there, by the way?
0041 nChannels = size(templateArray,2);
0042 
0043 
0044 % Read in the filter...
0045 [home,installeddir] = where_am_i();
0046 f = open([home '/' installeddir '/reduc/load/filters/load_filter_hp.mat']);
0047 
0048 
0049 % Removing noise diode events...
0050 d = noiseRemove(d);
0051 
0052 
0053 
0054 % Filtering the data below 1.1Hz...
0055 filtData = filter(f.Num, 1, d.antenna0.receiver.data);
0056 
0057 filtLen = length(f.Num);
0058 filtData = filtData(filtLen:end,:);
0059 templateArray = templateArray(ceil(filtLen/2):end-floor(filtLen/2),:);
0060 
0061 utcVec = d.antenna0.receiver.utc(ceil(filtLen/2):end-floor(filtLen/2));
0062 
0063 
0064 % Rejecting data based on standard deviations...
0065 dataLen = length(filtData);
0066 
0067 nSub = floor(dataLen / tN);
0068 subStart = floor((dataLen - nSub * tN) / 2) + 1;
0069 subEnd = dataLen - ceil((dataLen - nSub * tN) / 2);
0070 
0071 % disp(size(filtData(subStart:subEnd,:)));
0072 % disp(tN);
0073 % disp(nSub);
0074 % disp(nChannels);
0075 stdData = squeeze(std( ...
0076     reshape(filtData(subStart:subEnd,:), tN, nSub, nChannels)));
0077 stdUtc = mean(reshape(utcVec(subStart:subEnd), tN, nSub))';
0078 
0079 
0080 % Initialize the clip vector.
0081 clipVec = zeros(nSub,1);
0082 
0083 % Do the clipping...
0084 for l=1:nClip
0085     stdDataStd = std(stdData(clipVec == 0,:));
0086     stdDataMed = median(stdData(clipVec == 0,:));
0087 
0088     clipVec = max(abs(stdData - repmat(stdDataMed,nSub,1)) > ...
0089         clipThresh * repmat(stdDataStd,nSub,1),[],2);
0090 end
0091 
0092 % Lengthen the flagging vector.  Also, rejection via proximity.
0093 flagVec = interp1(stdUtc, +(smooth(clipVec,2) > 0), ...
0094     utcVec, 'nearest', 'extrap');
0095 
0096 
0097 % Apply the flags...
0098 filtData = filtData(flagVec == 0,:);
0099 templateArray = templateArray(flagVec == 0,:);
0100 dfLen = sum(flagVec == 0);
0101 
0102 % Now calculate the scaling factors for each channel.
0103 scaleVec = zeros(1,nChannels);
0104 for k=1:nChannels
0105 %     figure;
0106 %     semilogy(abs(fft(filtData(1:2^12,k))).^2);
0107 %     hold on;
0108     fitResults = [ones(dfLen,1) templateArray(:,k)] \ filtData(:,k);
0109 %     semilogy(abs(fft(fitResults(2) * templateArray(1:2^12,k))).^2,'r');
0110 %     disp(fitResults);
0111 %     input('adfasd');
0112     
0113     scaleVec(k) = fitResults(2);
0114 end
0115     
0116 end
0117