analyzeSpectrograms

 This script reads in data between start and finish. It plots and returns
 spectrogram and some telescope diagnostics.

 ogk, 21 Oct 2010

0001 function [time,freq,spectrogram,diagnostics] = analyzeSpectrograms(start,finish,chan)
0002 % This script reads in data between start and finish. It plots and returns
0003 % spectrogram and some telescope diagnostics.
0004 %
0005 % ogk, 21 Oct 2010
0006 
0007 % Parse start and finish:
0008 mjdstart = tstr2mjd(start);
0009 mjdstop  = tstr2mjd(finish);
0010 % to convert mjd to string: mjd2string(mjdtime)
0011 
0012 % Take the data in 2 hour segments:
0013 Nreps = floor((mjdstop-mjdstart)*24/2);
0014 
0015 if Nreps == 0
0016     Nreps = 1;
0017 end
0018 
0019 fprintf('Number of 2 hour segments to scan: %d\n',Nreps)
0020 
0021 
0022 time = [];
0023 freq = [];
0024 spectrogram = [];
0025 diagnostics = [];
0026 
0027 for k=1:Nreps
0028     stime = mjdstart + (k-1)*2/24;
0029     if k == Nreps
0030         ftime = mjdstop;
0031     else
0032         ftime = mjdstart + k*2/24;
0033     end
0034     fprintf('***** Reading scan %d of %d:\n',k,Nreps)
0035     disp(['      Start:  ' mjd2string(stime)])
0036     disp(['      Finish: ' mjd2string(ftime)])
0037      
0038     % Read in the data:
0039     try
0040         d = pipe_read(mjd2string(stime),mjd2string(ftime));   
0041         Nfft = 1024;
0042         dat = d.antenna0.receiver.switchData(:,chan);
0043         t = d.antenna0.receiver.utc;
0044         N = floor(size(dat,1)/Nfft); % number of spectra to get
0045         [freq,Y] = ffto(reshape(dat(1:N*Nfft),Nfft,N)-...
0046             repmat(mean(reshape(dat(1:N*Nfft),Nfft,N),1),Nfft,1),100);
0047         spectrogram = [spectrogram Y];
0048         time = [time mean(reshape(t(1:N*Nfft),Nfft,N),1)];
0049         ddiag = d.antenna0.thermal.dlpTemperatureSensors(:,1:7);
0050         ddiag = interp1(d.array.frame.utc,ddiag,d.antenna0.receiver.utc,'linear');
0051         ddiag = ddiag(1:N*Nfft,:);
0052         diagnostics = [diagnostics transpose(reshape(mean(reshape(ddiag,Nfft,N*7),1),N,7))];
0053         
0054     catch
0055         disp('***** BAD DATA; SKIPPING TO NEXT SET')
0056     end
0057 end
0058 
0059 figure
0060 subplot(2,1,1)
0061 imagesc(time-time(1),freq(freq>0),log10(abs(spectrogram(freq>0,:))),[-0.5 2])
0062 xlabel('Time [days]')
0063 ylabel('Frequency [Hz]')
0064 colormap('gray')
0065 subplot(2,1,2)
0066 % plot(time-time(1),diagnostics)
0067 % xlabel('Time [days]')
0068 % ylabel('DLP Temperature Sensors')
0069 
0070 plot(time-time(1),sum(abs(spectrogram(freq > 39 & freq < 41,:)),1),'k-',...
0071     time-time(1),sum(abs(spectrogram(freq > 35 & freq < 37,:)),1),'r-')
0072 xlabel('Time [days]')
0073 ylabel('Power')
0074 legend('39 Hz to 41 Hz','Ref (35 to 37 Hz)')
0075 
0076 end