From 21c5f04bb36ec90587778a77334c07ce7ab8c8a1 Mon Sep 17 00:00:00 2001
From: Ling Sun <featherlet.sun@gmail.com>
Date: Fri, 29 Jul 2016 08:13:11 +1000
Subject: [PATCH] adjust variable names, minor changes
---
binary_fstat.m => binaryfilter.m | 25 +++++++++++++------------
functions/argmax.m | 4 ++--
readfstats.m | 13 ++++++-------
search_O1.m | 25 +++++++++++++------------
viterbi_colFLT.m | 26 +++++++++++++++-----------
5 files changed, 49 insertions(+), 44 deletions(-)
rename binary_fstat.m => binaryfilter.m (68%)
diff --git a/binary_fstat.m b/binaryfilter.m
similarity index 68%
rename from binary_fstat.m
rename to binaryfilter.m
index bbe75be..9eb2950 100644
--- a/binary_fstat.m
+++ b/binaryfilter.m
@@ -1,5 +1,5 @@
-function [b,h] = binary_fstat(f,f0,P,a0)
-% Generate the Bessel-weighted matched filter (also generate Comb filter)
+function [b,h] = binaryfilter(f,f0,P,a0)
+% BINARYFILTER Generate the Bessel-weighted matched filter (also generate Comb filter)
%
% Inputs:
% f - Frequency bin list (i.e. hidden stats)
@@ -11,23 +11,24 @@ function [b,h] = binary_fstat(f,f0,P,a0)
% 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
-N=ceil(c)+100;
-n=-N:N;
-b0=besselj(n,c);
-M=length(f);
+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
-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
-b=zeros(M,1);
+b = zeros(M,1);
% Place the filter
-b(k)=abs(b0).^2;
+b(k) = abs(b0).^2;
% h is comb matched filter, not used in the search
-h=zeros(M,1);
-h(k)=1;
+h = zeros(M,1);
+h(k) = 1;
diff --git a/functions/argmax.m b/functions/argmax.m
index 309466d..9e0609f 100644
--- a/functions/argmax.m
+++ b/functions/argmax.m
@@ -1,2 +1,2 @@
-function j=argmax(L);
-[~,j]=max(L);
\ No newline at end of file
+function j = argmax(L);
+[~,j] = max(L);
\ No newline at end of file
diff --git a/readfstats.m b/readfstats.m
index 4203546..234ee2c 100644
--- a/readfstats.m
+++ b/readfstats.m
@@ -1,5 +1,5 @@
-function [f,X,N,T]=readfstats(dirName, startFreq, steps)
-% Read Fstat files
+function [f, X, N, T] = readfstats(dirName, startFreq, steps)
+% READFSTATS Read Fstat files
%
% Inputs:
% dirName - Directory of Fstat files
@@ -13,14 +13,13 @@ function [f,X,N,T]=readfstats(dirName, startFreq, steps)
% N - Number of bins
% T - Total number of steps
-d = dir(dirName);
-X = [];
-for n=1:steps
- fileName=['Fstat-' num2str(startFreq) '-' num2str(n-1) '.dat']
+X = [];
+for n = 1:steps
+ fileName = ['Fstat-' num2str(startFreq) '-' num2str(n-1) '.dat']
try
- fstatData=load([dirName fileName]);
+ fstatData = load([dirName fileName]);
X(:,n) = fstatData(:,7);
f = fstatData(:,1);
catch
diff --git a/search_O1.m b/search_O1.m
index 7ab37f1..7d560fc 100644
--- a/search_O1.m
+++ b/search_O1.m
@@ -1,5 +1,5 @@
function search_O1(freq,idx)
-% This script is for O1 Sco X-1 search.
+% SEARCH_O1 This script is for O1 Sco X-1 search.
% 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.
%
@@ -11,6 +11,7 @@ function search_O1(freq,idx)
% Outputs:
% Save output files to the output directory.
+
format long g
startFreq = str2num(freq) + str2num(idx)
@@ -21,7 +22,7 @@ outdir = '/home/ling.sun/O1/Viterbi/output/';
outfile = strcat(outdir,'Viterbi-', num2str(startFreq),'.dat')
% Set parameters
-a = 1.44
+a0median = 1.44
P = 68023.70496
steps = 13
@@ -30,32 +31,32 @@ tic
[f,X]=readfstats(dirName, startFreq, steps);
toc
-fmean=f(1)+0.5; % Take the mean value of each 1-Hz sub-band
-g=[-0.1:0.01:0.1]; % asini grid (+/-10% uncertainty, may search larger uncertainty)
-a0=(1+g)*a; % asini values being searched
+fmean = f(1)+0.5; % Take the mean value of each 1-Hz sub-band
+g = [-0.1:0.01:0.1]; % asini grid (+/-10% uncertainty, may search larger uncertainty)
+a0 = (1+g)*a0median; % asini values being searched
for m=1:length(g)
disp(sprintf('a0=%g',a0(m)))
+
tic
% Generate Bessel-weighted matched filter
- [b,h] = binary_fstat(f,fmean,P,a0(m));
+ [b,h] = binaryfilter(f,fmean,P,a0(m));
% Compute Bessel-weighted Fstat
for n=1:steps
- Y(:,n)=fftshift(ifft(fft(X(:,n)).*fft(b)));
+ Y(:,n) = fftshift(ifft(fft(X(:,n)).*fft(b)));
end
- [path0(m,:),delta,psi,sc(m)] = viterbi_colFLT(3, Y);
-
+ [path0(m,:), delta, psi, sc(m)] = viterbi_colFLT(3, Y);
toc
end
% Get the optimal path for all asini searched
-[maxsc,a0ind]=max(sc)
-path=path0(a0ind,:)
+[maxsc, a0ind] = max(sc)
+path = path0(a0ind,:)
% Save output file for each 1-Hz sub-band
-fid=fopen(outfile,'wt');
+fid = fopen(outfile,'wt');
fprintf(fid,'StartFreq\t%g\n',startFreq);
fprintf(fid,'Score\t%g\n',maxsc);
fprintf(fid,'Bin_a0\t%g\n',a0ind);
diff --git a/viterbi_colFLT.m b/viterbi_colFLT.m
index a10375c..7d85c20 100644
--- a/viterbi_colFLT.m
+++ b/viterbi_colFLT.m
@@ -1,5 +1,5 @@
function [path,delta,psi,score] = viterbi_colFLT(M, obslik)
-% VITERBI 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,
@@ -13,10 +13,13 @@ function [path,delta,psi,score] = viterbi_colFLT(M, obslik)
% 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
-
+
-if round(M/2)==M/2 M=M+1;end
-[Q,T]=size(obslik); % Q bins, T steps
+if round(M/2)==M/2
+ M = M+1;
+end
+
+[Q,T] = size(obslik); % Q bins, T steps
delta = zeros(Q,T);
psi = zeros(Q,T);
@@ -27,23 +30,24 @@ path = zeros(1,T);
% the max or argmax function to this matrix. This function actually
% compares all three possible paths for each step and the operation
% is efficient.
-
-t=1;
+
+cor = (0:Q-1)'-(M-1)/2; % correction term for each bin
+
+t = 1;
delta(:,t) = obslik(:,t);
-cor=(0:Q-1)'-(M-1)/2; % correction term for each bin
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;
+ 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
disp('backward')
% Only identify the optimal path
-[sc,path(T)]=max(delta(:,T));
-score=(sc-mean(delta(:,T)))/std(delta(:,T));
+[sc, path(T)] = max(delta(:,T));
+score = (sc-mean(delta(:,T)))/std(delta(:,T));
for t=T-1:-1:1
path(t) = psi(path(t+1),t+1);
--
GitLab