【粒子跟踪】现有代码解析（2）：track

track_XYBayesZBayes.mhttps://github.com/welsherlab/3DTrackSim/blob/main/track_XYBayesZBayes.m

一般情况。要设定均匀背景，需将 r 设为空矩阵 [] ，并在 bofr 中仅输入一个值。一般情况下的跟踪具有真实的阶段响应、不均匀的背景强度以及在跟踪过程中对活动对象的估计。固定 kt 扫描模式。采样时间 tau 限制为完整的 TAG 周期。将采样时间限制为 1.5 个 TAG 周期可能会导致像素高度采样出现不对称性，如果需要分别沿 xy 和 z 方向预测背景的话（可能需要这样做），这种对高和低像素的扫描可能会有用。

并非真的 1D 模型，必须使用来自 2D 贝叶斯查找表的中间值来对查找表进行沿 z 轴的缩放，否则代码将无法正常运行。使用真实的舞台响应和不均匀的背景强度来模拟 2D 跟踪。

函数定义

function [trkerr,len,tau,mleerr,photons,stg_out,part_out,posest_out,N,sbandests,aoe]=track_XYBayesZBayes(D,s,bofr,r,ki,kiz,N,ogtau,varargin)

输入参数：
- D：粒子的扩散系数，单位为 m²/s。
- s：粒子的亮度，单位为计数每秒（cps）。
- bofr：背景强度的 1D 向量，长度比 r 的长度多 1。背景强度根据 r 的值选择。
- r：背景强度变化的 1D 向量。
- ki：xy 轴的积分反馈参数。
- kiz：z 轴的积分反馈参数。
- N：轨迹的初始步数。
- ogtau：初始的时间间隔。
- varargin：可选参数。

输出参数：
- trkerr：粒子位置与舞台位置之间的绝对差异（xy、z、xyz）。
- len：轨迹的长度（步数）。
- tau：经过调整后的最终时间间隔。
- mleerr：粒子位置与估计位置之间的绝对差异（xy、z、xyz）。
- photons：每个时间间隔内观察到的光子数。
- stg_out：舞台位置（x、y、z）/激光扫描中心。
- part_out：粒子位置（x、y、z）。
- posest_out：相对于舞台中心的位置估计（x、y、z）。
- N：轨迹中的最大步数。
- sbandests：整个轨迹的信号、背景、信号估计和背景估计。
- aoe：逃逸轴编码。1 表示 x 轴，2 表示 y 轴，3 表示 z 轴，4 表示粒子达到完整持续时间。
例子：参数初始化

D = 6e-12; % 扩散系数 (m^2/s)
s = 5e6; % 粒子亮度 (cps)
bofr = [0, 1e6]; % 背景强度 (cps)
r = 2e-6; % 不均匀背景，背景强度在 2 微米处发生变化
ki = 0.05; % xy 轴的积分反馈参数
kiz = 0.04; % z 轴的积分反馈参数
N = 5e4; % 轨迹步数
ogtau = 20e-6; % 时间间隔 (s)

代码逻辑

输入参数的处理：
- 如果 r 是空矩阵 []，则表示背景是均匀的。
- bofr 是一个 1D 向量，其长度比 r 的长度多 1，用于定义不同位置的背景强度。
跟踪过程：
- 函数模拟粒子在 xy 平面上的运动，并考虑 z 轴上的运动。
- 使用贝叶斯方法（2D Bayes LUT）来估计粒子的位置。
- 考虑了现实平移台响应和不均匀背景强度的影响。
- 在跟踪过程中，动态调整时间间隔 tau，以确保跟踪的准确性和效率。
输出结果：
- trkerr 和 mleerr 分别表示粒子位置与平移台位置、估计位置之间的误差。
- photons 记录了每个时间间隔内观察到的光子数。
- stg_out 和 part_out 分别记录了平移台位置和粒子位置。
- posest_out 记录了相对于平移台中心的位置估计。
- sbandests 提供了信号、背景、信号估计和背景估计的详细信息。
- aoe 编码了粒子逃逸的轴向信息。

