Commit 21c5f04b authored by Ling Sun's avatar Ling Sun
Browse files

adjust variable names, minor changes

parent 72f727e6
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;
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
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
......
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);
......
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);
......
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