Commit 33d4a047 authored by Philip Roy Charlton's avatar Philip Roy Charlton

Only calculate PSD for the middle segment when needed for naive sigmas

parent f75d89ba
......@@ -509,7 +509,8 @@ for I = 1:params.numIntervalsTotal
for K=1:params.numTrials
if ~params.intermediate %cet------------------------------------------------
% clear arrays for PSD's
% Clear arrays for periodograms. These arrays are used to accumulate
% periodograms as we progress through the segments.
avg_data1evn = [];
avg_data1odd = [];
avg_data2evn = [];
......@@ -651,45 +652,41 @@ for I = 1:params.numIntervalsTotal
o2.deltaT, ...
o2.fbase, o2.phase);
% estimate power spectra for optimal filter - needed for naive sigmas
calPSD1 = estimatePowerSpectrum(r1(J), params.psd1.FFTLength, ...
params.psd1.Window, ...
params.psd1.OverlapLength, ...
params.psdMethod, ...
response1(J));
calPSD2 = estimatePowerSpectrum(r2(J), params.psd2.FFTLength, ...
params.psd2.Window, ...
params.psd2.OverlapLength, ...
params.psdMethod, ...
response2(J));
% calculate avg power spectra, ignoring middle segment if desired
params.midSegment = (params.numSegmentsPerInterval+1)/2;
if ( (params.ignoreMidSegment) & (J==params.midSegment) )
% do nothing
%fprintf('Ignoring middle segment\n');
% do nothing if this is the mid segment and we are not including it in the PSDs
else
%% Only accumulate periodograms on segments either side of the middle segment
%% Channel 1
% add periodograms for the current segment
[ pevn, podd ] = periodograms(r1(J).data, params.psd1.Window, params.psd1.OverlapLength, params.psd1.FFTLength, 1/r1(J).deltaT);
avg_data1evn = [ avg_data1evn, pevn ];
avg_data1odd = [ avg_data1odd, podd ];
%% Channel 2
[ pevn, podd ] = periodograms(r2(J).data, params.psd2.Window, params.psd2.OverlapLength, params.psd2.FFTLength, 1/r2(J).deltaT);
avg_data2evn = [ avg_data2evn, pevn ];
avg_data2odd = [ avg_data2odd, podd ];
end
if ( (params.writeNaiveSigmasToFiles) & (J==params.midSegment) )
% This calculates the "naive" theorerical variance, i.e.,
% estimate power spectra for optimal filter - needed for naive sigmas
calPSD1 = estimatePowerSpectrum(r1(J), params.psd1.FFTLength, ...
params.psd1.Window, ...
params.psd1.OverlapLength, ...
params.psdMethod, ...
response1(J));
calPSD2 = estimatePowerSpectrum(r2(J), params.psd2.FFTLength, ...
params.psd2.Window, ...
params.psd2.OverlapLength, ...
params.psdMethod, ...
response2(J));
% This calculates the "naive" theoretical variance ie.
% that calculated from the current segment without averaging
% over the whole interval.
% This is useful for the stationarity veto which excludes
% segments for which the naive sigma differs too much from
% the one calculated with the sliding PSD average.
% the one calculated with the sliding PSD.
[naiQ, result.naiVar, result.naiSensInt] ...
= calOptimalFilter(params.segmentDuration, params.gamma, ...
params.fRef, params.alphaExp, ...
......@@ -749,21 +746,18 @@ for I = 1:params.numIntervalsTotal
end % loop over segments J
%% This is the point to
%% - take mean or medians of periodograms
%% - average even & odd PSDs
%% - coarse-grain the PSDs
%% - calibrate PSDs with response functions
%% giving the representative PSD relevant for the middle segment
%% Channel 1
% This is the point to
% - take mean or median of periodograms
% - average even & odd PSDs
% - coarse-grain PSDs
% - calibrate PSDs with response functions
% giving the representative PSD relevant for the middle segment
result.calPSD1_avg = estimatePowerSpectrumFromPeriodograms(avg_data1evn, avg_data1odd, ...
params.psd1.FFTLength, ...
o1.fbase, 1/o1.deltaT, ...
params.psdMethod, ...
response1(J));
%% Channel 2
result.calPSD2_avg = estimatePowerSpectrumFromPeriodograms(avg_data2evn, avg_data2odd, ...
params.psd2.FFTLength, ...
o2.fbase, 1/o2.deltaT, ...
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment