% get data
start_time = 1107840506;
end_time = 1107841305;

g = GWData;
[darm, t, info] = g.fetch(start_time, end_time, 'H1:LSC-DARM_IN1_DQ');
fs = 1/t(2);

% aux_channels = {'H1:IMC-DOF_1_P_IN1_DQ', 'H1:IMC-DOF_1_Y_IN1_DQ', ...
%     'H1:IMC-DOF_2_P_IN1_DQ', 'H1:IMC-DOF_2_Y_IN1_DQ', ...
%     'H1:IMC-DOF_3_P_IN1_DQ', 'H1:IMC-DOF_3_Y_IN1_DQ', ...
%     'H1:IMC-DOF_4_P_IN1_DQ', 'H1:IMC-DOF_4_Y_IN1_DQ'};
% aux_channels = {'H1:ASC-AS_A_DC_PIT_OUT_DQ', 'H1:ASC-AS_A_DC_YAW_OUT_DQ', ...
%     'H1:ASC-AS_B_DC_PIT_OUT_DQ', 'H1:ASC-AS_B_DC_YAW_OUT_DQ'};
aux_channels = {'H1:ASC-AS_A_RF45_I_PIT_OUT_DQ', 'H1:ASC-AS_A_RF45_I_YAW_OUT_DQ', ...
                'H1:ASC-AS_A_RF36_I_PIT_OUT_DQ', 'H1:ASC-AS_A_RF36_I_YAW_OUT_DQ', ...
                'H1:ASC-AS_B_RF45_I_PIT_OUT_DQ', 'H1:ASC-AS_B_RF45_I_YAW_OUT_DQ', ...
                'H1:ASC-AS_B_RF36_I_PIT_OUT_DQ', 'H1:ASC-AS_B_RF36_I_YAW_OUT_DQ'};
% aux_channels = {'H1:ASC-REFL_A_RF45_I_PIT_OUT_DQ', 'H1:ASC-REFL_A_RF45_I_YAW_OUT_DQ', ...
%                 'H1:ASC-REFL_A_RF9_I_PIT_OUT_DQ', 'H1:ASC-REFL_A_RF9_I_YAW_OUT_DQ', ...
%                 'H1:ASC-REFL_B_RF45_I_PIT_OUT_DQ', 'H1:ASC-REFL_B_RF45_I_YAW_OUT_DQ', ...
%                 'H1:ASC-REFL_B_RF9_I_PIT_OUT_DQ', 'H1:ASC-REFL_B_RF9_I_YAW_OUT_DQ'};

[aux, t, info] = g.fetch(start_time, end_time, aux_channels);
fs_aux = 1./t(2);

% band pass and compute BLRMS
band = [200 300];    % band of interest
f_lowpass = 8;      % low pass frequency
f_out = 16;         % output frequency

[B,A] = butter(4, band/(fs/2), 'bandpass');
blrms = filtfilt(B,A,darm);
blrms = blrms.^2;
[B,A] = butter(4, f_lowpass/(fs/2), 'low');
blrms = filtfilt(B,A,blrms);
blrms = blrms(1:fs/f_out:end);

% low pass the aux channels
[B,A] = butter(4, f_lowpass/(fs_aux/2), 'low');
aux = filtfilt(B,A,aux);
aux = aux(1:fs_aux/f_out:end,:);

% remove initial transient (why is it there?)
blrms = blrms(f_out:end);
aux = aux(f_out:end,:);

% create 2d histograms
lmin = min(aux);
lmax = max(aux);
npt = 100;
h = {};
blrms_range = linspace(0, 1e-23, npt);

nxplot = 4;
nyplot = 2;

figure()
for i=1:size(aux,2)
    h{i} = hist2(aux(:,i), blrms, linspace(lmin(i), lmax(i), npt), blrms_range);
    ax(i) = subplot(nxplot, nyplot, i);
    pcolor(linspace(lmin(i), lmax(i), npt), blrms_range, h{i}');
    shading flat
    title(strrep(aux_channels{i}, '_', '\_'))
end
linkaxes(ax, 'y')


% t = (1:length(blrms))/f_out;
% ax = [];
% ax(1) = subplot(211);
% plot(t, blrms)
% ax(2) = subplot(212);
% plot(t, aux(:,2))

% compute the best reconstruction of the BLRMS using the aux channels
A = [aux, aux.^2, aux.^3, aux.^4, ones(size(aux,1),1)];
B = blrms/std(blrms);
fun = @(p) sum((A*p - B).^2);
p0 = zeros(size(A,2),1);
p = fminsearch(fun, p0, optimset('display', 'iter', 'maxfunevals', 1e6, 'maxiter', 1e6));

t = (1:length(B))/f_out;
figure()
plot(t, B, t, A*p)
xlabel('Time [s]')
ylabel('Rescaled BLRMS [a.u.]')
legend('BLRMS', 'Best fit')
title('Using ASC-REFL\_RF signals')

