Description of plot1overf

0001 function stat = plot1overf2(d, plotparams, field)
0002 %  function stat = plot1overf2(d, plotparams, field)
0003 %
0004 %   SJCM's adaptation of joe's makeWebpageSpectra.m.  Main differences are:
0005 %
0006 %  1.  flags data
0007 %  2.  selects scans to plot depending on how much data is unflagged --
0008 %  interpolates over flagged data if just a few
0009 %  3.  takes average spectra of all plots in given track.
0010 %  4.  plots only the average spectra
0011 %  5.  uses all channels, instead of limiting it to 6.
0012 %  6.  all 1/f plots on a single plot, instead of nchan plots
0013 %  7.  writes all fit data out into two files, one for everyscan and one for
0014 %      composite
0015 %
0016 %
0017 %
0018 % Make noise spectra plots for the data.
0019 %
0020 % We make plots of both the noise stares and the general data.
0021 % Inputs:
0022 %   d                     data
0023 
0024 % define variables for sampling
0025 sample_frequency = 100;
0026 samples_per_hour = sample_frequency * 60 * 60;
0027 hours_per_stare_plot = 1;
0028 hours_per_scan_plot = 0.5;
0029 min_stare_length = 200;
0030 min_scan_length = 10000;
0031 samples_between_scan_plots = samples_per_hour*hours_per_scan_plot;
0032 samples_between_stare_plots = samples_per_hour*hours_per_stare_plot;
0033 
0034 % plot the scans
0035 stareEndPoints = chooseEndPoints(d.index.noise_event.fast & ...
0036     (~d.index.noise.fast),min_stare_length,samples_between_stare_plots );
0037 if(plotparams.plot)
0038   figure(1);clf;
0039   setwinsize(gcf, 1000, 750);
0040 else
0041   setPlotDisplay(0);
0042   setwinsize(gcf, 1000, 750);
0043 end
0044 [stat.scan.fknee stat.scan.alpha stat.scan.white stat.scan.mjd stat.scan.errorFlags] = getStatsPlots(d,stareEndPoints, ...
0045       plotparams.save, field, 1);
0046 
0047   
0048 % plot the noise data
0049 scanEndPoints = chooseEndPoints(~d.index.noise_event.fast,min_scan_length,samples_between_scan_plots ...
0050     );
0051 if(plotparams.plot)
0052   figure(2);clf;
0053   setwinsize(gcf, 1000, 750);
0054 else
0055   setPlotDisplay(plotparams.plot);
0056   setwinsize(gcf, 1000, 750);
0057 end
0058 
0059 [stat.noise.fknee stat.noise.alpha stat.noise.white stat.noise.mjd stat.noise.errorFlags] = getStatsPlots(d,stareEndPoints, ...
0060       plotparams.save, field, 2);
0061   
0062 return;
0063 
0064 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0065 
0066 
0067 function [endPoints doPlots] = chooseEndPoints(index,min_length,samples_between_plots)
0068 
0069 % Find the start and end points of the scans.
0070 endPoints=flagEndpoints(index);
0071 %Cut to use only one long enough
0072 lengths=endPoints(:,2)-endPoints(:,1);
0073 longEnough=(lengths>min_length);
0074 endPoints=endPoints(longEnough,:);
0075 nScan = size(endPoints,1);
0076 
0077 return;
0078 
0079 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0080 
0081 function [noise_fknee noise_alpha noise_white mjd errorFlags] = getStatsPlots(d,stareEndPoints,gen, field, figNum)
0082 
0083 
0084 % we don't need to be calculating all the channels.
0085 switch size(d.antenna0.receiver.data,2)
0086   case 6
0087     chansToPlot = [1 6 2 3];
0088     
0089   case 10
0090     chansToPlot = [1 9 7 8];
0091     
0092   case 8
0093     chansToPlot = [1 8 6 7];
0094 end
0095 
0096 nChannel = length(chansToPlot);
0097 grid = [nChannel/2 2];  % subplot grid for this plot.
0098 sample_frequency = 100;
0099 
0100 % we'll cut nStare down depending on how well they do.
0101 nStare=size(stareEndPoints,1);
0102 noise_white=zeros(nStare,nChannel);
0103 noise_alpha=zeros(nStare,nChannel);
0104 noise_fknee=zeros(nStare,nChannel);
0105 mjd=zeros(nStare,1);
0106 errorFlags=zeros(nStare,nChannel);
0107 
0108 % apply the flags
0109 d = applyFlags(d);
0110 
0111 
0112 for stareIndex=1:nStare,
0113   
0114   stare_start_index=stareEndPoints(stareIndex,1);
0115   stare_end_index=stareEndPoints(stareIndex,2);
0116   
0117   stareRange1=stare_start_index+1:stare_end_index-1;
0118   
0119   start_utc = d.antenna0.receiver.utc(stare_start_index);
0120   start_string = mjd2string(start_utc);
0121   mjd(stareIndex) = start_utc;
0122 
0123   for channel=1:nChannel,
0124     stareData=d.antenna0.receiver.data(stareRange1,chansToPlot(channel));
0125     % check for the length of unflagged-ness
0126     lflag = length(find(isnan(stareData)));
0127     ll    = length(stareData);
0128     if(lflag > 0.03*ll)
0129       % too much flagged data, don't trust
0130       aa = floor((ll-1)/2);
0131       alpha1 = nan;
0132       fknee1 = nan;
0133       white1 = nan;
0134       fn1    = nan(aa,1);;
0135       measuredn1 = nan(aa,1);
0136       fittedn1 = nan(aa,1);
0137       errorFlag1 = nan;
0138     else
0139       % interpolate so we have continuous data stream
0140       goodStareData = stareData;
0141       goodLength = 1:ll;
0142       goodLength(isnan(stareData)) = [];
0143       goodStareData(isnan(stareData)) = [];
0144       stareData = interp1(goodLength, goodStareData, 1:ll, 'linear', ...
0145       'extrap');
0146       [alpha1 fknee1 white1 fn1 measuredn1 fittedn1 errorFlag1] = fitOneOverF(sample_frequency,stareData');
0147     end
0148     
0149     noise_white(stareIndex, channel) = white1;
0150     noise_alpha(stareIndex, channel) = alpha1;
0151     noise_fknee(stareIndex, channel) = fknee1;
0152     errorFlags(stareIndex, channel) = errorFlag1;
0153     fn1_all{stareIndex, channel} = fn1;
0154     meas1_all{stareIndex, channel} = measuredn1;
0155     fit1_all{stareIndex, channel} = fittedn1;
0156   end
0157 end
0158 
0159 % next we make our plot, first we need to get an average fn1 on which we
0160 % sample measuredn1 and fittedn1
0161 thisSize = [];
0162 for m=1:nChannel
0163   thisSize = [];
0164   for mm=1:size(fn1_all,1)
0165     thisSize = [thisSize; length(fn1_all{mm,m})];
0166   end
0167   f = find(thisSize == max(thisSize),1);
0168   fn1 = fn1_all{f,m};
0169   allMes = [];
0170   allFit = [];
0171   for mm=1:size(fn1_all,1)
0172     if(~isnan(fn1_all{mm,m}))
0173       if(length(fn1) == length(fn1_all{mm,m}))
0174     % no interpolation needed.
0175     allMeas(mm,:) = meas1_all{mm,m};
0176     allFit(mm,:)  = fit1_all{mm,m};
0177       else
0178     % interpolate on same data scale.
0179     thisMeas = interp1(fn1_all{mm,m}, meas1_all{mm,m}, fn1, 'linear', ...
0180         'extrap');
0181     thisFit = interp1(fn1_all{mm,m}, fit1_all{mm,m}, fn1, 'linear', ...
0182         'extrap');      
0183     allMeas(mm,:) = thisMeas;
0184     allFit(mm,:) = thisFit;
0185       end
0186     else
0187       allMeas(mm,:) = nan(size(fn1));
0188       allFit(mm,:)  = nan(size(fn1));
0189     end
0190   end
0191 
0192   % next we take an average
0193   %thisMeas = nanmean(allMeas, 2);
0194   %thisStd  = nanstd(allMeas, [], 2);
0195   %thisFit = nanmean(allFit,2);
0196   mean_white = nanmean(noise_white(:,m));
0197   mean_fknee = nanmean(noise_fknee(:,m));
0198   mean_alpha = nanmean(noise_alpha(:,m));
0199   
0200   
0201   % let's not take the average, let's just plot them all.
0202   subplot(grid(1), grid(2), m);
0203   loglog(fn1, allMeas);
0204   if(m==grid(1) || m==nChannel)
0205     xlabel('Frequency [Hz]');
0206   end
0207   ylabel('Power');
0208   switch(m)
0209     case 1
0210       title(sprintf('I1 [%e, %e, %e]',mean_alpha, mean_fknee, mean_white));
0211       
0212     case 2
0213       title(sprintf('I2 [%e, %e, %e]',mean_alpha, mean_fknee, mean_white));
0214 
0215     case 3
0216       title(sprintf('Q [%e, %e, %e]',mean_alpha, mean_fknee, mean_white));      
0217 
0218     case 4
0219       title(sprintf('U [%e, %e, %e]',mean_alpha, mean_fknee, mean_white));      
0220   end
0221 end
0222 
0223 numGoodScans = ~isnan(noise_fknee);
0224 numGoodScans = sum(numGoodScans,1);
0225 numGoodScans = median(numGoodScans);
0226 
0227 
0228 switch(figNum)
0229   case 1
0230     gtitle(sprintf('Power Spectra for %d Scans', numGoodScans));
0231     
0232   case 2
0233     gtitle(sprintf('Power Spectra for %d Noise Source Events', numGoodScans));
0234 end
0235 
0236     
0237 % if we're generating a plot, do so
0238 if(gen)
0239   % get the directory to write to
0240   maindir=getMainDir(d, field);
0241   dbclear if error
0242   set(gcf,'paperposition',[0 0 6.0 6.0])
0243   eval(sprintf('print -dpng -r200 %s/overf/fig%d.png;', ...
0244       maindir, figNum));
0245   dbstop if error
0246 end
0247 
0248 return;
0249
plot1overf

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