注意事项

时间间隔 tau 的限制：
- 时间间隔 tau 被限制为完整的 TAG 周期，以避免在像素高度采样中引入不对称性。
背景强度的处理：
- 背景强度在 r 指定的位置发生变化，这会影响粒子的跟踪过程。
贝叶斯方法的应用：
- 使用 2D Bayes LUT 来估计粒子的位置，需要确保中间值的正确使用，否则代码可能无法正常工作。

输入参数处理

%signal background estimation flag
p = inputParser;
addParameter(p,'sbest',true,@islogical)%optional input to toggle signal and bg estimation off. Defaults to true.
addParameter(p,'send',s,@isnumeric)%optional input to vary signal intensity over course of a traj using exponential decay. 
addParameter(p,'Dest',D,@isnumeric)%enables passing of misleading D for use in estimation while still generating trajectory with true diffusion coefficient. Defaults to truth. 
parse(p,varargin{:})

inputParser：用于解析输入参数。
addParameter：添加可选参数及其默认值和验证函数。
- 'sbest'：信号和背景估计标志，默认为 true。
- 'send'：允许在轨迹过程中使用指数衰减变化信号强度，默认为 s。
- 'Dest'：允许传递误导性的扩散系数 D，用于估计，而实际轨迹仍然使用真实的扩散系数，默认为 D。
parse：解析输入参数，确保它们符合预期。

输入验证

if p.Results.send > s
    error('Why would intensity increase over the course of a traj? Submit less stupid inputs.')
end
Dest = p.Results.Dest;

检查 send 是否大于 s，如果大于，则抛出错误，因为信号强度不应在轨迹过程中增加。
将 Dest 设置为解析后的值。

TAG 和 EOD 扫描参数

%TAG and EOD Scan Parameters.
d = 500e-9;       %Size of KT
f = 70e3;         %TAG lens frequency
amp = 1e-6;       %TAG max scan amplitude (amplitude 1 micron total scan size 2)
binspertag = 24;  %this makes subtau/z binning 2x finer than the bins in xiaochen's existing sim which has 14 bins per xy scan pos.

d：KT 的大小（500 纳米）。
f：TAG 镜头频率（70 kHz）。
amp：TAG 最大扫描振幅（1 微米，总扫描尺寸 2 微米）。
binspertag：每个 TAG 周期的 bin 数，这里设置为 24，比现有的模拟更精细。

输入检查

%check inputs
if (length(bofr)-1)~=length(r)
    error('input dimensions of r and background at given r are incompatible')
end
if length(r)>1
    if any(diff(r)<=0)
        error('radius of background flips must increase monotonically to avoid ambiguous definitions')
    end
end
if N>=2^31-1
    error('trajectory duration exceeds step numbers numerical precision, increase precision of N or shorten trajectory')
end

检查 bofr 和 r 的长度是否匹配，确保背景强度的定义与 r 的长度一致。
如果 r 的长度大于 1，检查 r 是否单调递增，避免背景强度的定义出现歧义。
检查 N 是否超过数值精度限制（2^31-1），如果超过，则抛出错误，建议增加 N 的精度或缩短轨迹。

处理空的 `r`

%if no radius for flip define arbitrary radius r and same background on
%both sides of flip
if isempty(r)
    bofr=[bofr,bofr];
    r=1;
end

如果 r 是空的，定义一个任意的半径 r，并设置相同的背景强度 bofr 在 r 的两侧。

调整时间间隔 τ和网格大小 gs

时间间隔 τ 的调整

%Round tau to nearest TAG period which introduces reproducible photon
%indexing/binning. 
tp = 1 / f; % TAG period
tau = tp * 1 * round(ogtau / tp / 1);

tp：TAG 周期，计算为 1/f。
tau：将初始时间间隔 τ 调整为最接近 TAG 周期的整数倍。

时间间隔 τ 的最小值限制

if tau <= 20e-6 % computation cannot compile/run in less than 20 so min bin size is 20
    tau = tp * 1 * ceil(ogtau / tp / 1);
end

