Commit 1081f013 authored by Ling Sun's avatar Ling Sun
Browse files

modify descriptions and comments

parent 21c5f04b
function j = argmax(L)
% ARGMAX Index j of the maximum value in L.
[~,j] = max(L);
\ No newline at end of file
function [b,h] = binaryfilter(f,f0,P,a0) 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: % Output:
% f - Frequency bin list (i.e. hidden stats) % b - Bessel-weighted matched filter
% f0 - Mean frequency for each sub-band % h - Comb matched filter (not used in the current search)
% 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)
% 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; c = 2*pi*a0*f0;
% Make it slightly wider than c % Take a wider range [-N N] for calculation
N = ceil(c)+100; N = ceil(c)+100;
n = -N:N; n = -N:N;
b0 = besselj(n,c); b0 = besselj(n,c);
M = length(f); 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); 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); b = zeros(M,1);
% Place the filter % Place the filter
b(k) = abs(b0).^2; 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 = zeros(M,1);
h(k) = 1; h(k) = 1;
function j = argmax(L);
[~,j] = max(L);
\ No newline at end of file
function [f, X, N, T] = readfstats(dirName, startFreq, steps) 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: % I/O Spec
% dirName - Directory of Fstat files % ========
% startFreq - Start frequency of the Fstat % Input:
% steps - Total number of steps % dirName - Directory of Fstat files
% Fstat file name should follow the format: Fstat-<start freq>-<step index>.dat % startFreq - Start frequency marking the Fstat file name
% steps - Total number of steps
% %
% Outputs: % Output:
% f - Frequency bins % f - Frequency bins
% X - F-stat values % X - F-stat values
% N - Number of bins % N - Total number of bins
% T - Total number of steps % T - Total number of steps
X = []; X = [];
...@@ -23,7 +30,7 @@ for n = 1:steps ...@@ -23,7 +30,7 @@ for n = 1:steps
X(:,n) = fstatData(:,7); X(:,n) = fstatData(:,7);
f = fstatData(:,1); f = fstatData(:,1);
catch catch
disp(n) disp(sprintf('Failed to read step %g',n))
end end
end end
......
...@@ -3,14 +3,15 @@ function search_O1(freq,idx) ...@@ -3,14 +3,15 @@ function search_O1(freq,idx)
% The search is separated into parallel 1-Hz sub-jobs. % 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. % The start frequency for each sub-job is the overal start frequency plus job index.
% %
% Inputs: % I/O Spec
% freq - Start frequency for all the parallel jobs % ========
% idx - Job index % Input:
% Start frequency for each sub-job: freq + idx % freq - Start frequency for all the parallel jobs
% idx - Condor job index
% Note: Start frequency for each sub-job: freq + idx
% %
% Outputs: % Output:
% Save output files to the output directory. % Save output files to the output directory.
format long g format long g
...@@ -28,7 +29,7 @@ steps = 13 ...@@ -28,7 +29,7 @@ steps = 13
% Read Fstats % Read Fstats
tic tic
[f,X]=readfstats(dirName, startFreq, steps); [f,X] = readfstats(dirName, startFreq, steps);
toc toc
fmean = f(1)+0.5; % Take the mean value of each 1-Hz sub-band fmean = f(1)+0.5; % Take the mean value of each 1-Hz sub-band
...@@ -46,7 +47,7 @@ for m=1:length(g) ...@@ -46,7 +47,7 @@ for m=1:length(g)
for n=1:steps for n=1:steps
Y(:,n) = fftshift(ifft(fft(X(:,n)).*fft(b))); Y(:,n) = fftshift(ifft(fft(X(:,n)).*fft(b)));
end end
% Viterbi tracking
[path0(m,:), delta, psi, sc(m)] = viterbi_colFLT(3, Y); [path0(m,:), delta, psi, sc(m)] = viterbi_colFLT(3, Y);
toc toc
end end
......
function [path,delta,psi,score] = viterbi_colFLT(M, obslik) 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: % I/O Spec
% M - One dimension of the colfilt matrix, in the current model M=3, % ========
% but keep it as an input for generic purpose % Input:
% obslik - Observation likelihood, i.e. Fstat % 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: % Output:
% path(n,t) - The nth best path, where path(n,1), path(n,2) ... path(n,T) stands for % path(t) - The optimal path, where path(1), path(2), ..., path(T) stands
% the frequency bin index for each step in the nth best path. % for the frequency bin index for each step.
% delta(j,t) - Max probability of the sequence of length t-1 and then going to state j % delta(j,t) - Max probability of the sequence till step t-1
% psi(j,t) - The best predecessor state, given that we ended up in state j at t % and then go to state j
% score - Number of standard deviations above the mean value of log likelihood % psi(j,t) - The best predecessor state, given that it ends up in state j at t
% of paths to all the states at the final step % 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 if round(M/2)==M/2
M = M+1; M = M+1;
...@@ -31,21 +34,22 @@ path = zeros(1,T); ...@@ -31,21 +34,22 @@ path = zeros(1,T);
% compares all three possible paths for each step and the operation % compares all three possible paths for each step and the operation
% is efficient. % 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; t = 1;
delta(:,t) = obslik(:,t); delta(:,t) = obslik(:,t);
% Recursion and Termination
disp('forward') disp('forward')
for t=2:T for t=2:T
delta(:,t) = colfilt(delta(:,t-1),[M 1],'sliding',@max) + obslik(:,t); delta(:,t) = colfilt(delta(:,t-1),[M 1],'sliding',@max) + obslik(:,t);
psi(:,t) = colfilt(delta(:,t-1),[M 1],'sliding',@argmax) + cor; psi(:,t) = colfilt(delta(:,t-1),[M 1],'sliding',@argmax) + cor;
end end
% Sort the paths % Optimal path backtracking
disp('backward') disp('backward')
% Only identify the optimal path
[sc, path(T)] = max(delta(:,T)); [sc, path(T)] = max(delta(:,T));
score = (sc-mean(delta(:,T)))/std(delta(:,T)); score = (sc-mean(delta(:,T)))/std(delta(:,T));
......
Supports Markdown
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