Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@hweeraso
hweeraso committed May 15, 2024
1 parent 4a40b5b commit 626b696
Show file tree
Hide file tree
Showing 30 changed files with 2,015 additions and 0 deletions.
59 changes: 59 additions & 0 deletions MATLAB Implementation/Figure02_matched_filter.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
clear all
close all
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code visualizes the matched filter
%
% Stanley Chan
% 10/16/2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

option.Display = 'off';

% Define space grid and time grid
T = 10;
dt = 1/256;
t = 0:dt:T-dt;
dx = 1/2048;
x = 0:dx:1-dx;
lengthT = T/dt;
lengthX = 1/dx;

% Define ground truth delay profile
tau = 4./(1+exp(-20*(x-0.5)))+4;
tau_slope = gradient(tau)/dx;

% Define pulse
sigma_t = 0.5; % spread of pulse
alpha = 1000; % per pixel flux
lambda_b = 0; % noise floor

tau0 = 5;
s = (1/sqrt(2*pi*sigma_t^2)) * exp( -(t-tau0).^2/(2*sigma_t^2) );
lambda = alpha*s + lambda_b;

% Monte Carlo simulation
M = poissrnd(alpha+lambda_b*T);
tj = generate_time_stamps(lambda,t,M);

tx_set = randi(lengthT, [1,10]);
lambda_hat = zeros(length(tx_set), lengthT);
for tidx = 1:length(tx_set)
tau = t( tx_set(tidx) );
lambda_hat(tidx, :) = alpha * 1/sqrt(2*pi*sigma_t^2) * exp( -(t-tau).^2/(2*sigma_t^2) ) + lambda_b;
end

lambda0 = alpha * 1/sqrt(2*pi*sigma_t^2) * exp( -(t-mean(tj)).^2/(2*sigma_t^2) ) + lambda_b;

figure;
[num, edges] = histcounts(tj,50);
h1 = bar(edges(1:end-1), num/50, 'EdgeColor', 'none'); hold on;
h2 = plot(t, lambda_hat/max(lambda_hat(:)),'LineWidth', 2,'color',[0.6 0.6 0.6]);
h3 = plot(t, lambda0/max(lambda0(:)),'LineWidth', 4,'color',[0.8 0.4 0.2]);
grid on;
legend([h1(1), h2(1), h3(1)], 'Samples', 'Candidates', 'Max Likelihood Fit', 'Location', 'NE');
xlabel('time, t');
set(gca,'FontWeight','bold','fontsize',14);
set(gcf, 'Position', [100, 100, 600, 300]);

116 changes: 116 additions & 0 deletions MATLAB Implementation/Figure03_space_time.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,116 @@
clear all
close all
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code helps visualize the following three functions
% s(t-tau) = Gaussian(t | tau, sigma_t^2)
% lambda(t) = alpha s(t-tau) + lambda_b
% lambda_tilde = h(x) conv lambda(t)
% lambda_bar = piecewise approx of lambda_tilde
%
% Stanley Chan
% 10/2/2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dx = 1/2048;
x = 0:dx:1-dx;

% %%% Shape 1
% tau = zeros(length(x),1);
% tau(x<0.3) = 5;
% tau( (x>=0.3) & (x<0.7)) = 8;
% tau( x>=0.7) = 6;
% tau_slope = gradient(tau)/dx;
% plot(x,tau);

% %%% Shape 2
% a = [6 0.01 0.55 1.5 1.2];
% tau = a(1)*legendreP(0,x) + a(2)*legendreP(1,x) + a(3)*legendreP(2,x)+ a(4)*legendreP(3,x)+a(5)*legendreP(4,x);
% tau_slope = gradient(tau)/dx;
% % plot(x,tau_slope);

%%% Shape 3
tau = 4./(1+exp(-20*(x-0.5)))+4;
tau_slope = gradient(tau)/dx;
% plot(x,tau);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T = 12;
dt = 1/256;
t = 0:dt:T-dt;
cols = T/dt;
rows = 1/dx;
sigma_t = 0.5;
lambda_b = 0;

N = 16;
dN = rows/N; % number of points per pixel
alpha0 = 16; % overall flux
alpha = alpha0/N; % per pixel flux