这段代码主要负责调整时间间隔 τ 和网格大小 gs，以确保模拟的准确性和计算的可行性。以下是对代码的详细解释：

时间间隔 τ 的调整

%Round tau to nearest TAG period which introduces reproducible photon
%indexing/binning. 
tp = 1 / f; % TAG period
tau = tp * 1 * round(ogtau / tp / 1);

tp：TAG 周期，计算为 f1。
tau：将初始时间间隔 τ 调整为最接近 TAG 周期的整数倍。

时间间隔 τ 的最小值限制

if tau <= 20e-6 % computation cannot compile/run in less than 20 so min bin size is 20
    tau = tp * 1 * ceil(ogtau / tp / 1);
end

如果调整后的 τ 小于 20 微秒，将其向上取整到最接近的 TAG 周期。

网格大小 gs 的初步定义

gsmax = 41; % Maximum grid size for FPGA implementation
if Dest ~= 0
    gs = ceil(1e-6 / 1.5 / sqrt(2 * Dest * tau));
else
    gs = gsmax;
end

gsmax：FPGA 实现的最大网格大小。
如果扩散系数 D 不为零，计算网格大小 gs，否则直接设置为 gsmax。

确保激光扫描位置落在网格上

% ensure laser scan positions fall exactly onto grid by fixing grid size to
% nearest multiple of 5
if rem(gs - 5, 4) == 0
elseif rem(gs - 5, 4) == 1
    gs = gs - 1;
elseif rem(gs - 5, 4) == 2
    gs = gs + 2;
elseif rem(gs - 5, 4) == 3
    gs = gs + 1;
end

调整 gs 使其为 5 的倍数，确保激光扫描位置恰好落在网格上。

重新计算 τ 以适应 FPGA 实现

% limit computation duration to accurate FPGA implementation max grid size by recalculating tau to the nearest us
if gs > gsmax
    % recalculate tau rounded to the nearest ms
    tau = round((1e-6 / 1.5 / gsmax)^2 / 2 / Dest, 6);
    % Define finalized modified grid size if necessary
    gs = ceil(1e-6 / 1.5 / sqrt(2 * Dest * tau));
    % ensure laser scan positions fall exactly onto grid
    if rem(gs - 5, 4) == 0
    elseif rem(gs - 5, 4) == 1
        gs = gs - 1;
    elseif rem(gs - 5,4 ) == 2
        gs = gs + 2;
    elseif rem(gs - 5, 4) == 3
        gs = gs + 1;
    end
else
    % do not change tau from the rounded to nearest TAG period value
end

如果计算出的 gs 超过 FPGA 实现的最大网格大小 gsmax，重新计算 τ。
重新计算 gs 并确保激光扫描位置落在网格上。

错误处理

if gs > gsmax
    error('correction by changing tau failed so the round vs ceil thing did matter rethink this implementation')
end
if gs < 5
    error('particle motion is too fast to be tracked without expanding scan area')
end

如果重新计算后的 gs 仍然超过 gsmax，抛出错误。
如果 gs 小于 5，表示粒子运动过快，无法在当前扫描区域内跟踪，抛出错误。

确定最终的时间间隔 τ

% finalize tau by rounding to nearest TAG Period
tau = tp * 1 * round(tau / tp / 1);

将 τ 调整为最接近 TAG 周期的整数倍。

调整轨迹步数 N

if tau ~= ogtau
    % recalculate N preserve to intended trajectory duration.
    N = int32(ceil(N * ogtau / tau));
    if tau < ogtau
%         disp(['Tau has been decreased to ' num2str(tau) ' N has been increased to ' num2str(N) ' to preserve intended trajectory duration'])
    else
%         disp(['Tau has been increased to ' num2str(tau) ' N has been decreased to ' num2str(N) ' to preserve intended trajectory duration'])
    end
end

如果 τ 发生变化，重新计算轨迹步数 N，以保持轨迹的总持续时间不变。

信号强度的定义

%after N has been recaulculated define signal intensity for traj.
if s == p.Results.send
    s(1:N) = s;
