set(0,'DefaultAxesFontSize',16)
startTime = 1257809856+100;
duration = 1000;
chanList1 = {'H1:CAL-PCALX_TX_PD_OUT_DQ','H1:CAL-PCALX_RX_PD_OUT_DQ','H1:CAL-PCALY_TX_PD_OUT_DQ','H1:CAL-PCALY_RX_PD_OUT_DQ','H1:CAL-DELTAL_EXTERNAL_DQ'};

segmentDuration = 10;
sampleRate = 16384;
windowFunction = hann(segmentDuration*sampleRate);
overlap = 0.5;
f0 = 395.1; % Hz

datastruct1 = get_data2(chanList1,'raw', startTime, duration);

[pxRxAvg,freq] = pwelch(double(datastruct1(2).data), windowFunction, segmentDuration*sampleRate*overlap, segmentDuration*sampleRate, sampleRate);
[pyRxAvg,freq] = pwelch(double(datastruct1(4).data), windowFunction, segmentDuration*sampleRate*overlap, segmentDuration*sampleRate, sampleRate);
[pDarmAvg,freq] = pwelch(double(datastruct1(5).data), windowFunction, segmentDuration*sampleRate*overlap, segmentDuration*sampleRate, sampleRate);
[xRx_yRxAvg,freq] = tfestimate(double(datastruct1(2).data), double(datastruct1(4).data), windowFunction, segmentDuration*sampleRate*overlap, segmentDuration*sampleRate, sampleRate);
[xRxyRxDiff,freq] = pwelch(double(datastruct1(2).data)-double(datastruct1(4).data), windowFunction, segmentDuration*sampleRate*overlap, segmentDuration*sampleRate, sampleRate);
[xRxyRxDiffCorrected,freq] = pwelch(1.0031*double(datastruct1(2).data)-double(datastruct1(4).data), windowFunction, segmentDuration*sampleRate*overlap, segmentDuration*sampleRate, sampleRate);
[xRxyRxCoh,freq] = mscohere(1.0031*double(datastruct1(2).data)-double(datastruct1(4).data), double(datastruct1(5).data), windowFunction, segmentDuration*sampleRate*overlap, segmentDuration*sampleRate, sampleRate);

% calibration into meter
sys = zpk(-2*pi*[], -2*pi*[1 1], prod(-2*pi*[1 1]));
deWhitenPCAL = squeeze(freqresp(sys, 2*pi*freq));

sys2 = zpk(-2*pi*[30 30 30 30 30 30], -2*pi*[0.3 0.3 0.3 0.3 0.3 0.3], 1e-12);
deWhitenCAL = squeeze(freqresp(sys2, 2*pi*freq));

cutFreq = freq==f0;

figure
semilogy(freq, sqrt(pDarmAvg).*abs(deWhitenCAL),'LineWidth',2)
hold on
semilogy(freq, sqrt(xRxyRxDiff).*abs(deWhitenPCAL),'LineWidth',2)
semilogy(freq, sqrt(xRxyRxDiffCorrected).*abs(deWhitenPCAL),'LineWidth',2)
grid on
xlim([394 396])
set(gca,'YScale','linear')
xlabel('Freq (Hz)')
ylabel('ASD m/ $\sqrt{Hz}$','interpreter','latex')
legend('darm','PCal (X - Y)','PCal (1.0031 * X - Y)')
