-
Notifications
You must be signed in to change notification settings - Fork 1
/
monochromaticLightAsGaussian.m
100 lines (73 loc) · 3.96 KB
/
monochromaticLightAsGaussian.m
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
function SPD = monochromaticLightAsGaussian(peak, FWHM, lambda)
% to demo how to reconstruct the used light in some experiment based on
% the peak wavelength and the half-bandwidth (HBW or FWHM)
% Petteri Teikari, INSERM, 2010, petteri.teikari@inserm.fr
%% user-defined parameters
% peak = 480; % nm
% FWHM = 10; % +/- nm
% sigma / deviation of the gaussian is defined as a function of FWHM
sigma = (FWHM*2) / 2.355; % e.g., http://cp.literature.agilent.com/litweb/pdf/5980-0746E.pdf
%% create a GAUSSIAN Spectral Power Distribution (SPD) with the above
%% parameters and lambda
f = gauss_distribution(lambda, peak, sigma);
SPD = f / max(f); % normalize
% Get the HBW / FWHM of the created light source, just for debug
% hbw = getHBW(lambda, f);
function f = gauss_distribution(x, mu, s)
%% subfunction for Gaussian distribution
% x x vector (i.e. wavelength)
% mu mean / peak nm
% sigma standard deviation
p1_1 = ((x - mu)/s) .^ 2;
p1 = -.5 * p1_1;
p2 = (s * sqrt(2*pi));
f = exp(p1) ./ p2;
function hbw = getHBW(x, y)
%% Subfunction to calculate the HBW of the created light
% output
% hbw structure containing various parameters
% inputs
% x wavelength of the light SPD
% y spectral irradiance/radiance/photon flux/etc. of the light
hbw_safeLimit = 0.03; % safe limit to find the peak
hbw_upsampleFactor = 10; % upsample factor for better spectral res.
hbw_interpMode = 'pchip'; % interpolation method for upsampling
% get the half-bandwidth of the spectrum (assuming that it
% is a monochromatic spectrum with a Gaussian form)
% normalize the spectrum
y = y / max(y);
% get a small portion around the peak for further interpolation
ind_min = find(y > hbw_safeLimit, 1, 'first');
ind_max = find(y > hbw_safeLimit, 1, 'last');
x_hbw = x(ind_min:ind_max);
y_hbw = y(ind_min:ind_max);
vectorSizeReduction = (length(x_hbw) / length(x)); % debug variable
% interpolate the vectors for better spatial resolution
noOfPoints = (length(x_hbw) * hbw_upsampleFactor);
x_interp = linspace(min(x_hbw), max(x_hbw), noOfPoints)';
pp = interp1(x_hbw, y_hbw, hbw_interpMode, 'pp');
y_interp = ppval(pp, x_interp);
% get the difference between the interpolated y and 0.5
% (brute force method) to find the points closest to 0.5
% (in other words for the hbw points), it is very likely
% that the y vector does not contain values EXACTLY 0.5000
y_diff = abs(y_interp - 0.5);
% sort the rows in the hope that first two values would
% correspond to the HBW points, this could be improved at
% some point adding more constraints if it seems that it
% starts to fail
[y_diff,index] = sortrows(y_diff);
index_1 = index(1); index_2 = index(2);
wavel_1 = x_interp(index_1); wavel_2 = x_interp(index_2);
hbw.min = min([wavel_1 wavel_2]);
hbw.max = max([wavel_1 wavel_2]);
[maxValue maxInd] = max(y_interp);
hbw.peak = x_interp(maxInd);
hbw.range = hbw.max - hbw.min;
hbw.interp_x = x_interp;
hbw.interp_y = y_interp;
hbw.wavel_1 = wavel_1;
hbw.wavel_2 = wavel_2;
hbw.index1 = index_1;
hbw.index2 = index_2;
hbw.peakRounded = round(hbw.peak);