else
    s(1:N) = s * ((p.Results.send / s).^(1 / double(N))).^(double(1:N));
end

如果 s 等于 p.Results.send，则信号强度在整个轨迹中保持不变。
如果 s 不等于 p.Results.send，则信号强度在轨迹中按照指数衰减变化。具体来说，信号强度从初始值 s 逐渐衰减到 p.Results.send，衰减过程在 N 步内完成。

计算每个 EOD bin 的 bin 数

binsperz = round(binspertag * tau / tp);

计算每个 EOD bin 的 bin 数，确保每个 TAG 周期内的 bin 数与时间间隔 τ 成正比。

初始化光子生成的参数

tagtau = tau / binsperz;
b = zeros(1, binsperz); %initialize b values for photon generation.

tagtau：每个 bin 的时间间隔。
b：初始化用于光子生成的参数数组。

生成 XY 和 Z 坐标网格

%Generate XY and Z coordinates, granularity along z is 2 x gs of xy to make
%voxels square making D universal instead of axially dependent.
coord = linspace(-0.5e-6, 0.5e-6, gs);
[x2D, y2D] = meshgrid(coord);
z1D = linspace(-amp, amp, 2 * gs - 1); %z axis scan is 2x xy so equivalent grid spacing will need 2x pts.

coord：生成 XY 平面上的坐标网格，范围从 -0.5 微米到 0.5 微米，共 gs 个点。
x2D 和 y2D：使用 meshgrid 生成二维坐标网格。
z1D：生成 Z 轴上的坐标，范围从 -amp 到 amp，共 2 * gs - 1 个点。Z 轴的分辨率是 XY 平面的两倍，以确保体素（voxels）是正方形，从而使扩散系数 D 在所有轴上都是通用的。

质量控制：检查 TAG 采样

%Quality control bins per z. sampling must have at least 8 pts/sine wave
%and meet or exceed bins along z axis.
if binsperz < tau / tp * 8
    error('TAG sampling is inadequate to reproduce properties of sine wave, TAU must be subdivided into a minimum of 8 points per period to ensure adequate sampling, increase points per subbin')
end

检查每个 TAG 周期内的 bin 数是否足够，以确保能够准确再现正弦波的特性。
如果 binsperz 小于 τ/tp×8，则抛出错误，提示需要增加每个子 bin 的点数。

先验概率

初始化先验概率

%Initialize priors
Pr_xy = single(ones(size(x2D)));       %2D
Pr_z = single(ones(size(z1D)));

Pr_xy：初始化 XY 平面上的先验概率，大小与 x2D 相同。
Pr_z：初始化 Z 轴上的先验概率，大小与 z1D 相同。

归一化先验概率

%Normalize Priors
Pr_xy = Pr_xy / sum(Pr_xy(:));       %2D
Pr_z = Pr_z / sum(Pr_z);

将先验概率归一化，使其总和为 1。

计算对数先验概率

%Log Prior
lPr_xy = log(Pr_xy);            %2D
lPr_z = log(Pr_z);

计算先验概率的对数，用于后续的贝叶斯估计。

计算扩散卷积核

%Calculate the diffusion convolution kernels along xy and z seperately
if Dest ~= 0
    %XY diffusion Kernel Generation
    difffull = exp(-coord.^2 / 4 / Dest / tau);
    c = find(difffull >= .1 * max(difffull));
    %check symmetry of coordinates
    mid = ceil(length(difffull) / 2);
    if mid - c(1) == c(end) - mid
        diffkernxy = difffull(c);
    else
        warning('Diffusion kernel is asymmetric along xy')
        diffkernxy = difffull(c);
    end

    %Z diffusion Kernel Generation
    difffull = exp(-z1D.^2 / 4 / Dest / tau);
    c = find(difffull >= .1 * max(difffull));
    %check symmetry of coordinates
    mid = ceil(length(difffull) / 2);
    if mid - c(1) == c(end) - mid
        diffkernz = difffull(c);
    else
        warning('Diffusion kernel is asymmetric along z')
        diffkernz = difffull(c);
    end

    %quality control diffusion coefficients
    if length(diffkernxy) ~= 3 %|| length(diffkernz) ~= 3
        error('diffusion kernel size is larger than 3 pts')
    end

    if abs(diffkernxy(1) - 0.3247) / 0.3247 > 0.5 %|| abs(diffkernz(1) - 0.3247) / 0.3247 > 0.5
        error('Tau modifications have led to deviation from targeted diffusion kernel by more than 50%, actual diffusion kernel value falls outside of range from 0.1623 to 0.4870')
    end
    % disp(['zaxis diffusion kernel is ' num2str(length(diffkernz)) ' units long'])
    clear c mid
