%-- Sign/delay/amplitude check

%-- Change this to 1 to make plots, or 0 to skip that
makeplots = 1;

%-- Filter parameters
ffreqlo = [  70, 70, 70 ];
ffreqhi = [ 260, 260, 260 ];
wfiles = { 'coherenttest1from15hz_1126257408.out', 'coherenttest1from15hz_1126257408.out', 'coherenttest1from15hz_1126257408.out' };
sigoffset = 8;
%-- Note: 'gpsstarts' times should be 8 seconds before the signal time
gpsstarts = [ 1128303091, 1128303224, 1128303391 ];
injtype = { 'PCAL', 'PCAL', 'PCAL' };
dur = [ 0.04, 0.04, 0.04 ];
hfiles = { '/archive/frames/O1/hoft/H1/H-H1_HOFT_C00-11283/H-H1_HOFT_C00-1128300544-4096.gwf' '/archive/frames/O1/hoft/H1/H-H1_HOFT_C00-11283/H-H1_HOFT_C00-1128300544-4096.gwf' '/archive/frames/O1/hoft/H1/H-H1_HOFT_C00-11283/H-H1_HOFT_C00-1128300544-4096.gwf' };

for ifile = 1:length(wfiles)

%-- Construct butterworth filter
[B,A] = butter(3,[ffreqlo(ifile),ffreqhi(ifile)]/8192);

%-- Read intended strain signal from waveform file, and filter it
model = load(wfiles{ifile})';
nsecs = length(model)/16384;
predicted = filtfilt(B,A,model);
reltimes = [(-sigoffset*16384) : ((nsecs-sigoffset)*16384-1) ] / 16384;

gpsstart = gpsstarts(ifile);

%-- Read in the right length of H1 data, and filter it
htime = gpsstart;
h1data = frgetvect(hfiles{ifile},'H1:GDS-CALIB_STRAIN',htime,nsecs)';
h1f = filtfilt(B,A,h1data);

%-- Select an appropriate sub-interval for calculating the corrrelation
sel = find( reltimes > -dur(ifile) & reltimes < -0.25*dur(ifile) );
showtimes = reltimes(sel);

%-- Find out what time delay makes the predicted curve match the data best
ihbest=0; hbest=0;
shifted = predicted;
for idelay=[0:100]
  hcorr = dot( shifted(sel), h1f(sel) );
  if abs(hcorr) > abs(hbest)
    ihbest = idelay; hbest = hcorr;
    hscale = sqrt( sum(h1f(sel).^2) / sum(shifted(sel).^2) );
  end

  shifted = [0,shifted([1:nsecs*16384-1])];  
end

%-- Select an appropriate sub-interval for plotting
psel = find( reltimes > -2*dur(ifile) & reltimes < dur(ifile) );
showtimes = reltimes(psel);

%-- Now plot
if makeplots == 1

h=plot(showtimes,h1f(psel),'r','LineWidth',2);
hold on
plot(showtimes,predicted(psel),'b')
hold off
set(gca,'XLim',[min(showtimes),max(showtimes)])
maxampl = max( max(abs(h1f(psel))), max(abs(predicted(psel))) );
set(gca,'YLim',[-1.3*maxampl,1.1*maxampl]);
title('H1 h(t) (red) and intended strain (blue) after bandpass')
xlabel(sprintf('time relative to %10d (s)',gpsstart+sigoffset));
ylabel('Absolute amplitude');
xlims = get(gca,'XLim'); ylims = get(gca,'YLim');

hdelay = ihbest/16.384;
stext = sprintf('Best match: scale model by %.2f with %d sample (%.2f ms) delay',hscale,ihbest,hdelay);
if hbest > 0
  text(0.95*xlims(1),0.85*ylims(1),stext);
else
  text(0.95*xlims(1),0.85*ylims(1),[stext,', inverted'])
end

h = text( 0.95*xlims(1), 0.8*ylims(2), injtype{ifile} );
set(h,'FontSize',24,'FontWeight','bold');

print('-dpng','-r200',sprintf('PCALtest_%02d.png',ifile));
ignore = input('Press RETURN to continue');

end  %-- end if makeplots

end  %-- end loop over waveform files
