% load all data, using GWData
gps = 1111190416+3600*9.5;
dt = 3600;
g = GWData;
% I have to split in four chunks, other MATLAb goes out of memory!
[data1,t,info] = g.fetch(gps, gps+dt, 'H1:CAL-DELTAL_EXTERNAL_DQ');
[data2,t,info] = g.fetch(gps+3600*1, gps+3600*1+dt, 'H1:CAL-DELTAL_EXTERNAL_DQ');
[data3,t,info] = g.fetch(gps+3600*2, gps+3600*2+dt, 'H1:CAL-DELTAL_EXTERNAL_DQ');
[data4,t,info] = g.fetch(gps+3600*3, gps+3600*3+dt, 'H1:CAL-DELTAL_EXTERNAL_DQ');
data = [data1; data2; data3; data4];
fs = 16384;
dt = 3600*4;
% decimate to 4kHz
data0 = decimate(data, 16384 / 4096);
fs = 4096;
% nothc out some known lines at low frequency
fnotch = [60, 120, 180, 56.75, 113.75, 170.5];
qpole = [100, 100, 100, 100, 100, 100];
qzero = [10000, 10000, 10000, 10000, 10000, 10000];
for i=1:numel(fnotch)
    num = [1, 2*pi*fnotch(i)/qzero(i), (2*pi*fnotch(i))^2];
    den = [1, 2*pi*fnotch(i)/qpole(i), (2*pi*fnotch(i))^2];
    [B,A] = bilinear(num, den, fs, fnotch(i));
    data0 = filtfilt(B,A, data0);
end
% compute spectrogram
np = fs*32;
[S,F,T] = spectrogram(data0, hanning(np), np/2, np, fs);
% track lines
fr166 = detectLine(F,S, [160,170]);
fr230 = detectLine(F,S, [230,240]);
fr250 = detectLine(F,S, [245,260]);
fr415 = detectLine(F,S, [410,430]);
fr970 = detectLine(F,S, [960,980]);
% save results
L(:,1) = T;
L(:,2) = fr166;
L(:,3) = fr230;
L(:,4) = fr250;
L(:,5) = fr415;
L(:,6) = fr970;
save -ascii wandering_lines.txt L