s = zeros(rows,cols);
for i=1:rows
s(i,:) = (1/sqrt(2*pi*sigma_t^2)) * exp( -(t-tau(i)).^2/(2*sigma_t^2) );
end
lambda = alpha*s + lambda_b;

figure;
surf(t,x,lambda,'EdgeColor','none');
grid minor;
view(-35,70);
xlabel('time');
ylabel('space');
title('Continous space $\lambda(x,t)$' , 'Interpreter', 'latex')
set(gca,'FontWeight','bold','fontsize',14);

sigma_x = dN*dx;
h = fspecial('gaussian',[6*dN+1, 1], sigma_x/dx);
lambda_tilde = imfilter(lambda, h, 'replicate');
figure;
surf(t,x,lambda_tilde,'EdgeColor','none');
grid minor;
view(-35,70);
xlabel('time');
ylabel('space');
title('Pulse broadening due to oblique surfaces' , 'Interpreter', 'latex')
set(gca,'FontWeight','bold','fontsize',14);


lambda_bar = zeros(rows,cols);
for n=1:N
lambda_bar(dN*(n-1)+[1:dN],:) = repmat(mean( lambda_tilde(dN*(n-1)+[1:dN], :)), [dN, 1]);
end
figure;
surf(t,x,lambda_bar,'EdgeColor','none');
grid minor;
view(-35,70);
xlabel('time');
ylabel('space');
title('Discritized space $\bar{\lambda}(x,t)$' , 'Interpreter', 'latex')
set(gca,'FontWeight','bold','fontsize',14);



% %%%%%%%%%%%%%%%%
% % Incorrect Approximation
% %%%%%%%%%%%%%%%%
% lambda_bar = zeros(rows,cols);
% for n=1:N
% tau_h = imfilter(tau,h,'replicate');
% tau_bar = mean(tau_h(dN*(n-1)+[1:dN]));
% s_bar = (1/sqrt(2*pi*sigma_t^2)) * exp( -(t-tau_bar).^2/(2*sigma_t^2) );
% lambda_bar(dN*(n-1)+[1:dN],:) = repmat(alpha*s_bar+lambda_b,[dN,1]);
% end
%
% figure;
% surf(t,x,lambda_bar,'EdgeColor','none');
% grid minor;
% view(-35,70);
% xlabel('time');
% ylabel('space');
% set(gca,'FontWeight','bold','fontsize',14);
% set(gcf, 'Position', [100, 100, 800, 400]);
% export_fig Figure02_lambda_tau_bar.png -transparent -RGB;
154 changes: 154 additions & 0 deletions MATLAB Implementation/Figure05_1D_mse.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
clear all
close all
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code verifies Stanley's prediction
% (Special Case when lambda_b = 0)
%
% MSE = bias + var
% where
% bias = c^2 / 12N^2
% var = (sigma_x^2 c^2 + sigma_t^2)N/alpha0
%
% Remark: If lambda_b != 0, then the estimate
% is not a simple average.
%
% Stanley Chan
% 10/4/2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nset = [4 8 16 32 64 128 256];

mse_num = zeros(1,length(Nset));
mse_theory = zeros(1,length(Nset));
var_num = zeros(1,length(Nset));
var_theory = zeros(1,length(Nset));
bias_num = zeros(1,length(Nset));
bias_theory = zeros(1,length(Nset));

for nidx = 1:length(Nset)
fprintf('nidx = %3g \n', nidx);
% Define space grid and time grid
T = 10;
dt = 1/256;
t = 0:dt:T-dt;
dx = 1/2048;
x = 0:dx:1-dx;
lengthT = T/dt;
lengthX = 1/dx;

% Define ground truth delay profile
tau = 4./(1+exp(-20*(x-0.5)))+4;
tau_slope = gradient(tau)/dx;

% Define pulse
sigma_t = 0.5; % spread of pulse
N = Nset(nidx); % number of pixels
dN = lengthX/N; % number of points per pixel
alpha0 = 10000; % overall flux
alpha = alpha0/N; % per pixel flux
lambda_b = 0; % noise floor
s = zeros(lengthX,lengthT);
for i=1:lengthX
s(i,:) = (1/sqrt(2*pi*sigma_t^2)) * exp( -(t-tau(i)).^2/(2*sigma_t^2) );
end
lambda = alpha*s + lambda_b;