end

XY 扩散卷积核生成

difffull：计算整个坐标范围内的扩散卷积核。
c：找到扩散卷积核中值大于最大值 10% 的点。
mid：计算坐标范围的中点。
对称性检查：检查扩散卷积核是否对称。如果对称，则提取对称部分作为扩散卷积核；如果不对称，则发出警告。
diffkernxy：提取对称部分作为 XY 方向的扩散卷积核。

Z 扩散卷积核生成

difffull：计算整个 Z 轴范围内的扩散卷积核。
c：找到扩散卷积核中值大于最大值 10% 的点。
mid：计算坐标范围的中点。
对称性检查：检查扩散卷积核是否对称。如果对称，则提取对称部分作为扩散卷积核；如果不对称，则发出警告。
diffkernz：提取对称部分作为 Z 方向的扩散卷积核。

质量控制

检查扩散卷积核的大小：如果扩散卷积核的大小超过 3 个点，抛出错误。
检查扩散卷积核的值：如果扩散卷积核的第一个值与目标值（0.3247）的偏差超过 50%，抛出错误。目标值范围为 0.1623 到 0.4870。

初始化和生成粒子轨迹

初始化输出变量

% assorted output initializations
stg_out = zeros(N, 3); % 舞台位置输出
photons = zeros(N, 1); % 每步观察到的光子数
posest_out = zeros(N, 3); % 位置估计输出
stageCom = zeros(N, 3); % 舞台命令输出
sb_est = zeros(N, 2); % 信号和背景估计输出
mnb_act = zeros(N, 1); % 最大邻域亮度活动输出

初始化了多个输出变量，这些变量将用于存储模拟过程中生成的数据。

积分时间初始化

% Integration time
errint_x = 0;
errint_y = 0;
errint_z = 0;

初始化了积分时间变量，这些变量可能用于积分误差计算。

定义跟踪开始前的等待步数

% Define steps to wait before turning on tracking
stepstowait = 0;

定义了在开始跟踪之前需要等待的步数。这里设置为 0，表示立即开始跟踪。

生成粒子轨迹

% Generate Trajectory
% Generate steps in x (dx) and y (dy)
dx = sqrt(2 * D * tau) * randn(N, 1)';
dy = sqrt(2 * D * tau) * randn(N, 1)';
dz = sqrt(2 * D * tau) * randn(N, 1)';

使用扩散系数 D 和时间间隔 τ 生成粒子在 x、y 和 z 方向上的随机步长。这些步长遵循高斯分布，其标准差由扩散系数和时间间隔决定。

% Generate XY trajectory
x_part = cumsum(dx)';
y_part = cumsum(dy)';
z_part = cumsum(dz)';
part_out = [x_part, y_part, z_part];

通过累积求和生成粒子在 XY 平面和 Z 轴上的轨迹。
part_out 存储了粒子的完整轨迹。

设置点扩散函数（PSF）大小

% PSF size
sigma_xy = 128e-9; % XY 平面上的 PSF 尺寸
sigma_z = 718e-9; % Z 轴上的 PSF 尺寸
% PSF Measurements empirical Chen Zhang
w = sigma_xy^2 * eye(2, 2); % XY 平面上的 PSF 协方差矩阵
wz = (sigma_z)^2; % Z 轴上的 PSF 方差

定义了 XY 平面和 Z 轴上的点扩散函数（PSF）的大小。
w 是 XY 平面上的 PSF 协方差矩阵，wz 是 Z 轴上的 PSF 方差。

