Commit b380dd5f authored by Ling Sun's avatar Ling Sun
Browse files

add more comments

parent d92d1364
function [b,h] = binary_fstat(f,f0,P,a0)
% Phase - Ph = f0*(t+a0*sin(2*pi*t/P));
% Frequency - F = f0+2*pi*a0*f0/P*cos(2*pi*t/P);
% Fstat width = [-c c]/P where c --- frequency magnitude (Carson's rule)
% Generate the Bessel-weighted matched filter (also generate Comb filter)
%
% 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)
% 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
......
% Output: Frequency bins, F-stat values, number of bins, total steps
% Fstat name format - Fstat-<start freq>-<step index>.dat
function [f,X,N,T]=readfstats(dirName, startFreq, steps)
% Read Fstat files
%
% 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
%
% Outputs:
% f - Frequency bins
% X - F-stat values
% N - Number of bins
% T - Total number of steps
d = dir(dirName);
X = [];
......
function search_O1(freq,idx)
% 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.
%
% Inputs:
% freq - Start frequency for all the parallel jobs
% idx - Job index
% Start frequency for each sub-job: freq + idx
%
% Outputs:
% Save output files to the output directory.
format long g
startFreq = str2num(freq) + str2num(idx)
% Set directories and file name
dirName = '/home/ling.sun/O1/Viterbi/fstat/';
outdir = '/home/ling.sun/O1/Viterbi/output/';
outfile = strcat(outdir,'Viterbi-', num2str(startFreq),'.dat')
% Set parameters
a = 1.44
P = 68023.70496
steps = 13
% Read Fstats
tic
[f,X]=readfstats(dirName, startFreq, steps); % Fstats
[f,X]=readfstats(dirName, startFreq, steps);
toc
fmean=f(1)+0.5;
g=[-0.1:0.01:0.1]; % asini grid
a0=(1+g)*a;
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
for m=1:21
tic
[b,h] = binary_fstat(f,fmean,P,a0(m)); % Bessel-weighted matched filter
% Generate Bessel-weighted matched filter
[b,h] = binary_fstat(f,fmean,P,a0(m));
% Compute Bessel-weighted Fstat
for n=1:steps
Y(:,n)=fftshift(ifft(fft(X(:,n)).*fft(b))); % Compute Bessel-weighted Fstat
Y(:,n)=fftshift(ifft(fft(X(:,n)).*fft(b)));
end
[path0(m,:),delta,psi,sc(m)] = viterbi_colFLT(3, Y);
%D(:,m)=delta(:,end);
toc
end
......
Markdown is supported
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