% Define spatial spread
sigma_x = dN*dx/sqrt(12);
h = ones(dN,1)/dN;
lambda_tilde = imfilter(lambda, h, 'replicate');

% Convolved pulse lambda(x,t)
lambda_bar_n = zeros(N,lengthT); % coefficients
lambda_bar = zeros(lengthX,lengthT); % full function
for n=1:N
lambda_bar_n(n,:) = lambda_tilde(dN*(n-1)+dN/2,:);
lambda_bar(dN*(n-1)+[1:dN],:) = repmat(lambda_bar_n(n,:),[dN,1]);
end

% Approximated delay function tau(x)
taubar_n = zeros(1,N); % coefficients
taubar = zeros(1,lengthX); % full function
for n=1:N
taubar_n(n) = mean(tau(dN*(n-1)+[1:dN]));
taubar(dN*(n-1)+[1:dN]) = taubar_n(n)*ones(1,dN);
end

% Monte Carlo simulation
max_trial = 100;
tauhat_n = zeros(max_trial,N);
tauhat = zeros(max_trial,lengthX);
for trial = 1:max_trial
for n=1:N
M = poissrnd(alpha+lambda_b*T);
tj = generate_time_stamps(lambda_bar_n(n,:),t,M);
tauhat_n(trial,n) = mean(tj);
tauhat(trial, dN*(n-1)+[1:dN]) = tauhat_n(trial,n)*ones(1,dN);
end
end

% Numerical Values
bias_num(nidx) = sum((taubar - tau).^2*dx);
var_num_trial = zeros(1, max_trial);
mse_num_trial = zeros(1, max_trial);
for trial = 1:max_trial
var_num_trial(trial) = mean( (tauhat(trial,:)-taubar).^2 );
mse_num_trial(trial) = mean( (tauhat(trial,:)-tau).^2 );
end
var_num(nidx) = mean(var_num_trial);
mse_num(nidx) = mean(mse_num_trial);
end

Nset1 = 4:4:256;
for nidx = 1:length(Nset1)
N = Nset1(nidx);
dN = floor(lengthX/N);
rows1 = N*dN;
c = zeros(1,N);
for n=1:N
c(n) = mean( tau_slope(dN*(n-1)+[1:dN]) );
end
var_theory(nidx) = mean(( c.^2 * sigma_x^2 + sigma_t^2 )*N/alpha0);
bias_theory(nidx) = mean((c.^2)/(12*N^2));
mse_theory(nidx) = bias_theory(nidx) + var_theory(nidx);

var_theory2(nidx) = mean(( sigma_t^2 )*N/alpha0);
bias_theory2(nidx) = mean((c.^2)/(12*N^2));
mse_theory2(nidx) = bias_theory2(nidx) + var_theory2(nidx);
end

% figure;
% semilogy(Nset, var_num,'o', 'LineWidth', 2); hold on;
% semilogy(Nset1, var_theory, 'LineWidth', 2);
% grid on;
% legend('numerical', 'theory');
% xlabel('N, number of pixels');
% ylabel('variance');
% ylabel('Variance');
% set(gca,'FontWeight','bold','fontsize',14);
% set(gcf, 'Position', [100, 100, 600, 400]);
% %
% figure;
% semilogy(Nset, bias_num,'o', 'LineWidth', 2); hold on;
% semilogy(Nset1, bias_theory, 'LineWidth', 2);
% grid on;
% legend('numerical', 'theory');
% xlabel('N, number of pixels');
% ylabel('bias');
% ylabel('Bias');
% set(gca,'FontWeight','bold','fontsize',14);
% set(gcf, 'Position', [100, 100, 600, 400]);

figure;
area(Nset1, mse_theory, 'FaceColor',[0.6 0.6 0.6]); hold on;
set (gca, 'Yscale', 'log');
semilogy(Nset, mse_num,'o', 'LineWidth', 4, 'Color', [0.3 0.4 0.9]);
semilogy(Nset1, mse_theory, 'LineWidth', 4, 'Color', [0.8 0.4 0.2]);
grid on;
% legend('Infeasible', 'Simulated MSE (1D)', 'Theoretical MSE (1D)');
xlabel('N, number of pixels');
ylabel('MSE');
set(gca,'FontWeight','bold','fontsize',14);

0 comments on commit 626b696

Please sign in to comment.