%% set the GPS time for a period with jitter feed forward OFF
gps = 1160186417;   % start time
dt = 600;           % duration

%% get data 
conn = nds2.connection('h1nds1', 8088);
d = conn.fetch(gps, gps+dt, {'H1:LSC-DARM_IN1_DQ', 'H1:PSL-DBB_QPD_1QY_OUT_DQ'});
darm = d(1).getData();
q1y  = d(2).getData();

%% compute coherence and transfer function
fs = d(1).getLength()/dt;
np = fs;    % change here the frequency resolution, typically 1 Hz is fine
[c,fr] = mscohere(q1y, darm, hanning(np), np/2, np,fs);
[t,fr] = tfestimate(q1y, darm, hanning(np), np/2, np,fs);

%% plot coherence
figure(1)
semilogx(fr, c)
xlabel('Frequency [Hz]', 'FontSize', 18)
ylabel('Coherence', 'FontSize', 18)
set(gca, 'FontSize', 18)
grid on

%% fit transfer function
% Set the coherence threshold and the frequency range. You might have to
% tweak those to get a reliable fit
threshold = 0.05;
fmin = 80;
fmax = 3000;

% Set the order of the transfer function. 2nd or 4th order used to work 
% just fine
order = 4;

% fit in the given range
idx = fr>fmin & fr<fmax & c>threshold;
a = fit_zpk(fr(idx), t(idx), order, 1./abs(t(idx)));

% plot the result
figure(2)
ax(1) = subplot(211);
fresponse = @(a,fr) reshape(freqresp(a, 2*pi*fr), size(fr));
loglog(fr, abs(t), fr, abs(fresponse(a, fr)))
legend('Data', 'Fit')
xlabel('Frequency [Hz]', 'FontSize', 18)
ylabel('Abs.', 'FontSize', 18)
set(gca, 'FontSize', 18)
grid on
ax(2) = subplot(212);
semilogx(fr, angle(t), fr, angle(fresponse(a, fr)))
xlabel('Frequency [Hz]', 'FontSize', 18)
ylabel('Phase [rad.]', 'FontSize', 18)
set(gca, 'FontSize', 18)
grid on
linkaxes(ax, 'x')

%% print out the result. You might want to check for high Q poles and 
% zeroes that are sometimes added at the edges of the fit band. You can
% normally remove them without problems.
fprintf(1, '\n\nZeros [Hz]:\n');
disp(a.z{1}/2/pi)
fprintf(1, 'Poles [Hz]:\n');
disp(a.p{1}/2/pi)
