%-- Sign/delay/amplitude check

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

%-- Filter parameters
ffreqlo = [  70, 70, 70, 70 ];
ffreqhi = [ 260, 260, 260, 260 ];
wfiles = { 'H1-HWINJ_CBC-1126257410-12.txt', 'H1-HWINJ_CBC-1126257410-12.txt', 'H1-HWINJ_CBC-1126257410-12.txt', 'H1-HWINJ_CBC-1126257410-12.txt' };
gpsstarts = [ 1127074635, 1127074768, 1127075230, 1127075737 ];
dur = [ 0.04, 0.04, 0.04, 0.04 ];
hfiles = { '/archive/frames/O1/hoft/H1/H-H1_HOFT_C00-11270/H-H1_HOFT_C00-1127071744-4096.gwf', '/archive/frames/O1/hoft/H1/H-H1_HOFT_C00-11270/H-H1_HOFT_C00-1127071744-4096.gwf', '/archive/frames/O1/hoft/H1/H-H1_HOFT_C00-11270/H-H1_HOFT_C00-1127071744-4096.gwf', '/archive/frames/O1/hoft/H1/H-H1_HOFT_C00-11270/H-H1_HOFT_C00-1127071744-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
reltimes = [(-6*16384) : (6*16384-1) ] / 16384;
model = load(wfiles{ifile})';
predicted = filtfilt(B,A,model);

gpsstart = gpsstarts(ifile);

%-- Read in 12 seconds of H1 data, and filter it
htime = gpsstart;
h1data = frgetvect(hfiles{ifile},'H1:GDS-CALIB_STRAIN',htime,12)';
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:12*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+6));
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

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

end  %-- end if makeplots

end  %-- end loop over waveform files
