calculate_mains

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

 function [mains_factors indices] = calculate_mains(d, useSwitch)

   To apply a correction factor to remove the effect
   of fluctuations picked up from mains at 60Hz sin-wave 
   sampled at 100Hz --> 20Hz saw tooth like effect
   Do this on each of the 6 channels

   switch set to 1 means you calculate on teh switch data.
    act

    mod by sjcm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [mains_factors indices] = calculate_mains(d, useSwitch)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % function [mains_factors indices] = calculate_mains(d, useSwitch)
0006 %
0007 %   To apply a correction factor to remove the effect
0008 %   of fluctuations picked up from mains at 60Hz sin-wave
0009 %   sampled at 100Hz --> 20Hz saw tooth like effect
0010 %   Do this on each of the 6 channels
0011 %
0012 %   switch set to 1 means you calculate on teh switch data.
0013 %    act
0014 %
0015 %    mod by sjcm
0016 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0017 
0018 numSamples = 1000;
0019 maxHarmonic = 10;
0020 if(useSwitch)
0021   chan = 1:24;
0022 else
0023   chan = 1:6;
0024 end
0025 
0026 % check if there are 1000 samples of clean data
0027 if(any(d.index.blank.fast))
0028   % we have a clean set on a blank sky.  let us make sure there are at least
0029   % 1000 unflagged samples
0030   [dblank indOrig]   = cutObs(d, 'blank', 'only');
0031   forig = find(indOrig);
0032   res_utc  = dblank.antenna0.receiver.utc;
0033   if(useSwitch)
0034     res_data = dblank.antenna0.receiver.switchData;
0035     res_data(dblank.flags.fast) = nan;
0036   else
0037     res_data = dblank.antenna0.receiver.data;    
0038     res_data(dblank.flags.fast) = nan;
0039   end
0040 else
0041   % in this case, we do not have a blank sky obs
0042   forig    = find(d.antenna0.receiver.utc);
0043   res_utc  = d.antenna0.receiver.utc;
0044   if(useSwitch)
0045     res_data = d.antenna0.receiver.switchData;
0046     res_data(d.flags.switch) = nan;
0047   else
0048     res_data = d.antenna0.receiver.data;
0049     res_data(d.flags.fast) = nan;
0050   end
0051 end
0052 
0053 
0054 % look for at least 1000 entries that are not flagged in each channel
0055 for mm=1:length(chan)
0056   indFlag = isnan(res_data(:,mm));
0057 
0058   [si(mm) ei(mm)] = findIndices(indFlag);
0059 
0060   numSamples(mm) = ei(mm) - si(mm);
0061 end
0062 
0063   
0064 % indices from the full data set
0065 fsi = forig(si);
0066 fei = forig(ei);
0067 
0068 indices = [fsi fei];
0069 
0070 sampLen = fei - fsi + 1;
0071 
0072 for chanNumber = 1:length(chan) 
0073 
0074   display(sprintf('calculate_mains:: Channel %d - Calculating for %d samples', chanNumber, ...
0075       sampLen(chanNumber)));
0076 
0077   % pick out the sample data for chanNumber
0078   residual_data = res_data(si(chanNumber):ei(chanNumber),chanNumber);
0079   
0080   % Create a template mains signal - 1000 samples
0081   % Make the same length as data
0082   x = (res_utc(si(chanNumber):ei(chanNumber)))*24*60*60;
0083   
0084   for harmonic=1:maxHarmonic
0085     
0086     short_data = residual_data;
0087     short_flags = logical(zeros(size(residual_data)));
0088     
0089     % Now loop over phi and find best fit via correlation (dot_product)
0090     
0091     for j=1:100
0092       freq=(-0.1+(j*0.002));
0093       for i=1:300
0094     phi=(2*pi)*i/300;
0095     y=sin(   (x*2*pi*( (freq+60)*harmonic) )+phi );
0096     product(j,i) = (dot(y(~short_flags),short_data(~short_flags))); % should not be necessary.
0097     phase(i)=phi;
0098     f(j)=(harmonic*(freq+60));
0099       end
0100     end
0101     %Find the maximum combined value of phase and freq
0102     [C,I]=max(max(product));
0103     [L,M]=max(product(:,I));
0104     
0105     best_phase = phase(I);
0106     best_freq  = f(M);
0107     best_corr  = product(M,I);
0108     
0109     best_fit_sin_one = sin((x*2*pi*(best_freq))+best_phase);
0110     
0111     % Now find the best fit amplitude
0112     
0113     short_zero_data = short_data - mean(short_data(~short_flags));
0114     
0115     %First sine - mains at ~60Hz
0116     
0117     amplitude_function = @(A) ( std(short_zero_data(~short_flags) - (A*best_fit_sin_one(~short_flags))));
0118     
0119     
0120     [A,f_val] = fminsearch(amplitude_function, [1 ]);
0121     
0122     best_amp_one = A(1);
0123     
0124     best_fit_sin_one = best_fit_sin_one * best_amp_one;
0125     
0126     residual_data = short_zero_data  - best_fit_sin_one;
0127     
0128     % Now plot the power spectrum of the residuals
0129     residual_data(short_flags)=0;
0130     t_sec = (d.antenna0.receiver.utc(1:1000) - d.antenna0.receiver.utc(1))*24*60*60;
0131     Fs = 1/(t_sec(2)-t_sec(1));
0132     L=length(residual_data);
0133     NFFT = 2^nextpow2(L);
0134     
0135     y_power_spec = (fft(residual_data,NFFT)/L);
0136     x_power_spec = Fs/2*linspace(0,1,NFFT/2+1);
0137     
0138     
0139     if harmonic==1;
0140       short_zero_data(short_flags)=0;
0141       t_sec = (d.antenna0.receiver.utc(1:1000) - d.antenna0.receiver.utc(1))*24*60*60;
0142       Fs = 1/(t_sec(2)-t_sec(1));
0143       L=length(short_zero_data);
0144       NFFT = 2^nextpow2(L);
0145       
0146       y_power_spec = (fft(short_zero_data,NFFT)/L);
0147       x_power_spec = Fs/2*linspace(0,1,NFFT/2+1);
0148     end
0149     
0150     mains_factors(chanNumber,harmonic,:) = [best_freq best_phase best_amp_one];
0151     
0152   end
0153 end
0154 
0155 
0156 return;
0157 
0158 
0159 
0160 
0161 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0162 function [si ei] = findIndices(indFlag)
0163 
0164 %  find the indices where the observations start/stop
0165 stop =0;
0166 index = 1;
0167 indFlag = double(indFlag);
0168 while(stop==0)
0169   s = find(indFlag == 0,1);
0170   
0171   if(isempty(s))
0172     stop = 1;
0173   else
0174     indFlag(1:s) = 2;
0175     e = find(indFlag == 1, 1);
0176     if(isempty(e))
0177       e = length(indFlag);
0178       stop = 1;
0179     end
0180     indFlag(1:e) = 2;
0181 
0182     si(index) = s;
0183     ei(index) = e;
0184     index = index+1;
0185   end
0186 end
0187 
0188 diffVals = ei-si;
0189 f = find(diffVals == max(diffVals));
0190 
0191 
0192 si = si(f);
0193 ei = ei(f);
0194 
0195 if( (ei - si) > 1000)
0196   ei = si+1000-1;
0197 end
0198 
0199 
0200 
0201 return;