diff --git a/argmax.m b/argmax.m
new file mode 100644
index 0000000000000000000000000000000000000000..740bd40286e42f3d602222ffeec7fce94b78b70b
--- /dev/null
+++ b/argmax.m
@@ -0,0 +1,4 @@
+function j = argmax(L)
+% ARGMAX Index j of the maximum value in L.
+
+[~,j] = max(L);
\ No newline at end of file
diff --git a/binaryfilter.m b/binaryfilter.m
index 9eb29509777e8cd16dc95679790341b0faa6385f..97111096bc67ac1130c6e245f1b8e7d47409df86 100644
--- a/binaryfilter.m
+++ b/binaryfilter.m
@@ -1,34 +1,41 @@
function [b,h] = binaryfilter(f,f0,P,a0)
-% BINARYFILTER Generate the Bessel-weighted matched filter (also generate Comb filter)
+% BINARYFILTER Generate Bessel-weighted matched filter
+% (also generate Comb filter)
+%
+% Note
+% ----
+% Fstat power is concentrated at [-c c]/P where c = 2*pi*a0*f0.
+% Use fixed f0 (mean value of the sub-band) to setup matched filter.
+%
+% I/O Spec
+% ========
+% Input:
+% f - Frequency bin list (i.e. hidden stats)
+% f0 - Mean frequency for each sub-band
+% P - Orbit period
+% a0 - Projected orbital semimajor axis (i.e. asini)
%
-% Inputs:
-% f - Frequency bin list (i.e. hidden stats)
-% f0 - Mean frequency for each sub-band
-% P - Orbit period
-% a0 - Projected orbital semimajor axis (i.e. asini)
-%
-% Outputs:
-% b - Bessel-weighted matched filter
-% h - Comb filter (not used in the current search)
-
+% Output:
+% b - Bessel-weighted matched filter
+% h - Comb matched filter (not used in the current search)
-% Fstat power concentrated at [-c c]/P where c = 2*pi*a0*f0
-% f0 is fixed to setup matched filter for each 1-Hz sub-band
c = 2*pi*a0*f0;
-% Make it slightly wider than c
+% Take a wider range [-N N] for calculation
N = ceil(c)+100;
n = -N:N;
b0 = besselj(n,c);
M = length(f);
-% Identify the indices to place the matched filter in the 1-Hz band
-% k indicates the interpolated indices at query points n/P using nearest interpolation, 1 is the extrapval for out of range values
+
+% Identify the indices to place the matched filter in the sub-band.
+% k indicates the interpolated indices at query points n/P,
+% using nearest interpolation. Set extrapval to 1 for out of range values.
k = interp1(f-mean(f),1:M,n/P,'nearest',1);
-% Initialise the 1-Hz band matched filter
+% Initialise the matched filter for the whole sub-band
b = zeros(M,1);
% Place the filter
b(k) = abs(b0).^2;
-% h is comb matched filter, not used in the search
+% Set comb matched filter, which can be used for a C-stat search
h = zeros(M,1);
h(k) = 1;
diff --git a/functions/argmax.m b/functions/argmax.m
deleted file mode 100644
index 9e0609f3b6608ce6f3e81f497efcf22f92229f19..0000000000000000000000000000000000000000
--- a/functions/argmax.m
+++ /dev/null
@@ -1,2 +0,0 @@
-function j = argmax(L);
-[~,j] = max(L);
\ No newline at end of file
diff --git a/readfstats.m b/readfstats.m
index 234ee2cdbf44c44dde92b0c1a7b119caee9bc10c..321cf9d29be2bc0e786bcd427eba3429c5ecfbce 100644
--- a/readfstats.m
+++ b/readfstats.m
@@ -1,17 +1,24 @@
function [f, X, N, T] = readfstats(dirName, startFreq, steps)
-% READFSTATS Read Fstat files
+% READFSTATS Read Fstat files for all steps
+%
+% Note
+% ----
+% Fstat file name should follow the format:
+% Fstat-<start freq>-<step index>.dat
+% <step index> start from 0.
%
-% Inputs:
-% dirName - Directory of Fstat files
-% startFreq - Start frequency of the Fstat
-% steps - Total number of steps
-% Fstat file name should follow the format: Fstat-<start freq>-<step index>.dat
+% I/O Spec
+% ========
+% Input:
+% dirName - Directory of Fstat files
+% startFreq - Start frequency marking the Fstat file name
+% steps - Total number of steps
%
-% Outputs:
-% f - Frequency bins
-% X - F-stat values
-% N - Number of bins
-% T - Total number of steps
+% Output:
+% f - Frequency bins
+% X - F-stat values
+% N - Total number of bins
+% T - Total number of steps
X = [];
@@ -23,7 +30,7 @@ for n = 1:steps
X(:,n) = fstatData(:,7);
f = fstatData(:,1);
catch
- disp(n)
+ disp(sprintf('Failed to read step %g',n))
end
end
diff --git a/search_O1.m b/search_O1.m
index 7d560fc5863b8e33b07e43bc5726e8a054c1e6a2..4456b61c85bc7bf2ee3876a49e4e09de98b1adf9 100644
--- a/search_O1.m
+++ b/search_O1.m
@@ -3,14 +3,15 @@ function search_O1(freq,idx)
% The search is separated into parallel 1-Hz sub-jobs.
% The start frequency for each sub-job is the overal start frequency plus job index.
%
-% Inputs:
-% freq - Start frequency for all the parallel jobs
-% idx - Job index
-% Start frequency for each sub-job: freq + idx
+% I/O Spec
+% ========
+% Input:
+% freq - Start frequency for all the parallel jobs
+% idx - Condor job index
+% Note: Start frequency for each sub-job: freq + idx
%
-% Outputs:
-% Save output files to the output directory.
-
+% Output:
+% Save output files to the output directory.
format long g
@@ -28,7 +29,7 @@ steps = 13
% Read Fstats
tic
-[f,X]=readfstats(dirName, startFreq, steps);
+[f,X] = readfstats(dirName, startFreq, steps);
toc
fmean = f(1)+0.5; % Take the mean value of each 1-Hz sub-band
@@ -46,7 +47,7 @@ for m=1:length(g)
for n=1:steps
Y(:,n) = fftshift(ifft(fft(X(:,n)).*fft(b)));
end
-
+ % Viterbi tracking
[path0(m,:), delta, psi, sc(m)] = viterbi_colFLT(3, Y);
toc
end
diff --git a/viterbi_colFLT.m b/viterbi_colFLT.m
index 7d85c207643a37ebdb628ceb0c83a4cd3780d70a..6a56ec1619c53b275da129c7f78094719aa2abbd 100644
--- a/viterbi_colFLT.m
+++ b/viterbi_colFLT.m
@@ -1,19 +1,22 @@
function [path,delta,psi,score] = viterbi_colFLT(M, obslik)
-% VITERBI_COLFLT Find the most-probable (Viterbi) path through the HMM states (frequency bins).
+% VITERBI_COLFLT Find the most-probable (Viterbi) path
+% through the HMM states (frequency bins).
%
-% Inputs:
-% M - One dimension of the colfilt matrix, in the current model M=3,
-% but keep it as an input for generic purpose
-% obslik - Observation likelihood, i.e. Fstat
+% I/O Spec
+% ========
+% Input:
+% M - One dimension of the colfilt matrix, in the current model M=3.
+% Keep it as an input for generic purpose.
+% obslik - Observation likelihood, i.e. Fstat
%
-% Outputs:
-% path(n,t) - The nth best path, where path(n,1), path(n,2) ... path(n,T) stands for
-% the frequency bin index for each step in the nth best path.
-% delta(j,t) - Max probability of the sequence of length t-1 and then going to state j
-% psi(j,t) - The best predecessor state, given that we ended up in state j at t
-% score - Number of standard deviations above the mean value of log likelihood
-% of paths to all the states at the final step
-
+% Output:
+% path(t) - The optimal path, where path(1), path(2), ..., path(T) stands
+% for the frequency bin index for each step.
+% delta(j,t) - Max probability of the sequence till step t-1
+% and then go to state j
+% psi(j,t) - The best predecessor state, given that it ends up in state j at t
+% score - Number of standard deviations above the mean value of log
+% likelihood of paths to all the states at the final step
if round(M/2)==M/2
M = M+1;
@@ -31,21 +34,22 @@ path = zeros(1,T);
% compares all three possible paths for each step and the operation
% is efficient.
-cor = (0:Q-1)'-(M-1)/2; % correction term for each bin
+% Set correction term for each bin
+cor = (0:Q-1)'-(M-1)/2;
+% Initialisation
t = 1;
delta(:,t) = obslik(:,t);
+% Recursion and Termination
disp('forward')
for t=2:T
delta(:,t) = colfilt(delta(:,t-1),[M 1],'sliding',@max) + obslik(:,t);
psi(:,t) = colfilt(delta(:,t-1),[M 1],'sliding',@argmax) + cor;
end
-% Sort the paths
+% Optimal path backtracking
disp('backward')
-
-% Only identify the optimal path
[sc, path(T)] = max(delta(:,T));
score = (sc-mean(delta(:,T)))/std(delta(:,T));