预计算部分伽马值以生成查找表（LUT）

预计算部分伽马值以生成 LUT

% precalculate partial gamma to generate LUT
partialgamma_xy = single(zeros(length(coord), length(coord), 25));
for i = 1:25
    [xlaser, ylaser] = knightsTour(i, d);
    partialgamma_xy(:, :, i) = exp(-((x2D - xlaser).^2 ./ 2 ./ sigma_xy ./ sigma_xy) - ((y2D - ylaser).^2 ./ 2 ./ sigma_xy ./ sigma_xy));
end

partialgamma_xy：预计算部分伽马值，用于生成 XY 平面上的查找表（LUT）。
knightsTour：函数 knightsTour 生成激光在 XY 平面上的位置，这些位置遵循骑士巡逻模式。
xlaser 和 ylaser：激光在 XY 平面上的位置。
partialgamma_xy：计算每个激光位置的高斯分布，用于后续的光子生成。

预计算 Z 轴上的部分伽马值

% partial gamma along z LUT
% calculate p(n) photons for all grid positions assuming laser exists at
% all grid positions
partialgamma_z = zeros(binsperz, length(z1D), 2);
for k = 1:2
    laserPosz = tagLens((1 + (k - 1) * binsperz : k * binsperz), f, tagtau, amp);
    partialgamma_z(:, :, k) = exp(-((z1D - laserPosz').^2 ./ 2 ./ sigma_z ./ sigma_z));
end

partialgamma_z：预计算部分伽马值，用于生成 Z 轴上的查找表（LUT）。
tagLens：函数 tagLens 生成激光在 Z 轴上的位置，这些位置遵循正弦波模式。
laserPosz：激光在 Z 轴上的位置。
partialgamma_z：计算每个激光位置的高斯分布，用于后续的光子生成。

初始化舞台位置

% Stage position
prev_stgxy = [0; 0];
prev_stgz = 0;

prev_stgxy：初始化 XY 平面上的舞台位置。
prev_stgz：初始化 Z 轴上的舞台位置。

初始化光子特征

% Initialize photon characterization for signal and background estimation
pixA = ismember(mod((1:N) - 1, 25) + 1, 23);
pixF = ismember(mod((1:N) - 1, 25) + 1, [1, 3, 5, 7]);
testimate = 15e-3; % rolling average of 15 ms worth of data goes into s and b estimation

binsestimate = int32(testimate / tau); % number of bins that are equal to 15 ms of data
numA = sum(pixA(1:binsestimate));
numF = sum(pixF(1:binsestimate));
pixA = [];
pixF = [];

pixA 和 pixF：定义用于信号和背景估计的光子特征。
testimate：定义用于信号和背景估计的时间窗口（15 毫秒）。
binsestimate：计算时间窗口内的 bin 数。
numA 和 numF：计算在时间窗口内满足特定条件的光子数。

加载权重和脉冲响应数据

% original cal dist weight (σXY = 0.038 µm, σZ = 0.15 µm).
% weight = [0.5905, 1; 1.9853e-7, 1];
weight = [0.44750577, 1; 1.0185348e-05, 1];

% Load impulse response data
measuredImpulseResponse_stg = dlmread('Impulse_Response.txt');
load('galvoimpulseresponse.mat', 'y');
measuredImpulseResponse_glvo = y;

if tau == 20e-6
    % do nothing undecimated Impulse Response Function is Entirely adequate
elseif rem(tau / 20e-6, 1) == 0
    % if integer relationship decimate measured impulse response to
    % appropriate sampling rate
    measuredImpulseResponse_stg = decimate(measuredImpulseResponse_stg, tau / 20e-6);
    measuredImpulseResponse_glvo = decimate(measuredImpulseResponse_glvo, tau / 20e-6);
else
    % interpolate to nearest us sampling then decimate
    measuredImpulseResponse_stg = interp(measuredImpulseResponse_stg, 20);
    measuredImpulseResponse_stg = decimate(measuredImpulseResponse_stg, double(int32(round(tau, 6) / 1e-6)));
    measuredImpulseResponse_glvo = interp(measuredImpulseResponse_glvo, 20);
    measuredImpulseResponse_glvo = decimate(measuredImpulseResponse_glvo, double(int32(round(tau, 6) / 1e-6)));
end

neIR = length(measuredImpulseResponse_stg);
if neIR < 1000
    warning('Fourier Convolution May no Longer Be Temporally Efficient/Faster than traditional conv')
end

weight：定义用于信号和背景估计的权重。
measuredImpulseResponse_stg 和 measuredImpulseResponse_glvo：加载舞台和振镜的脉冲响应数据。
decimate 和 interp：根据时间间隔 τ 调整脉冲响应数据的采样率。
neIR：计算调整后的脉冲响应数据的长度。
warning：如果脉冲响应数据的长度小于 1000，发出警告，提示傅里叶卷积可能不再比传统卷积更高效。

预计算脉冲响应的 FFT

% precalculate fft of measured Impulse Response
Ly = length(measuredImpulseResponse_stg) * 2 - 1;  % Length of the impulse response
ftzerpad = pow2(nextpow2(Ly));    % Find smallest power of 2 that is > Ly
fftExpImpResp = fft(measuredImpulseResponse_stg, ftzerpad); % Fast Fourier transform
lengthImpResp = length(measuredImpulseResponse_stg);

Ly：脉冲响应的长度。
ftzerpad：找到大于 Ly 的最小 2 的幂，用于零填充。
fftExpImpResp：计算脉冲响应的快速傅里叶变换（FFT）。
lengthImpResp：脉冲响应的实际长度。

粒子跟踪模拟

光子的生成

初始化循环变量

% Now simulate tracking
for k = int32(1:N)

k：当前时间步长，从 1 到 N。

获取激光位置

% First, get simulated observed number of photons based on the particle
% and laser positions
[xlaser, ylaser] = knightsTour(k, d);   % xy laser positions
laserPosxy = [xlaser; ylaser];

knightsTour：函数 knightsTour 生成激光在 XY 平面上的位置，这些位置遵循骑士巡逻模式。
xlaser 和 ylaser：激光在 XY 平面上的坐标。
laserPosxy：将激光的 XY 坐标组合成一个 2x1 的向量。

获取 Z 轴上的激光位置

% pull tag phase for subbins
laserPosz = tagLens((1 + (k - 1) * binsperz : k * binsperz), f, tagtau, amp);

tagLens：函数 tagLens 生成激光在 Z 轴上的位置，这些位置遵循正弦波模式。
laserPosz：激光在 Z 轴上的位置，每个时间步长有 binsperz 个子 bin。

检查 Z 轴上的分辨率

if k == 1 % on first loop check relative spacing.
    edges = linspace(-1e-6, 1e-6, 3); % only need two bins high and low n, less bins less math.
    if mean(diff(unique(laserPosz))) >= mean(diff(edges))
        warning('photons are being generated at a lower z resolution than they are being assigned to positions for bayesian estimation')
    end
    [nbinperk, ~, bin] = histcounts(laserPosz, edges);
    bin = bin'; % indexes must be column for accumarray.
    nbinperk = nbinperk'; % pretty sure I need this to be a column also for consistency
    laserPosz_binned = accumarray(bin, laserPosz) ./ nbinperk;
    partialgamma_z = exp(-(z1D - laserPosz_binned).^2 ./ 2 ./ sigma_z / sigma_z);
end

edges：定义 Z 轴上的两个边界。
nbinperk：每个 bin 中的光子数。
bin：每个光子所属的 bin。
laserPosz_binned：将激光位置分到不同的 bin 中。
partialgamma_z：计算每个 bin 的高斯分布。

计算激光位置到初始舞台位置的距离

% radius of laser position from initial stage position 0,0 to determine
% background intensity
laser_r = sqrt(sum((prev_stgxy + laserPosxy).^2) + laserPosz.^2);

laser_r：激光位置到初始舞台位置的距离。

获取粒子位置

% particlePosxy for given step in correct matrix format
particlePosxy = [x_part(k); y_part(k)];

particlePosxy：粒子在 XY 平面上的位置。

信号和背景的估计

确定背景强度

% Determine actual bg at given laser positions
for zscn = 1:binsperz
    if laser_r(zscn) > r(end)
        % if laser radius is outside last point of bg increase, set to value
        % at furthest radius
        b(zscn) = bofr(end);
    else
        % otherwise set background to appropriate level
        b(zscn) = bofr(find(laser_r(zscn) < r, 1, 'first'));
    end
end
mnb_act(k) = mean(b);

b：根据激光位置计算背景强度。
mnb_act：存储每个时间步长的平均背景强度。

计算观察到的光子数

% calculate observed number of photons
nseq = poissrnd(tagtau .* emRate((particlePosxy - prev_stgxy), (z_part(k) - prev_stgz), laserPosxy, laserPosz, b, s(k), w, wz)); % number of photons
nseq = nseq';
n = sum(nseq);
photons(k) = n; % store total photons for output
ktindx = mod(k - 1, 25) + 1;

nseq：根据粒子和激光位置生成观察到的光子数。
photons：存储每个时间步长的总光子数。
ktindx：当前时间步长在骑士巡逻模式中的索引。

将光子数添加到特定像素

% if photons belong to pixels of interest append them to photon string
if ismember(ktindx, 23)
    pixA = [pixA, photons(k)];
elseif ismember(ktindx, [1, 3, 5, 7])
    pixF = [pixF, photons(k)];
end

pixA 和 pixF：存储特定像素的光子数。

估计信号和背景水平

% estimate signal and bg levels based on observed photons
if k <= binsestimate
    % some given estimate. Can be accurate or could be not. We will make
    % it accurate. But at less than 15 ms of data that's hard to deal
    % with
    [calc_sb] = [s(k); mnb_act(k)];
else
    % above 15 ms of available data estimate using only the most recent
    % 15 ms worth of data.
    % calculate index value for end of given photon type string
    % then select most recent desired number of photons for a given sb
    % estimation time window calculated early
    A = length(pixA);
    photonsA = pixA(A - numA + 1:A);
    
    F = length(pixF);
    photonsF = pixF(F - numF + 1:F);
    
    [calc_sb] = weight \ [sum(photonsA) / (length(photonsA) * tau); sum(photonsF) / (length(photonsF) * tau)];
end
if p.Results.sbest == false
    % give perfect info
    [calc_sb] = [s(k); mnb_act(k)];
end

calc_sb：估计信号和背景水平。calc_sb(1)：表示信号（signal）的估计值；calc_sb(2)：表示背景（background）的估计值。
pixA 和 pixF：用于信号和背景估计的光子数。
weight：权重矩阵，用于加权最小二乘估计。

检查信号和背景估计

% if signal estimate is negative or 0 crash
if calc_sb(1) <= 0
    warning('negative or 0 estimated signal')
    break
end

% if bg is negative replace it with 0
if calc_sb(2) < 0
    calc_sb(2) = 0;
end

% store signal/background estimates for a given step along with the
% actual bg level
sb_est(k, :) = calc_sb';

calc_sb：确保信号和背景估计值为正。
sb_est：存储每个时间步长的信号和背景估计值。

使用贝叶斯方法估计粒子位置

2D(XY)位置

% use 2D Bayesian to Estimate Particle Position Along X and Y
gammaxy = calc_sb(1) * partialgamma_xy(:, :, ktindx) + calc_sb(2);
log_n_factorial = stirling(n);
lP_n_xy = n * log(gammaxy * tau) - gammaxy * tau - log_n_factorial;
lPoxy = lPr_xy + lP_n_xy;
lPoxy = lPoxy - max(max(lPoxy));
lPoxy = lPoxy - log(sum(sum(exp(lPoxy))));
if Dest ~= 0
    lPr_xy = log(conv2(diffkernxy, diffkernxy, exp(lPoxy), 'same'));
    lPr_xy(isinf(lPr_xy)) = -745;
else
    lPr_xy = lPoxy;
end