Commit efdadc90 authored by Philip Roy Charlton's avatar Philip Roy Charlton

Merge branch 'O3Development' into 'master'

O3 development

See merge request stochastic/stochastic!1
parents 1b892333 24edace2
......@@ -120,7 +120,7 @@ function params = setDefaultParams(paramsFile, startTime, jobDuration)
params.signalType = 'const';
% Default index of power law used for signal type const
parames.powerIndex = 0;
params.powerIndex = 0;
% value of Omega_gw(f_Ref) for simulated SB signal
params.simOmegaRef1 = 0;
......
......@@ -677,7 +677,7 @@ for I = 1:params.numIntervalsTotal
% 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.
[naiQ, result.naiVar, naiSensInt] ...
[naiQ, result.naiVar, result.naiSensInt] ...
= calOptimalFilter(params.segmentDuration, params.gamma, ...
params.fRef, params.alphaExp, ...
calPSD1, calPSD2, ...
......
......@@ -85,6 +85,7 @@ function writeOutputToMatFile(params, I, K, segmentStartTime, result)
if (params.doAllSkyComparison & params.writeSensIntsToFiles)
% write sensitivity integrands to output files
params.matFile.sensIntAllSky(I, K) = result.sensIntAllSky;
params.matFile.naiSensInt(I, K) = result.naiSensInt;
end
% Q
......
%
% coarseGrain --- coarse grain a frequency-series
%
% coarseGrain(x, flowy, deltaFy, Ny) returns a frequency-series structure
% coarse-grained to the frequency values f = flowy + deltaFy*[0:Ny-1].
%
% coarseGrain also returns the frequency indices of the lowest and highest
% frequency bins of x that overlap with y (0 <= index1 <= length(x);
% 1 <= index2 <= length(x)+1) and the fractional contributions from these
% frequency bins (note that these indices start from 0 and and are 1 less
% than the corresponding Matlab indices). The fractional contribution is
% the fraction of the fine bin that overlaps with any part of the coarse
% series.
%
% The method used is to treat the x-values
%
% x(k) = x.flow + (k - 1)*x.deltaF
%
% as bin centres and the value x.data(k) as the average value over the
% the bin ie.
%
% (1/x.deltaF)*(integral from the lower edge of the bin to upper edge)
%
% The coarse graining can then be performed by finding the integral
% from the start of the fine series just below the coarse series (via
% the cumulative sum), interpolating to the coarse series, and taking
% differences to recover the average value for each coarse bin.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y, index1, index2, frac1, frac2]=coarseGrain_new(x, flowy, deltaFy, Ny)
% Set the metadata for the coarse-grained sequence
y.flow = flowy;
y.deltaF = deltaFy;
if (isfield(x, 'symmetry'))
y.symmmetry = x.symmetry;
else
y.symmetry = 0;
end;
% Length of coarse series
Nx = length(x.data);
% Lower edge of the first bin of the fine series
xLowerEdge = x.flow - 0.5*x.deltaF;
% Upper edge of the last bin of the fine series
xUpperEdge = x.flow + (Nx - 0.5)*x.deltaF;
% yi(k) is the lower edge of bin k for the coarse-grained series
yi = y.flow + y.deltaF*[-0.5:1:Ny-0.5].';
% Lower edge of the first bin of the coarse series
yLowerEdge = yi(1);
% Upper edge of the last bin of the coarse series
yUpperEdge = yi(end);
% Error checking
if (Nx <= 0)
error('Length of fine series is invalid');
end;
if (Ny <= 0)
error('Length of coarse series is invalid');
end;
if (x.deltaF <= 0)
error('Spacing of fine series is invalid');
end;
if (y.deltaF <= 0)
error('Spacing of coarse series is invalid');
end;
if (y.deltaF < x.deltaF)
error('Frequency spacing of coarse series is smaller than fine series');
end;
if (yLowerEdge < xLowerEdge)
error('Start of coarse series is less than start of fine series');
end;
if (yUpperEdge > xUpperEdge)
error('End of coarse series is more than end of fine series');
end;
% xlow is the index of the last bin whose lower edge is <= the
% lower edge of the coarse-grained sequence, that is
% x(xlow) <= y(1) < x(xlow+1)
xlow = 1 + floor((yLowerEdge - xLowerEdge)/x.deltaF);
% xhi is the index of the last bin whose upper edge is >= the
% lower edge of the coarse-grained sequence, that is
% x(xhi-1) < y(end) <= x(xhi)
xhi = ceil((yUpperEdge - xLowerEdge)/x.deltaF);
% xi is the array of frequencies of the lower edge of each bin, that is,
% xi(k) is the lower edge of bin k, which is
% x.flow + (k-1)*x.deltaF - 0.5*x.deltaF = x.flow + (k-1.5)*x.deltaF
% This is only calculated for the bins that the coarse series overlaps with,
% that is, bins [xlow:xhi]
xi = x.flow + [(xlow-1.5):1:(xhi-0.5)].'*x.deltaF;
% Integrate the original function so that the value fi(k) is
% the integral from the lower edge of the fine series xi(1)
% to xi(k). The 0 is inserted so that the fi(1) is 0, as it should be,
% so that we can interpolate the integral to the coarse series
fi = [0; x.deltaF*cumsum(x.data(xlow:xhi)) ];
% Interpolate the integrated function using the lower edges of each bin in
% the coarse series as the ordinates
y.data = interp1(xi, fi, yi, 'linear');
% Take the difference to get the integral over each bin of the coarse series,
% that is, y.data(k) becomes the integral from y(k) to y(k+1)
y.data = (1/y.deltaF)*diff(y.data);
% For compatability with the original coarseGrain, we also return the
% frequency index of the first and last bins of the fine series that
% overlap with the coarse series (these are the array index - 1)
index1 = xlow - 1;
index2 = xhi - 1;
% Fraction of the end bins of the fine series that overlap with the
% coarse series
frac1 = (xi(2) - yLowerEdge)/x.deltaF;
frac2 = (yUpperEdge - xi(end-1))/x.deltaF;
return;
......@@ -14,7 +14,7 @@ segment=15, trial =1, start sec = 1126627729, CC stat = -1.111019e-04, theor sig
segment=16, trial =1, start sec = 1126627759, CC stat = 1.848009e-04, theor sigma = 5.036686e-04
segment=17, trial =1, start sec = 1126627789, CC stat = -6.454592e-04, theor sigma = 4.764881e-04
segment=18, trial =1, start sec = 1126627819, CC stat = 1.081851e-03, theor sigma = 4.761411e-04
segment=19, trial =1, start sec = 1126627849, CC stat = -4.415706e-04, theor sigma = 4.827310e-04
segment=19, trial =1, start sec = 1126627849, CC stat = -4.415707e-04, theor sigma = 4.827310e-04
segment=20, trial =1, start sec = 1126627879, CC stat = -5.005905e-04, theor sigma = 4.671271e-04
segment=21, trial =1, start sec = 1126627909, CC stat = 4.480586e-04, theor sigma = 4.535013e-04
segment=22, trial =1, start sec = 1126627939, CC stat = 2.599749e-04, theor sigma = 4.543806e-04
......
......@@ -14,7 +14,7 @@ segment=15, trial =1, start sec = 1126627729, CC stat = -1.111019e-04, theor sig
segment=16, trial =1, start sec = 1126627759, CC stat = 1.848009e-04, theor sigma = 5.036686e-04
segment=17, trial =1, start sec = 1126627789, CC stat = -6.454592e-04, theor sigma = 4.764881e-04
segment=18, trial =1, start sec = 1126627819, CC stat = 1.081851e-03, theor sigma = 4.761411e-04
segment=19, trial =1, start sec = 1126627849, CC stat = -4.415706e-04, theor sigma = 4.827310e-04
segment=19, trial =1, start sec = 1126627849, CC stat = -4.415707e-04, theor sigma = 4.827310e-04
segment=20, trial =1, start sec = 1126627879, CC stat = -5.005905e-04, theor sigma = 4.671271e-04
segment=21, trial =1, start sec = 1126627909, CC stat = 4.480586e-04, theor sigma = 4.535013e-04
segment=22, trial =1, start sec = 1126627939, CC stat = 2.599749e-04, theor sigma = 4.543806e-04
......
......@@ -14,7 +14,7 @@ segment=15, trial =1, start sec = 1126627729, CC stat = -1.111019e-04, theor sig
segment=16, trial =1, start sec = 1126627759, CC stat = 1.848009e-04, theor sigma = 5.036686e-04
segment=17, trial =1, start sec = 1126627789, CC stat = -6.454592e-04, theor sigma = 4.764881e-04
segment=18, trial =1, start sec = 1126627819, CC stat = 1.081851e-03, theor sigma = 4.761411e-04
segment=19, trial =1, start sec = 1126627849, CC stat = -4.415706e-04, theor sigma = 4.827310e-04
segment=19, trial =1, start sec = 1126627849, CC stat = -4.415707e-04, theor sigma = 4.827310e-04
segment=20, trial =1, start sec = 1126627879, CC stat = -5.005905e-04, theor sigma = 4.671271e-04
segment=21, trial =1, start sec = 1126627909, CC stat = 4.480586e-04, theor sigma = 4.535013e-04
segment=22, trial =1, start sec = 1126627939, CC stat = 2.599749e-04, theor sigma = 4.543806e-04
......
......@@ -6,7 +6,7 @@ segment=7, trial =1, start sec = 1126626195, CC stat = 2.825762e-04, theor sigma
segment=8, trial =1, start sec = 1126626225, CC stat = -3.154353e-04, theor sigma = 4.581993e-04
segment=9, trial =1, start sec = 1126626255, CC stat = 2.459836e-04, theor sigma = 4.594404e-04
segment=10, trial =1, start sec = 1126626285, CC stat = 7.509449e-04, theor sigma = 4.547512e-04
segment=11, trial =1, start sec = 1126626315, CC stat = 9.446346e-06, theor sigma = 4.539911e-04
segment=11, trial =1, start sec = 1126626315, CC stat = 9.446347e-06, theor sigma = 4.539911e-04
segment=12, trial =1, start sec = 1126626345, CC stat = -4.455225e-04, theor sigma = 4.421551e-04
segment=13, trial =1, start sec = 1126626375, CC stat = 3.781912e-04, theor sigma = 4.383439e-04
segment=14, trial =1, start sec = 1126626405, CC stat = 4.334577e-04, theor sigma = 4.400128e-04
......@@ -19,18 +19,18 @@ segment=20, trial =1, start sec = 1126626585, CC stat = -6.365942e-05, theor sig
segment=21, trial =1, start sec = 1126626615, CC stat = 4.523975e-04, theor sigma = 4.494876e-04
segment=22, trial =1, start sec = 1126626645, CC stat = 1.622646e-04, theor sigma = 4.452459e-04
segment=23, trial =1, start sec = 1126626675, CC stat = -2.868838e-04, theor sigma = 4.333874e-04
segment=24, trial =1, start sec = 1126626705, CC stat = 3.698967e-04, theor sigma = 4.352248e-04
segment=24, trial =1, start sec = 1126626705, CC stat = 3.698968e-04, theor sigma = 4.352248e-04
segment=25, trial =1, start sec = 1126626735, CC stat = 5.270730e-04, theor sigma = 4.448423e-04
segment=26, trial =1, start sec = 1126626765, CC stat = 5.086034e-04, theor sigma = 4.650624e-04
segment=27, trial =1, start sec = 1126626795, CC stat = -1.988094e-04, theor sigma = 4.585238e-04
segment=28, trial =1, start sec = 1126626825, CC stat = 2.549785e-04, theor sigma = 4.492477e-04
segment=29, trial =1, start sec = 1126626855, CC stat = 5.085108e-04, theor sigma = 4.494080e-04
segment=30, trial =1, start sec = 1126626885, CC stat = 3.134164e-04, theor sigma = 4.586114e-04
segment=31, trial =1, start sec = 1126626915, CC stat = -7.783852e-06, theor sigma = 4.631874e-04
segment=31, trial =1, start sec = 1126626915, CC stat = -7.783851e-06, theor sigma = 4.631874e-04
segment=32, trial =1, start sec = 1126626945, CC stat = -3.925304e-04, theor sigma = 4.598011e-04
segment=33, trial =1, start sec = 1126626975, CC stat = -4.352466e-04, theor sigma = 4.790235e-04
segment=34, trial =1, start sec = 1126627005, CC stat = 3.984703e-04, theor sigma = 4.748823e-04
segment=35, trial =1, start sec = 1126627035, CC stat = 9.681774e-05, theor sigma = 4.552185e-04
segment=35, trial =1, start sec = 1126627035, CC stat = 9.681773e-05, theor sigma = 4.552185e-04
segment=36, trial =1, start sec = 1126627065, CC stat = 4.286636e-04, theor sigma = 4.500477e-04
segment=37, trial =1, start sec = 1126627095, CC stat = 1.291056e-04, theor sigma = 4.710300e-04
segment=38, trial =1, start sec = 1126627125, CC stat = 6.644900e-04, theor sigma = 4.750279e-04
......@@ -51,13 +51,13 @@ segment=15, trial =1, start sec = 1126627729, CC stat = -1.111019e-04, theor sig
segment=16, trial =1, start sec = 1126627759, CC stat = 1.848009e-04, theor sigma = 5.036686e-04
segment=17, trial =1, start sec = 1126627789, CC stat = -6.454592e-04, theor sigma = 4.764881e-04
segment=18, trial =1, start sec = 1126627819, CC stat = 1.081851e-03, theor sigma = 4.761411e-04
segment=19, trial =1, start sec = 1126627849, CC stat = -4.415706e-04, theor sigma = 4.827310e-04
segment=19, trial =1, start sec = 1126627849, CC stat = -4.415707e-04, theor sigma = 4.827310e-04
segment=20, trial =1, start sec = 1126627879, CC stat = -5.005905e-04, theor sigma = 4.671271e-04
segment=21, trial =1, start sec = 1126627909, CC stat = 4.480586e-04, theor sigma = 4.535013e-04
segment=22, trial =1, start sec = 1126627939, CC stat = 2.599749e-04, theor sigma = 4.543806e-04
segment=23, trial =1, start sec = 1126627969, CC stat = -4.062091e-04, theor sigma = 4.388201e-04
segment=3, trial =1, start sec = 1126643603, CC stat = 5.597523e-04, theor sigma = 4.145244e-04
segment=4, trial =1, start sec = 1126643633, CC stat = 4.548432e-05, theor sigma = 4.289550e-04
segment=4, trial =1, start sec = 1126643633, CC stat = 4.548431e-05, theor sigma = 4.289550e-04
segment=5, trial =1, start sec = 1126643663, CC stat = -2.131570e-04, theor sigma = 4.390990e-04
segment=6, trial =1, start sec = 1126643693, CC stat = -9.293389e-04, theor sigma = 4.239447e-04
segment=7, trial =1, start sec = 1126643723, CC stat = 1.020191e-03, theor sigma = 4.138207e-04
......@@ -66,7 +66,7 @@ segment=9, trial =1, start sec = 1126643783, CC stat = -7.305027e-05, theor sigm
segment=10, trial =1, start sec = 1126643813, CC stat = 4.326493e-04, theor sigma = 4.137791e-04
segment=11, trial =1, start sec = 1126643843, CC stat = -1.080139e-04, theor sigma = 4.066223e-04
segment=12, trial =1, start sec = 1126643873, CC stat = -5.622439e-04, theor sigma = 4.073220e-04
segment=13, trial =1, start sec = 1126643903, CC stat = -6.418782e-05, theor sigma = 4.183537e-04
segment=13, trial =1, start sec = 1126643903, CC stat = -6.418781e-05, theor sigma = 4.183537e-04
segment=14, trial =1, start sec = 1126643933, CC stat = -2.308930e-04, theor sigma = 4.203017e-04
segment=15, trial =1, start sec = 1126643963, CC stat = 6.916928e-04, theor sigma = 4.116834e-04
segment=16, trial =1, start sec = 1126643993, CC stat = -2.682871e-04, theor sigma = 4.062399e-04
......@@ -94,8 +94,8 @@ segment=37, trial =1, start sec = 1126644623, CC stat = -4.210991e-04, theor sig
segment=38, trial =1, start sec = 1126644653, CC stat = -1.480465e-05, theor sigma = 4.061845e-04
segment=39, trial =1, start sec = 1126644683, CC stat = 1.791653e-04, theor sigma = 4.177789e-04
segment=40, trial =1, start sec = 1126644713, CC stat = -5.655936e-05, theor sigma = 4.181445e-04
segment=41, trial =1, start sec = 1126644743, CC stat = 6.392874e-05, theor sigma = 4.092614e-04
segment=42, trial =1, start sec = 1126644773, CC stat = 1.122695e-05, theor sigma = 4.149807e-04
segment=41, trial =1, start sec = 1126644743, CC stat = 6.392873e-05, theor sigma = 4.092614e-04
segment=42, trial =1, start sec = 1126644773, CC stat = 1.122694e-05, theor sigma = 4.149807e-04
segment=43, trial =1, start sec = 1126644803, CC stat = -1.060125e-05, theor sigma = 4.130287e-04
segment=44, trial =1, start sec = 1126644833, CC stat = 2.608625e-04, theor sigma = 4.083948e-04
segment=45, trial =1, start sec = 1126644863, CC stat = 2.561162e-04, theor sigma = 4.105651e-04
......
......@@ -16,7 +16,7 @@ writeSensIntsToFiles true
writeCoherenceToFiles false
writeCalPSD1sToFiles false
writeCalPSD2sToFiles false
writeResultsToScreen false
writeResultsToScreen true
% high-pass parameters
useCascadeFilter1 true
......
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
......@@ -12,7 +12,13 @@ duration = 60; % sec
tlow = 0; % gps start time
fsample = 4096; % Hz
gammaLM_coeffsPath='/Users/joeromano/src/matapps/src/searches/stochastic/CrossCorr/';
% Needs to point to the CrossCorr/mats directory
gammaLM_coeffsPath='../mats';
% Simfiles path. These large files are part of the sgwb repository
% Change to point to your local installation
simfilesPath = '/home/charlton/Devel/sgwb/trunk/S5/input/simfiles';
Lmax = 5;
isSpH = true;
......@@ -26,15 +32,15 @@ g2lm = calGammaLM(gammaLM_coeffsPath, 'LL', Lmax, numFreqs, flow, deltaF);
glm = calGammaLM(gammaLM_coeffsPath, 'HL', Lmax, numFreqs, flow, deltaF);
% load signal power spectrum H(f) from file
Hf = load('/Users/joeromano/src/sgwb/S5/input/simfiles/HSource_1e_0.txt');
Hf = load([ simfilesPath '/' 'HSource_1e_0.txt']);
intLog = 1;
% read in a skymap
filename = '/Users/joeromano/src/sgwb/S5/input/simfiles/monopole_plmComplex_L5.txt'
filename = '/Users/joeromano/src/sgwb/S5/input/simfiles/dipole_plmComplex_L5.txt'
filename = '/Users/joeromano/src/sgwb/S5/input/simfiles/quadrupole_plmComplex_L5.txt'
filename = '/Users/joeromano/src/sgwb/S5/input/simfiles/octupole_plmComplex_L5.txt'
filename = '/Users/joeromano/src/sgwb/S5/input/simfiles/pointSource_ra6_dec45_plmComplex_L5.txt'
filename = [ simfilesPath '/' 'monopole_plmComplex_L5.txt' ];
filename = [ simfilesPath '/' 'dipole_plmComplex_L5.txt' ];
filename = [ simfilesPath '/' 'quadrupole_plmComplex_L5.txt' ];
filename = [ simfilesPath '/' 'octupole_plmComplex_L5.txt' ];
filename = [ simfilesPath '/' 'pointSource_ra6_dec45_plmComplex_L5.txt' ];
fileType = 2;
numBins = 1;
InjectAsSpH = true;
......
......@@ -9,6 +9,11 @@ function test_coarseGrain()
reps = 10;
pass = true;
if (exist('test_coarseGrain.out'))
delete('test_coarseGrain.out');
end;
diary('test_coarseGrain.out');
% The first set of tests checks that the values returned by the new function have
% the expected properties.
%
......@@ -266,6 +271,8 @@ function test_coarseGrain()
fprintf('FAIL\n');
end;
diary off;
return;
function pass=test_general(f0_x, df_x, Nx, f0_y, df_y, Ny)
......@@ -287,7 +294,7 @@ function pass=test_coarse(f0_x, df_x, Nx, f0_y, df_y, Ny)
x = generate_data(f0_x, df_x, Nx);
[y, index1, index2, frac1, frac2] = coarseGrain_new(x, f0_y, df_y, Ny);
[y, index1, index2, frac1, frac2] = coarseGrain(x, f0_y, df_y, Ny);
% Check the characteristics of the coarse-grained series
......@@ -359,7 +366,7 @@ function pass=test_v3(f0_x, df_x, Nx, f0_y, df_y, Ny)
x = generate_data(f0_x, df_x, Nx);
[y, index1, index2, frac1, frac2] = coarseGrain_new(x, f0_y, df_y, Ny);
[y, index1, index2, frac1, frac2] = coarseGrain(x, f0_y, df_y, Ny);
[y1, i1, i2, f1, f2] = coarseGrain_v3(x, f0_y, df_y, Ny);
if (length(y.data) == length(y1.data))
......@@ -432,7 +439,7 @@ function pass=test_time(f0_x, df_x, Nx, f0_y, df_y, Ny, reps)
tic;
for k = 1:reps
[y, index1, index2, frac1, frac2] = coarseGrain_new(x, f0_y, df_y, Ny);
[y, index1, index2, frac1, frac2] = coarseGrain(x, f0_y, df_y, Ny);
end;
t = toc;
fprintf('New time = %f secs/use\n', t/reps);
......
** Testing f0_x = 1.00, df_x = 1.00, Nx = 1, f0_y = 1.00, df_y = 1.00, Ny = 1
PASS coarse length, expected 1, received 1
PASS lower edge 0.500000 <= 0.500000 < 1.500000
PASS upper edge 0.500000 < 1.500000 <= 1.500000
PASS lower fraction, received 1.000000, expected 1.000000
PASS upper fraction, received 1.000000, expected 1.000000
** Testing f0_x = 1.00, df_x = 1.00, Nx = 2, f0_y = 1.00, df_y = 1.00, Ny = 2
PASS coarse length, expected 2, received 2
PASS lower edge 0.500000 <= 0.500000 < 1.500000
PASS upper edge 1.500000 < 2.500000 <= 2.500000
PASS lower fraction, received 1.000000, expected 1.000000
PASS upper fraction, received 1.000000, expected 1.000000
** Testing f0_x = 1.00, df_x = 1.00, Nx = 8, f0_y = 1.00, df_y = 1.00, Ny = 8
PASS coarse length, expected 8, received 8
PASS lower edge 0.500000 <= 0.500000 < 1.500000
PASS upper edge 7.500000 < 8.500000 <= 8.500000
PASS lower fraction, received 1.000000, expected 1.000000
PASS upper fraction, received 1.000000, expected 1.000000
** Testing f0_x = 1.00, df_x = 1.00, Nx = 2, f0_y = 1.00, df_y = 1.00, Ny = 1
PASS coarse length, expected 1, received 1
PASS lower edge 0.500000 <= 0.500000 < 1.500000
PASS upper edge 0.500000 < 1.500000 <= 1.500000
PASS lower fraction, received 1.000000, expected 1.000000
PASS upper fraction, received 1.000000, expected 1.000000
** Testing f0_x = 1.00, df_x = 1.00, Nx = 2, f0_y = 2.00, df_y = 1.00, Ny = 1
PASS coarse length, expected 1, received 1
PASS lower edge 1.500000 <= 1.500000 < 2.500000
PASS upper edge 1.500000 < 2.500000 <= 2.500000
PASS lower fraction, received 1.000000, expected 1.000000
PASS upper fraction, received 1.000000, expected 1.000000
** Testing f0_x = 1.00, df_x = 1.00, Nx = 8, f0_y = 2.00, df_y = 1.00, Ny = 5
PASS coarse length, expected 5, received 5
PASS lower edge 1.500000 <= 1.500000 < 2.500000
PASS upper edge 5.500000 < 6.500000 <= 6.500000
PASS lower fraction, received 1.000000, expected 1.000000
PASS upper fraction, received 1.000000, expected 1.000000
** Testing f0_x = 1.00, df_x = 1.00, Nx = 8, f0_y = 1.25, df_y = 1.00, Ny = 5
PASS coarse length, expected 5, received 5
PASS lower edge 0.500000 <= 0.750000 < 1.500000
PASS upper edge 5.500000 < 5.750000 <= 6.500000
PASS lower fraction, received 0.750000, expected 0.750000
PASS upper fraction, received 0.250000, expected 0.250000
** Testing f0_x = 0.30, df_x = 0.01, Nx = 245760, f0_y = 0.42, df_y = 0.25, Ny = 6824
PASS coarse length, expected 6824, received 6824
PASS lower edge 0.295833 <= 0.295833 < 0.304167
PASS upper edge 1706.287500 < 1706.295833 <= 1706.295833
PASS lower fraction, received 1.000000, expected 1.000000
PASS upper fraction, received 1.000000, expected 1.000000
** Testing f0_x = 0.00, df_x = 0.01, Nx = 245760, f0_y = 342.12, df_y = 0.25, Ny = 6824
PASS coarse length, expected 6824, received 6824
PASS lower edge 341.995833 <= 341.995833 < 342.004167
PASS upper edge 2047.987500 < 2047.995833 <= 2047.995833
PASS lower fraction, received 1.000000, expected 1.000000
PASS upper fraction, received 1.000000, expected 1.000000
** Testing f0_x = 0.00, df_x = 0.01, Nx = 245760, f0_y = 20.00, df_y = 0.25, Ny = 6824
PASS coarse length, expected 6824, received 6824
PASS lower edge 19.870833 <= 19.875000 < 19.879167
PASS upper edge 1725.870833 < 1725.875000 <= 1725.879167
PASS lower fraction, received 0.500000, expected 0.500000
PASS upper fraction, received 0.500000, expected 0.500000
** Testing f0_x = 0.00, df_x = 0.01, Nx = 245760, f0_y = 20.00, df_y = 0.03, Ny = 54592
PASS coarse length, expected 54592, received 54592
PASS lower edge 19.979167 <= 19.984375 < 19.987500
PASS upper edge 1725.979167 < 1725.984375 <= 1725.987500
PASS lower fraction, received 0.375000, expected 0.375000
PASS upper fraction, received 0.625000, expected 0.625000
** Testing f0_x = 1.00, df_x = 1.00, Nx = 0, f0_y = 1.00, df_y = 1.00, Ny = 0
PASS caught exception 'Length of fine series is invalid'
** Testing f0_x = 1.00, df_x = 1.00, Nx = 8, f0_y = 1.00, df_y = 0.50, Ny = 8
PASS caught exception 'Frequency spacing of coarse series is smaller than fine series'
** Testing f0_x = 1.00, df_x = 1.00, Nx = 8, f0_y = 0.50, df_y = 1.00, Ny = 8
PASS caught exception 'Start of coarse series is less than start of fine series'
** Testing f0_x = 1.00, df_x = 1.00, Nx = 8, f0_y = 1.50, df_y = 1.00, Ny = 8
PASS caught exception 'End of coarse series is more than end of fine series'
** Testing f0_x = 0.00, df_x = 1.00, Nx = 8, f0_y = 0.25, df_y = 1.00, Ny = 6
PASS coarse length, expected 6, received 6
PASS with max difference 0.0000e+00
PASS with index1 = 0 and old index1 = 0
PASS with index2 = 6 and old index2 = 6
PASS with frac1 = 0.75 and old frac1 = 0.75
PASS with frac2 = 0.25 and old frac2 = 0.25
** Testing f0_x = 0.00, df_x = 0.01, Nx = 245760, f0_y = 20.00, df_y = 0.25, Ny = 6824
PASS coarse length, expected 6824, received 6824
PASS with max difference 3.8524e-06
PASS with index1 = 2385 and old index1 = 2385
PASS with index2 = 207105 and old index2 = 207105
PASS with frac1 = 0.50 and old frac1 = 0.50
PASS with frac2 = 0.50 and old frac2 = 0.50
** Testing f0_x = 0.00, df_x = 0.01, Nx = 245760, f0_y = 20.00, df_y = 0.03, Ny = 54592
PASS coarse length, expected 54592, received 54592
PASS with max difference 2.6319e-05
PASS with index1 = 2398 and old index1 = 2398
PASS with index2 = 207118 and old index2 = 207118
PASS with frac1 = 0.38 and old frac1 = 0.38
PASS with frac2 = 0.62 and old frac2 = 0.62
** Testing time for f0_x = 0.00, df_x = 0.01, Nx = 245760, f0_y = 20.00, df_y = 0.25, Ny = 6824
V3s time = 0.017427 secs/use
New time = 0.005767 secs/use
** Testing time for f0_x = 0.00, df_x = 0.01, Nx = 245760, f0_y = 20.00, df_y = 0.03, Ny = 54592
V3s time = 0.122165 secs/use
New time = 0.006693 secs/use
PASS: all
This diff is collapsed.
[Warning: MATLAB has disabled some advanced graphics rendering features by switching to software OpenGL. For more information, click <a href="matlab:opengl('problems')">here</a>.]
Not applying delta sigma cut or user cuts.
Final point estimate = 1.057090e-06
Final error bar = 7.206430e-07
......
close all;
clear;
% Input files path. These large files are part of the sgwb repository
% Change to point to your local installation
inputFilesPath = '/home/charlton/Devel/sgwb/trunk/S4/input';
f=(00:0.25:2000)';
source='/home/sballmer/stochastic/sgwb/S4/input/radiometer/sourceDirection.txt';
HS='/home/sballmer/stochastic/sgwb/S4/input/radiometer/Hf.txt';
source=[ inputFilesPath '/' 'radiometer/sourceDirection.txt' ];
HS=[ inputFilesPath '/' 'radiometer/Hf.txt' ];
target=[0 45];
target='/home/sballmer/stochastic/sgwb/S4/input/radiometer/sourceDirection.txt';
target=[ inputFilesPath '/' 'radiometer/sourceDirection.txt' ];
H ='/home/sballmer/stochastic/sgwb/S4/input/radiometer/Hf.txt';
Hg='/home/sballmer/stochastic/sgwb/S4/input/radiometer/Hgaussian.txt';
H =[ inputFilesPath '/' 'radiometer/Hf.txt' ];
Hg=[ inputFilesPath '/' 'radiometer/Hgaussian.txt' ];
calH1='/home/sballmer/stochastic/sgwb/S4/input/calibration/S04-H1-CAL-DERRRESPONSE-793099715.mat';
calL1='/home/sballmer/stochastic/sgwb/S4/input/calibration/S04-L1-CAL-DERRRESPONSE-792579613.mat';
calH1=[ inputFilesPath '/' 'calibration/S04-H1-CAL-DERRRESPONSE-793099715.mat' ];
calL1=[ inputFilesPath '/' 'calibration/S04-L1-CAL-DERRRESPONSE-792579613.mat' ];
s=load(source);
......
function makeDeltaSigmaCutsInPostProc(jobfile, stochOutputPrefix, dsc, doOverlap, segmentDuration, deltaF, outname, alpha)
%
% This script goes through the *_naivesigmas* output files for all of the
% jobs processed by stochastic.m, calculates a user-specified delta sigma
% cut (taking into account bias factor), and writes to file
% sGc 02/15/16
%
% Input
% - jobfile = jobfile used for particular run of stochastic.m
% - stochOutputPrefix = prefix for output of stochastic.m
% - dsc = desired delta sigma cut
% - doOverlap = 0 for no stochastic.m overlap of time segments (non-standard),
% 1 for overlap
% - segmentDuration = segment duration used for stochastic.m run
% - deltaF = frequency bin size used in stochastic.m
% - outname = name of output bad GPS times based on delta sigma cut
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%USER-SPECIFIED INPUT%%%
%Input
%jobfile='/home/sgwynne.crowder/sgwb/O1_zerolag/C02/input/jobfiles/JOB-FILE-1126623617-1136649617-master_C02-events-removed.dat';
%stochOutputPrefix='/home/sgwynne.crowder/sgwb/O1_zerolag/C02/output/a3/a3';
%Desired delta sigma cut
%dsc=0.2;
%Info for calculating bias factor for naive and sliding sigmas (info about bias factor on ilog & Stefan Ballmer's thesis)
%doOverlap=1;
%segmentDuration=192;
%deltaF=1/32;
if ~exist('alpha','var')
alpha=0;
end
fref=1;
if doOverlap
segs=segmentDuration*deltaF*2-1;
else
segs=segmentDuration*deltaF;
end
nn=2*9/11*segs; %9/11 for Welch factor
bf_ss=nn/(nn-1); %sliding sigma bias factor
nn=9/11*segs;
bf_ns=nn/(nn-1); %naive sigma bias factor
%bf_ss=1; %trying getting rid of bias factor
%bf_ns=1; %trying getting rid of bias factor
%Output
%outname='/home/sgwynne.crowder/sgwb/O1_zerolag/C02/output/a0_postproc/badGPSTimes_a3.txt';
%%%%%%%%%%%%%%%%%%%%%%%%%%
%Parameters
jobs=load(jobfile);
num_jobs=size(jobs,1);
%Process
fid=fopen(outname,'w+');
for ii=1:num_jobs %step through jobs
%sigmasfilename=[stochOutputPrefix '_naivesigmas.job' num2str(ii) '.trial1.dat'];
disp(['Processing job ' num2str(ii) ' of ' num2str(num_jobs)]);
sigmasfilename=[stochOutputPrefix '.job' num2str(ii) '.mat']
if exist(sigmasfilename)
tmp=load(sigmasfilename); %column 1 = start sec, column 2 = naive sigma, column 3 = sliding sigma
if ~isfield(tmp,'segmentStartTime')
disp('No data for this job')
continue
end
naiSensInt=tmp.naiSensInt;
sensInt=tmp.sensInt;
segmentStartTime=tmp.segmentStartTime;
% define frequency array
flow=naiSensInt(1).flow;
deltaF=naiSensInt(1).deltaF;
nbins=length(naiSensInt(1).data);
fmax=flow+(nbins-1)*deltaF;
freqs=transpose(flow:deltaF:fmax);
% Define H(f)=((f/fref)^alpha
Hf=(freqs/fref).^alpha;
% Rescale the naive snsiviity integrand and integrate to get naive sigma
nsegs=length(naiSensInt)
naiSigmaAlpha=zeros(nsegs,1);
ccSigmaAlpha=zeros(nsegs,1);
for jj=1:nsegs
naiSensIntAlpha=naiSensInt(jj).data.*Hf.^2;
naiSigmaAlpha(jj)=1/sqrt(sum(naiSensIntAlpha*deltaF));
sensIntAlpha=sensInt(jj).data.*Hf.^2;
ccSigmaAlpha(jj)=1/sqrt(sum(sensIntAlpha*deltaF));
end
sigmas=[segmentStartTime,naiSigmaAlpha,ccSigmaAlpha];
[m,n]=size(sigmas);
if m~=0
for jj=1:m %write times exceeding delta sigma cut to file
if abs(sigmas(jj,3)*bf_ss-sigmas(jj,2)*bf_ns)/(sigmas(jj,3)*bf_ss)>dsc %apply bias factor
%fprintf(fid,'%i %e %e\n',sigmas(jj,1),sigmas(jj,2),sigmas(jj,3));
fprintf(fid,'%i\n',sigmas(jj,1));
end
end
end
end
clear naiSigmaAlpha;
clear ccSigmaAlpha;
clear freqs;
clear Hf;
clear tmp;
end %jobs
fclose(fid)
return;
% Parameters which are the same for all runs of makeDeltaSigmaCutsInPostProc
jobfile='../input/jobs/JOB-FILE-1238166018-1239458471.dat';
dsc=0.2;
doOverlap=1;
segmentDuration=192;
deltaF=0.03125;
% For a=-5
stochOutputPrefix_aminus5='../output/a0/a0/H1L1';
outname_aminus5='./tmp_a-5.txt';
makeDeltaSigmaCutsInPostProc_naiSensInt(jobfile,stochOutputPrefix_aminus5,dsc,doOverlap,segmentDuration,deltaF,outname_aminus5,-5);
% For a=0
stochOutputPrefix_a0='../output/a0/a0/H1L1';
outname_a0='./tmp_a0.txt';
makeDeltaSigmaCutsInPostProc_naiSensInt(jobfile,stochOutputPrefix_a0,dsc,doOverlap,segmentDuration,deltaF,outname_a0,0);
% For a=3
stochOutputPrefix_a3='../output/a0/a0/H1L1';
outname_a3='./tmp_a3.txt';
makeDeltaSigmaCutsInPostProc_naiSensInt(jobfile,stochOutputPrefix_a3,dsc,doOverlap,segmentDuration,deltaF,outname_a3,3);
% Super cut
outname_supercut='./badGPSTimes_naiSensInt.dat'
produceSuperCut(outname_aminus5,outname_a0,outname_a3,outname_supercut);
% cleanup
delete(outname_aminus5,outname_a0,outname_a3);
function [y, index1, index2, frac1, frac2]=coarseGrain(x, flowy, deltaFy, Ny)
%
% coarseGrain --- coarse grain a frequency-series
%
% coarseGrain(x, flowy, deltaFy, Ny) returns a frequency-series structure
% coarse-grained to the frequency values f = flowy + deltaFy*[0:Ny-1].
% coarseGrain(x, flowy, deltaFy, Ny) accepts a frequency-series structure
% (usually a PSD) and returns a frequency-series structure which has been
% coarse-grained to the frequency values f = flowy + deltaFy*[0:Ny-1].
%
% coarseGrain also returns the indices of the lowest and highest frequency
% bins of x that overlap with y (0 <= index1 <= length(x);
% coarseGrain also returns the frequency indices of the lowest and highest
% frequency bins of x that overlap with y (0 <= index1 <= length(x);
% 1 <= index2 <= length(x)+1) and the fractional contributions from these
% frequency bins.
% frequency bins (note that these indices start from 0 and and are 1 less
% than the corresponding Matlab indices). The fractional contribution is
% the fraction of the fine bin that overlaps with any part of the coarse
% series.
%
% Routine written by Joseph D. Romano.
% Contact Joseph.Romano@astro.cf.ac.uk
% The method used is to treat the x-values
%
% $Id: coarseGrain.m,v 1.4 2005-02-24 14:36:27 jromano Exp $
% xc(k) = x.flow + (k - 1)*x.deltaF
%
% as bin centres and the value x.data(k) as the average value over the
% the bin ie.
%
% (1/x.deltaF)*(integral of x.data from the lower edge of the bin to upper edge)
%
% The coarse graining can then be performed by finding the integral
% from the start of the fine series just below the coarse series (via
% the cumulative sum), interpolating to the coarse series, and taking
% differences to recover the average value of the function for each coarse bin.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y, index1, index2, frac1, frac2]=coarseGrain_new(x, flowy, deltaFy, Ny)
% extract metadata
Nx = length(x.data);