%-- Sign/delay/amplitude check

%-- Change this to 1 to make plots for each frequency, or 0 to skip that
makeplots = 1

%-- Filter parameters
ffreqlo = [  50, 125, 17, 350 ];
ffreqhi = [ 200, 475, 24, 475 ];
wfiles = { 'hwinj_1126270499_2_H1.txt', 'hwinj_1126280499_2_H1.txt', 'hwinj_1126300499_2_H1.txt', 'hwinj_1126310499_2_H1.txt' };
halfdur = [ 0.1, 0.1, 0.15, 0.1 ];
hfiles = { '/archive/frames/ER8/hoft/H1/H-H1_HOFT_C00-11262/H-H1_HOFT_C00-1126268928-4096.gwf', '/archive/frames/ER8/hoft/H1/H-H1_HOFT_C00-11262/H-H1_HOFT_C00-1126277120-4096.gwf', '/archive/frames/ER8/hoft/H1/H-H1_HOFT_C00-11262/H-H1_HOFT_C00-1126297600-4096.gwf', '/archive/frames/ER8/hoft/H1/H-H1_HOFT_C00-11263/H-H1_HOFT_C00-1126309888-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 = [(-1*16384) : (1*16384-1) ] / 16384;
model = load(wfiles{ifile})';
predicted = filtfilt(B,A,model);

%-- Infer GPS time of start of waveform from the filename
gpsstart = sscanf(wfiles{ifile},'hwinj_%d_');

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

%-- Select an appropriate sub-interval for calculating the corrrelation
sel = find( abs(reltimes) < halfdur(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:32767])];  
end

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

%-- Now plot
if makeplots == 1

%%subplot(2,1,1)
subplot(1,1,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+1));
ylabel('Absolute amplitude');
xlims = get(gca,'XLim'); ylims = get(gca,'YLim');

%%text(0.95*xlims(1),0.8*ylims(2),sprintf('%.1f Hz',f0));
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

if 0
subplot(2,1,2)
h=plot(showtimes,h1data(psel),'Color',[0,.9,0],'LineWidth',2);
hold on
plot(showtimes,model(psel),'b')
hold off
set(gca,'XLim',[min(showtimes),max(showtimes)])
maxampl = max( max(abs(h1data(psel))), max(abs(model(psel))) );
set(gca,'YLim',[-1.3*maxampl,1.1*maxampl]);
title('Raw H1 h(t) data (green) and intended strain (blue)')
xlabel('time (s)'); ylabel('Absolute amplitude');
xlims = get(gca,'XLim'); ylims = get(gca,'YLim');
end %%- end 'if 0'

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

end  %-- end if makeplots

end  %-- end loop over waveform files
