fbLSC-library/g2Fit.m at main · BUNPC/fbLSC-library · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
function [tau_c,beta] = g2Fit(dt,g2,varargin)
%g2Fit Fits g2 curve to get that tau value
%   Inputs:
%       dt = time difference
%       g2 = g2 curve
%   Outputs:
%       tau_c
%       beta

p = inputParser;
addRequired(p,'dt');
addRequired(p,'g2');
addParameter(p,'Beta',0,@isnumeric);
addParameter(p,'g2Thr',0.7,@isnumeric);
parse(p,dt,g2,varargin{:});
beta0 = p.Results.Beta;
g2thr_ratio = p.Results.g2Thr;

if size(dt,2) ~= 1
    dt = dt';
end

% get the dimensions right
if size(g2,1) ~= numel(dt)
    g2 = g2';
end

%
if beta0 == 0
    fixed_beta = false;
else
    fixed_beta = true;
end

good_tau = dt >= 1e-8; % anything less than 10 ns probably has artifacts with current SPADs
dt = dt(good_tau);
g2 = g2(good_tau,:);
config_num = size(g2,2);

tau_c = nan(config_num,1);
beta = nan(config_num,1);

if fixed_beta
    fun = @(x,tau)(beta0*exp(-tau*x)+1);
    options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
    options.Display = 'off';
    lb = 0;
    ub = [];
else
    fun = @(x,tau)(x(1)*exp(-tau*x(2))+1);
    options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
    options.Display = 'off';
    lb = [0 0];
    ub = [];
end

options.FunctionTolerance = 1e-8;

% get reference beta even if no fixed beta
if ~fixed_beta
    x0 = [max(mean(g2,2)) - 1,0.001];
    [x,~,~,~,~]=lsqcurvefit(fun,x0,dt,mean(g2,2),lb,ub,options);
    beta0 = x(1);

    % only select the g2 that did not decay much
    thr = fun(x,dt)-1 >= beta0*g2thr_ratio;
    thr = find(thr); thr = thr(end);
else
    x0 = 0.001;
    [x,~,~,~,~]=lsqcurvefit(fun,x0,dt,mean(g2,2),lb,ub,options);

    thr = fun(x,dt)-1 >= beta0*g2thr_ratio;
    thr = find(thr); thr = thr(end);
end


for config=1:config_num
    g2_config = g2(:,config);

    if fixed_beta
        x0 = 0.001;
    else
        x0 = [max(g2_config) - 1,0.001];
    end

    g2_config = g2_config(1:thr);
    [x,~,~,exitflag(config),~]=lsqcurvefit(fun,x0,dt(1:thr),g2_config,lb,ub,options);

    if fixed_beta
        tau_c(config)=1./x;
        beta(config)=beta0;
    else
        tau_c(config)=1./x(2);
        beta(config)=x(1);
    end

%     y = log(g2_config - 1);
%     good_ind = imag(y) < 100*eps;
%     y = y(good_ind);
%     x = log(dt);

%     fun = @(x,tau)(x(1)*exp(-exp(tau)*x(2))+1);
%     [b,~,~,exitflag(config),~]=lsqcurvefit(fun,x0,x,g2_config,lb,ub,options);
%     tau(config)=1./b(2);
%     beta(config)=b(1);
%     [b,bint] = regress(y,[ones(size(y)) dt(good_ind)]);
end
end