Skip to content

Commit 40ce1aa

Browse files
JoonasJLJoonasJL
committedFeb 19, 2025·
KF functions for the class implementation
1 parent 7e6dbc1 commit 40ce1aa

16 files changed

+635
-0
lines changed
 

‎+plugins/+ClassKF/EnKF.m

+61
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
function [z_inverse] = EnKF(m, A, P, Q, L, R, timeSteps, number_of_frames, n_ensembles, q_given)
2+
%ENKF Summary of this function goes here
3+
%x_ensemble = mvnrnd(zeros(size(m)), 100* ones(size(m,1)), n_ensembles)';
4+
x_ensemble = mvnrnd(zeros(size(m)), P, n_ensembles)';
5+
z_inverse = cell(0);
6+
if not(q_given)
7+
q_values = Q;
8+
end
9+
h = zef_waitbar(0, 'EnKF Filtering');
10+
for f_ind = 1:number_of_frames
11+
zef_waitbar(f_ind/number_of_frames,h,...
12+
['EnKF Filtering ' int2str(f_ind) ' of ' int2str(number_of_frames) '.']);
13+
f = timeSteps{f_ind};
14+
if not(q_given)
15+
Q = diag(q_values(:,f_ind));
16+
end
17+
w = mvnrnd(zeros(size(m)), Q, n_ensembles)';
18+
% Forecasts
19+
20+
x_f = A * x_ensemble + w;
21+
C = cov(x_f');
22+
correlationLocalization = true;
23+
if correlationLocalization
24+
T = corrcoef(x_f');
25+
% explain How to find 0.05
26+
T(abs(T) < 0.05) = 0;
27+
C = C .* T;
28+
end
29+
v = mvnrnd(zeros(size(R,1),1), R, n_ensembles);
30+
31+
% method to calculate resolution D
32+
method = '3';
33+
if(method == '1')
34+
P_sqrtm = sqrtm(C);
35+
B = L * P_sqrtm;
36+
G = B' / (B * B' + R);
37+
w_t = 1 ./ sum(G.' .* B, 1)';
38+
D = w_t .* inv(P_sqrtm);
39+
elseif(method == '2')
40+
% complexity O(n^3)
41+
[Ur,Sr,Vr] = svd(C);
42+
Sr = diag(Sr);
43+
RNK = sum(Sr > (length(Sr) * eps(single(Sr(1)))));
44+
SIR = Vr(:,1:RNK) * diag(1./sqrt(Sr(1:RNK))) * Ur(:,1:RNK)'; % square root
45+
P_sqrtm = Vr(:,1:RNK) * diag(sqrt(Sr(1:RNK))) * Ur(:,1:RNK)';
46+
B = L * P_sqrtm;
47+
G = B' / (B * B' + R);
48+
w_t = 1 ./ sum(G.' .* B, 1)';
49+
D = w_t .* SIR;
50+
else
51+
D = speye(size(C));
52+
end
53+
% Update
54+
K = C * L' / (L * C * L' + R);
55+
x_ensemble = x_f + K *(f + v' - L*x_f);
56+
% x_ensemble = x_ensemble';
57+
mean_x = mean(x_ensemble,2);
58+
z_inverse{f_ind} = D*mean_x;
59+
end
60+
close(h);
61+
end

‎+plugins/+ClassKF/Q_quantities.m

+35
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
function [sigma, phi, B, C, D] = Q_quantities(P, m, G, y)
2+
3+
T = size(m,2) - 1;
4+
5+
sigma = zeros(size(P{1}));
6+
for k = 1:T
7+
sigma = sigma + P{k+1} + m{k+1}*m{k+1}';
8+
end
9+
sigma = sigma* 1/T;
10+
11+
phi = zeros(size(P{1}));
12+
for k = 1:T
13+
phi = phi + P{k} + m{k}*m{k}';
14+
end
15+
phi = 1/T * phi;
16+
17+
B = zeros(size(y{1},1), size(m{1}',2));
18+
for k = 1:T
19+
B = B + y{k+1} * m{k+1}';
20+
end
21+
B = 1/T * B;
22+
23+
C = zeros(size(P{1},1),size( G{1}', 2));
24+
for k = 1:T
25+
C = P{k+1} * G{k}' + m{k+1} * m{k}';
26+
end
27+
C = 1/T * C;
28+
29+
D = zeros(size(y,1));
30+
for k = 1:T
31+
D = y{k+1} * y{k+1}';
32+
end
33+
D = D * 1/T;
34+
35+
end

‎+plugins/+ClassKF/RTS_smoother.m

+34
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
function [P_s_store, m_s_store, G_store] = RTS_smoother(P_store, z_inverse, A, Q, number_of_frames)
2+
P_s_store = cell(0);
3+
m_s_store = cell(0);
4+
G_store = cell(0);
5+
h = zef_waitbar(0,'Smoothing');
6+
for f_ind = number_of_frames:-1:1
7+
zef_waitbar(1 - f_ind/number_of_frames,h, ['Smoothing ' int2str(number_of_frames -f_ind) ' of ' int2str(number_of_frames) '.']);
8+
9+
P = P_store{f_ind};
10+
m = z_inverse{f_ind};
11+
% if A is Identity
12+
if (isdiag(A) && all(diag(A) - 1) < eps)
13+
P_ = P + Q;
14+
m_ = m;
15+
G = P / P_;
16+
else
17+
P_ = A * P * A' + Q;
18+
m_ = A * m;
19+
G = (P * A') / P_;
20+
end
21+
if f_ind == number_of_frames
22+
m_s = m;
23+
P_s = P;
24+
else
25+
m_s = m + G * (m_s - m_);
26+
P_s = P + G * (P_s - P_) * G';
27+
end
28+
P_s_store{f_ind} = P_s;
29+
G_store{f_ind} = G;
30+
m_s_store{f_ind} = m_s;
31+
end
32+
33+
close(h);
34+
end

‎+plugins/+ClassKF/class_kf_predict.m

+11
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
function [m, P] = class_kf_predict(KFclassObj)
2+
% kf_predict is the prediction step of kalman filter
3+
% Skip multiplications if A is Indentity
4+
if (isdiag(KFclassObj.state_transition_model_A) && all(diag(KFclassObj.state_transition_model_A) - 1) < eps)
5+
P = KFclassObj.prev_step_posterior_cov + KFclassObj.evolution_cov;
6+
else
7+
% Basic kalman prediction steps
8+
m = KFclassObj.state_transition_model_A * self.prev_step_reconstruction;
9+
P = KFclassObj.state_transition_model_A * KFclassObj.prev_step_posterior_cov * KFclassObj.state_transition_model_A' + KFclassObj.evolution_cov;
10+
end
11+
end

‎+plugins/+ClassKF/connectivity_matrix.m

+28
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
function [A] = connectivity_matrix(source_positions, K, weighted_avg)
2+
% Connectivity matrix
3+
%SMOOTHING_MATRIX Summary of this function goes here
4+
% Detailed explanation goes here
5+
% TODO: weighted avg
6+
if nargin < 3
7+
weighted_avg = false;
8+
end
9+
10+
MdlKDT = KDTreeSearcher(source_positions);
11+
[IDX, D] = knnsearch(MdlKDT, source_positions, 'K', K);
12+
A = sparse(size(source_positions,1), size(source_positions,1));
13+
14+
row_idx = 1:size(A,1);
15+
16+
for n_neighbor = 1:K
17+
n_neighbor_ind = IDX(:,n_neighbor);
18+
idx = sub2ind(size(A),row_idx' , n_neighbor_ind);
19+
A(idx) = 1/K;
20+
% if use weighted avg. D has distances, so it needs adaptation
21+
%A(idx) = D(:,n_neighbor);
22+
end
23+
24+
% z is indexed as x_1,y_1,z_1,x_2,y_2,z_2,... so on.
25+
% A has the factor for each point, but factors need to expand to each x,y
26+
% and z component
27+
A = kron(A,eye(3));
28+
end

‎+plugins/+ClassKF/find_evolution_prior.m

+37
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
function [q] = find_evolution_prior(L, f_data, likelihood_std, evolution_prior_db,evolution_mode)
2+
% q to evolution prior
3+
4+
switch evolution_mode
5+
case 1
6+
%Since R = likelihood_std^2*eye, we have
7+
% A_noise = (1/(p*A_signal))*A_signal => SNR = p^2*A_signal^2 =
8+
% (q||L||^2)^2/E[||noise||^2]
9+
% => E[dy^2] -> sqrt(E[dy^2])
10+
f = sqrt(mean(diff(f_data').^2,2))*10^(evolution_prior_db/20);
11+
f = [f;f(end)];
12+
q = transpose((1-likelihood_std^2)*f./repelem(sum(reshape(sum(L.^2),3,[])),3));
13+
case 2
14+
%case 1 but spatial sensitivity is averaged
15+
f = sqrt(mean(diff(f_data').^2,2))*10^(evolution_prior_db/20);
16+
f = [f;f(end)];
17+
q = transpose((1-likelihood_std^2)*f/mean(repelem(sum(reshape(sum(L.^2),3,[])),3)));
18+
case 3
19+
%Mathematically quarantees a good tracking but numerically instable
20+
%in 2024
21+
[~,S,V] = svd(L,"econ");
22+
f = diff(f_data')';
23+
f = 10^(evolution_prior_db/20)*sum(f.^2,2)./sum(f_data.^2,2);
24+
S = max((diag(S).^2),likelihood_std^2/(1-likelihood_std^2));
25+
q = (V.*(f./S)')*V';
26+
case 4
27+
%averaged signal space contribution
28+
S = svd(L);
29+
f = diff(f_data')';
30+
f = sum(f.^2,2)./sum(f_data.^2,2);
31+
S = max(S.^2,likelihood_std^2/(1-likelihood_std^2));
32+
q = mean(f./S)*10^(evolution_prior_db/20);
33+
case 5
34+
q = time_step*(svds(L,1).^(2)/sum(L(:).^2))*10^(evolution_prior_db/20);
35+
end
36+
end
37+

‎+plugins/+ClassKF/kalman_filter.m

+27
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
function [P_store,z_inverse] = kalman_filter(m,P,A,Q,L,R, timeSteps ,number_of_frames, smoothing, q_given)
2+
P_store = cell(0);
3+
z_inverse = cell(0);
4+
if not(q_given)
5+
q_values = Q;
6+
end
7+
h = zef_waitbar(0, 'Filtering');
8+
for f_ind = 1: number_of_frames
9+
zef_waitbar(f_ind/number_of_frames,h,...
10+
['Filtering ' int2str(f_ind) ' of ' int2str(number_of_frames) '.']);
11+
f = timeSteps{f_ind};
12+
if not(q_given)
13+
Q = diag(q_values(:,f_ind));
14+
end
15+
% Prediction
16+
%[m,P] = zeffiro.plugins.Kalman.m.kf_predict(m, P, A, Q);
17+
[m,P] = kf_predict(m, P, A, Q);
18+
% Update
19+
%[m, P] = zeffiro.plugins.Kalman.m.kf_update(m, P, f, L, R);
20+
[m, P] = kf_update(m, P, f, L, R);
21+
if (smoothing == 2)
22+
P_store{f_ind} = gather(P);
23+
end
24+
z_inverse{f_ind} = gather(m);
25+
end
26+
close(h);
27+
end

‎+plugins/+ClassKF/kalman_filter_sLORETA.m

+26
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
function [P_store,z_inverse] = kalman_filter_sLORETA(m,P,A,Q,L,R, timeSteps ,number_of_frames, smoothing, q_given)
2+
P_store = cell(0);
3+
z_inverse = cell(0);
4+
h = zef_waitbar(0, 'Filtering');
5+
if not(q_given)
6+
q_values = Q;
7+
end
8+
9+
for f_ind = 1: number_of_frames
10+
zef_waitbar(f_ind/number_of_frames,h,...
11+
['Filtering ' int2str(f_ind) ' of ' int2str(number_of_frames) '.']);
12+
f = timeSteps{f_ind};
13+
if not(q_given)
14+
Q = diag(q_values(:,f_ind));
15+
end
16+
% Prediction
17+
[m,P] = kf_predict(m, P, A, Q);
18+
% Update
19+
[m, P, ~, D] = kf_sL_update(m, P, f, L, R);
20+
z_inverse{f_ind} = gather(D*m);
21+
if (smoothing == 2)
22+
P_store{f_ind} = gather(P);
23+
end
24+
end
25+
close(h);
26+
end

‎+plugins/+ClassKF/kf_predict.m

+11
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
function [m, P] = kf_predict(m,P,A,Q)
2+
% kf_predict is the prediction step of kalman filter
3+
% Skip multiplications if A is Indentity
4+
if (isdiag(A) && all(diag(A) - 1) < eps)
5+
P = P + Q;
6+
else
7+
% Basic kalman prediction steps
8+
m = A * m;
9+
P = A * P * A' + Q;
10+
end
11+
end

‎+plugins/+ClassKF/kf_sL_update.m

+32
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
function [m, P, K, D] = kf_sL_update(m,P,y,H,R)
2+
% Resolution matrix
3+
method = '1';
4+
if(method == '1')
5+
P_sqrtm = sqrtm(P);
6+
B = H * P_sqrtm;
7+
G = B' / (B * B' + R);
8+
w_t = 1 ./ sum(G.' .* B, 1)';
9+
D = w_t .* inv(P_sqrtm);
10+
elseif(method == '2')
11+
[Ur,Sr,Vr] = svd(P);
12+
Sr = diag(Sr);
13+
RNK = sum(Sr > (length(Sr) * eps(single(Sr(1)))));
14+
SIR = Vr(:,1:RNK) * diag(1./sqrt(Sr(1:RNK))) * Ur(:,1:RNK)'; % square root
15+
P_sqrtm = Vr(:,1:RNK) * diag(sqrt(Sr(1:RNK))) * Ur(:,1:RNK)';
16+
B = H * P_sqrtm;
17+
G = B' / (B * B' + R);
18+
w_t = 1 ./ sum(G.' .* B, 1)';
19+
D = w_t .* SIR;
20+
end
21+
% kf_update is the update step of kalman filter
22+
v = y-H*m;
23+
24+
S = H*P*H'+R;
25+
K = (P*H')/S; % /S is same as *inv(S) but faster and more accurate
26+
27+
m = m+K*v;
28+
P = P-K*S*K';
29+
30+
31+
end
32+

‎+plugins/+ClassKF/kf_sL_update_approx.m

+33
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
function [m, P, K, D] = kf_sL_update_approx(m,P,y,H,R)
2+
% Approximating the inverse of square root matrix
3+
N = 5; M = 1;
4+
Z = eye(length(m));
5+
Y = P;
6+
invY = Z;
7+
invZ = Z;
8+
for n = 1:N
9+
for k = 1:M
10+
invY = 2*invY-invY*Y*invY;
11+
invZ = 2*invZ-invZ*Z*invZ;
12+
end
13+
Y = 0.5*(Y+invZ);
14+
Z = 0.5*(Z+invY);
15+
end
16+
P_sqrtm_right = Z;
17+
B = H * P;
18+
K = B*P_sqrtm_right;
19+
G = K' / (B * H' + R);
20+
w_t = 1 ./ sum(G.' .* K, 1)';
21+
D = w_t .* P_sqrtm_right;
22+
% kf_update is the update step of kalman filter
23+
v = y-H*m;
24+
25+
S = H*P*H'+R;
26+
K = (P*H')/S; % /S is same as *inv(S) but faster and more accurate
27+
28+
m = m+K*v;
29+
P = P-K*S*K';
30+
31+
32+
end
33+

‎+plugins/+ClassKF/kf_update.m

+11
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
function [m, P, K] = kf_update(m,P,y,H,R)
2+
% kf_update is the update step of kalman filter
3+
v = y-H*m;
4+
% downsample? H, P, R before S
5+
S = H*P*H'+R;
6+
K = (P*H')/S; % /S is same as *inv(S) but faster and more accurate
7+
% upsample? K, P
8+
m = m+K*v;
9+
P = P-K*S*K';
10+
% upsample P
11+
end

‎+plugins/+ClassKF/zef_KF.m

+122
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
function [zef] = zef_KF(zef, q_value)
2+
% Optimal q_value as parameter
3+
%% Initial parameters
4+
snr_val = zef.inv_snr;
5+
pm_val = zef.inv_prior_over_measurement_db;
6+
amplitude_db = zef.inv_amplitude_db;
7+
pm_val = pm_val - amplitude_db;
8+
std_lhood = 10^(-snr_val/20);
9+
sampling_freq = zef.inv_sampling_frequency;
10+
high_pass = zef.inv_low_cut_frequency;
11+
low_pass = zef.inv_high_cut_frequency;
12+
number_of_frames = zef.number_of_frames;
13+
source_direction_mode = zef.source_direction_mode;
14+
source_directions = zef.source_directions;
15+
source_positions = zef.source_positions;
16+
time_step = zef.inv_time_3;
17+
18+
%% Reconstruction identifiers
19+
reconstruction_information.tag = 'Kalman';
20+
reconstruction_information.inv_time_1 = zef.inv_time_1;
21+
reconstruction_information.inv_time_2 = zef.inv_time_2;
22+
reconstruction_information.inv_time_3 = zef.inv_time_3;
23+
reconstruction_information.sampling_freq = zef.inv_sampling_frequency;
24+
reconstruction_information.low_pass = zef.inv_high_cut_frequency;
25+
reconstruction_information.high_pass = zef.inv_low_cut_frequency;
26+
reconstruction_information.number_of_frames = zef.number_of_frames;
27+
reconstruction_information.source_direction_mode = zef.source_direction_mode;
28+
reconstruction_information.source_directions = zef.source_directions;
29+
reconstruction_information.snr_val = zef.inv_snr;
30+
reconstruction_information.pm_val = zef.inv_prior_over_measurement_db;
31+
32+
%%
33+
[L,n_interp, procFile] = zef_processLeadfields(zef);
34+
35+
%get ellipse filteres full measurement data. f_data: "sensors" x "time points"
36+
[f_data] = zef_getFilteredData(zef);
37+
timeSteps = arrayfun(@(x) zef_getTimeStep(f_data, x, zef), 1:number_of_frames, 'UniformOutput', false);
38+
39+
z_inverse_results = cell(0);
40+
%% CALCULATION STARTS HERE
41+
% m_0 = prior mean
42+
m = zeros(size(L,2), 1);
43+
44+
%[theta0] = zef_find_gaussian_prior(snr_val-pm_val,L,size(L,2),zef.normalize_data,0);
45+
46+
% Transition matrix is Identity matrix
47+
%P = eye(size(L,2)) * theta0;
48+
49+
P = (1-std_lhood^2)*10.^(-(snr_val-pm_val)/10)*diag(mean(var(f_data(:,1:zef.kf_number_of_noise_steps),0,2))./repelem(sum(reshape(sum(L.^2),3,[])),3));
50+
51+
A = eye(size(L,2));
52+
q_given_flag = false;
53+
% If q_value given in the function call
54+
if nargin > 1
55+
if max(size(q_value)) == 1
56+
Q = q_value*eye(size(L,2));
57+
elseif min(size(q_value)) == 1
58+
Q = diiag(q_value(:));
59+
else
60+
Q = q_value;
61+
end
62+
q_given_flag = true;
63+
else
64+
zef_init_gaussian_prior_options;
65+
evolution_prior_db = zef.inv_evolution_prior;
66+
q_value = find_evolution_prior(L, f_data, std_lhood, evolution_prior_db, zef.kf_evolution_prior_mode);
67+
if zef.kf_evolution_prior_mode <= 2
68+
%time-dependent models
69+
Q = q_value;
70+
else
71+
q_given_flag = true;
72+
if zef.kf_evolution_prior_mode == 3
73+
%full matrix model
74+
Q = q_value;
75+
else
76+
%iid models
77+
Q = q_value*eye(size(L,2));
78+
end
79+
end
80+
end
81+
reconstruction_information.Q = q_value;
82+
clear q_value
83+
% std_lhood
84+
R = std_lhood^2 * eye(size(L,1));
85+
86+
Q_Store = cell(0);
87+
88+
%% KALMAN FILTER
89+
filter_type = zef.filter_type;
90+
smoothing = zef.kf_smoothing;
91+
if filter_type == 1
92+
[P_store, z_inverse] = kalman_filter(m,P,A,Q,L,R,timeSteps, number_of_frames, smoothing,q_given_flag);
93+
elseif filter_type == 2
94+
n_ensembles = zef.number_of_ensembles;
95+
z_inverse = EnKF(m,A,P,Q,L,R,timeSteps,number_of_frames, n_ensembles,q_given_flag);
96+
elseif filter_type == 3
97+
[P_store, z_inverse] = kalman_filter_sLORETA(m,P,A,Q,L,R,timeSteps, number_of_frames, smoothing,q_given_flag);
98+
elseif filter_type == 4
99+
[P_store, z_inverse] = kalman_filter(m,P,A,Q,L,R,timeSteps, number_of_frames, smoothing,q_given_flag);
100+
P_old = eye(size(L,2)) * theta0;
101+
H = L * sqrtm(P_old);
102+
W = inv(sqrtm(diag(diag(H'*inv(H*H' + R)*H))));
103+
z_inverse = cellfun(@(x) W*x, z_inverse, 'UniformOutput', false);
104+
105+
end
106+
107+
%% RTS SMOOTHING
108+
109+
if (smoothing == 2)
110+
[~, m_s_store, ~] = RTS_smoother(P_store, z_inverse, A, Q, number_of_frames);
111+
z_inverse = m_s_store;
112+
end
113+
114+
%% POSTPROCESSING
115+
[z] = zef_postProcessInverse(z_inverse, procFile);
116+
%normalize the reconstruction so that the highest value is equal to 1
117+
[z] = zef_normalizeInverseReconstruction(z);
118+
%% CALCULATION ENDS HERE
119+
zef.reconstruction_information = reconstruction_information;
120+
zef.reconstruction = z;
121+
122+
end

‎+plugins/+ClassKF/zef_kf_open_window.m

+108
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,108 @@
1+
%This is the startup function for KF app.
2+
%Beamformer, inverse_tools, zef_KF_start
3+
function zef = zef_kf_open_window(zef)
4+
5+
if nargin == 0
6+
zef = evalin('base','zef');
7+
end
8+
9+
zef.KF = zef_kf_app;
10+
11+
%_ Initial values _
12+
if isfield(zef,'inv_snr')
13+
zef.KF.inv_snr.Value = num2str(zef.inv_snr);
14+
else
15+
zef.KF.inv_snr.Value = '30';
16+
end
17+
18+
if isfield(zef,'inv_sampling_frequency')
19+
zef.KF.inv_sampling_frequency.Value = num2str(zef.inv_sampling_frequency);
20+
else
21+
zef.KF.inv_sampling_frequency.Value = '1025';
22+
end
23+
24+
if isfield(zef,'inv_low_cut_frequency')
25+
zef.KF.inv_low_cut_frequency.Value = num2str(zef.inv_low_cut_frequency);
26+
else
27+
zef.KF.inv_low_cut_frequency.Value = '7';
28+
end
29+
30+
if isfield(zef,'inv_high_cut_frequency')
31+
zef.KF.inv_high_cut_frequency.Value = num2str(zef.inv_high_cut_frequency);
32+
else
33+
zef.KF.inv_high_cut_frequency.Value = '9';
34+
end
35+
36+
if isfield(zef,'inv_time_1')
37+
zef.KF.inv_time_1.Value = num2str(zef.inv_time_1);
38+
else
39+
zef.KF.inv_time_1.Value = '0';
40+
end
41+
42+
if isfield(zef,'inv_time_2')
43+
zef.KF.inv_time_2.Value = num2str(zef.inv_time_2);
44+
else
45+
zef.KF.inv_time_2.Value = '0';
46+
end
47+
48+
if isfield(zef,'number_of_frames')
49+
zef.KF.number_of_frames.Value = num2str(zef.number_of_frames);
50+
else
51+
zef.KF.number_of_frames.Value = '0';
52+
end
53+
54+
if isfield(zef,'inv_time_3')
55+
zef.KF.inv_time_3.Value = num2str(zef.inv_time_3);
56+
else
57+
zef.KF.inv_time_3.Value = '0';
58+
end
59+
60+
if isfield(zef,'filter_type')
61+
zef.KF.filter_type.Value = num2str(zef.filter_type);
62+
else
63+
zef.KF.filter_type.Value = '1';
64+
end
65+
66+
if isfield(zef,'number_of_ensembles')
67+
zef.KF.number_of_ensembles.Value = num2str(zef.number_of_ensembles);
68+
else
69+
zef.KF.number_of_ensembles.Value = '100';
70+
end
71+
72+
if isfield(zef,'normalize_data')
73+
zef.KF.normalize_data.Value = num2str(zef.normalize_data);
74+
else
75+
zef.KF.normalize_data.Value = 1;
76+
end
77+
78+
%set parameters if saved in ZI:
79+
%(Naming concept: zef.KF."field" = zef."field")
80+
zef_props = properties(zef.KF);
81+
for zef_i = 1:length(zef_props)
82+
if isfield(zef,zef_props{zef_i})
83+
zef.KF.(zef_props{zef_i}).Value = num2str(zef.(zef_props{zef_i}));
84+
elseif isprop(zef.KF.(zef_props{zef_i}),'Value')
85+
zef.(zef_props{zef_i}) = str2num(zef.KF.(zef_props{zef_i}).Value);
86+
end
87+
end
88+
clear zef_props zef_i
89+
90+
%_ Functions _
91+
zef.KF.inv_snr.ValueChangedFcn = 'zef.inv_snr = str2num(zef.KF.inv_snr.Value);';
92+
zef.KF.kf_evolution_prior_mode.ValueChangedFcn = 'zef.kf_evolution_prior_mode = str2num(zef.KF.kf_evolution_prior_mode.Value);';
93+
zef.KF.normalize_data.ValueChangedFcn = 'zef.normalize_data = zef.KF.normalize_data.Value;';
94+
zef.KF.inv_sampling_frequency.ValueChangedFcn = 'zef.inv_sampling_frequency = str2num(zef.KF.inv_sampling_frequency.Value);';
95+
zef.KF.inv_low_cut_frequency.ValueChangedFcn = 'zef.inv_low_cut_frequency = str2num(zef.KF.inv_low_cut_frequency.Value);';
96+
zef.KF.inv_high_cut_frequency.ValueChangedFcn = 'zef.inv_high_cut_frequency = str2num(zef.KF.inv_high_cut_frequency.Value);';
97+
zef.KF.inv_time_1.ValueChangedFcn = 'zef.inv_time_1 = str2num(zef.KF.inv_time_1.Value);';
98+
zef.KF.inv_time_2.ValueChangedFcn = 'zef.inv_time_2 = str2num(zef.KF.inv_time_2.Value);';
99+
zef.KF.number_of_frames.ValueChangedFcn = 'zef.number_of_frames = str2num(zef.KF.number_of_frames.Value);';
100+
zef.KF.inv_time_3.ValueChangedFcn = 'zef.inv_time_3 = str2num(zef.KF.inv_time_3.Value);';
101+
zef.KF.kf_smoothing.ValueChangedFcn = 'zef.kf_smoothing = str2num(zef.KF.kf_smoothing.Value);';
102+
zef.KF.StartButton.ButtonPushedFcn = 'zef = zef_KF(zef);';
103+
zef.KF.ApplyButton.ButtonPushedFcn = 'zef_props = properties(zef.KF); for zef_i = 1:length(zef_props); if isprop(zef.KF.(zef_props{zef_i}),''Value''); zef.(zef_props{zef_i}) = str2num(zef.KF.(zef_props{zef_i}).Value); end; end; clear zef_props zef_i;';
104+
zef.KF.CloseButton.ButtonPushedFcn = 'delete(zef.KF);';
105+
zef.KF.filter_type.ValueChangedFcn = 'zef.filter_type = str2num(zef.KF.filter_type.Value);';
106+
107+
%set fonts
108+
set(findobj(zef.KF.UIFigure.Children,'-property','FontSize'),'FontSize',zef.font_size);

‎+plugins/+ClassKF/zef_kf_sLORETA_OLD.m

+8
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
function [P, L, Q] = zef_kf_sLORETA_OLD(L, Q, std_lhood, theta0)
2+
% Sparse
3+
D2 = inv(sparse(diag(diag(L' / (L * L' + std_lhood^2 / theta0 * eye(size(L,1))) * L))));
4+
P = theta0 * D2;
5+
L = L / chol(D2);
6+
Q = D2*Q;
7+
end
8+

‎+plugins/+ClassKF/zef_kf_start.m

+51
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
%This is the startup script for KF app. One must add this as launch
2+
%script to zeffiro_plugins file:
3+
%Beamformer, inverse_tools, zef_KF_start
4+
5+
zef.KF = zef_kf_app;
6+
7+
%_ Initial values _
8+
zef.KF.inv_snr.Value = '30';
9+
zef.KF.inv_sampling_frequency.Value = '1025';
10+
zef.KF.inv_low_cut_frequency.Value = '7';
11+
zef.KF.inv_high_cut_frequency.Value = '9';
12+
zef.KF.inv_time_1.Value = '0';
13+
zef.KF.inv_time_2.Value = '0';
14+
zef.KF.number_of_frames.Value = '0';
15+
zef.KF.inv_time_3.Value = '0';
16+
zef.KF.filter_type.Value = '1';
17+
zef.KF.number_of_ensembles.Value = '100';
18+
zef.KF.kf_number_of_noise_steps.Value = '4';
19+
20+
21+
%set parameters if saved in ZI:
22+
%(Naming concept: zef.KF."field" = zef."field")
23+
zef_props = properties(zef.KF);
24+
for zef_i = 1:length(zef_props)
25+
if isfield(zef,zef_props{zef_i})
26+
zef.KF.(zef_props{zef_i}).Value = num2str(zef.(zef_props{zef_i}));
27+
elseif isprop(zef.KF.(zef_props{zef_i}),'Value')
28+
zef.(zef_props{zef_i}) = str2num(zef.KF.(zef_props{zef_i}).Value);
29+
end
30+
end
31+
clear zef_props zef_i
32+
33+
%_ Functions _
34+
zef.KF.inv_snr.ValueChangedFcn = 'zef.inv_snr = str2num(zef.KF.inv_snr.Value);';
35+
zef.KF.inv_sampling_frequency.ValueChangedFcn = 'zef.inv_sampling_frequency = str2num(zef.KF.inv_sampling_frequency.Value);';
36+
zef.KF.inv_low_cut_frequency.ValueChangedFcn = 'zef.inv_low_cut_frequency = str2num(zef.KF.inv_low_cut_frequency.Value);';
37+
zef.KF.inv_high_cut_frequency.ValueChangedFcn = 'zef.inv_high_cut_frequency = str2num(zef.KF.inv_high_cut_frequency.Value);';
38+
zef.KF.inv_time_1.ValueChangedFcn = 'zef.inv_time_1 = str2num(zef.KF.inv_time_1.Value);';
39+
zef.KF.inv_time_2.ValueChangedFcn = 'zef.inv_time_2 = str2num(zef.KF.inv_time_2.Value);';
40+
zef.KF.number_of_frames.ValueChangedFcn = 'zef.number_of_frames = str2num(zef.KF.number_of_frames.Value);';
41+
zef.KF.kf_number_of_noise_steps.ValueChangedFcn = 'zef.kf_number_of_noise_steps = str2num(zef.KF.kf_number_of_noise_steps.Value);';
42+
zef.KF.number_of_ensembles.ValueChangedFcn = 'zef.number_of_ensembles = str2num(zef.KF.number_of_ensembles.Value);';
43+
zef.KF.inv_time_3.ValueChangedFcn = 'zef.inv_time_3 = str2num(zef.KF.inv_time_3.Value);';
44+
zef.KF.kf_smoothing.ValueChangedFcn = 'zef.kf_smoothing = str2num(zef.KF.kf_smoothing.Value);';
45+
zef.KF.StartButton.ButtonPushedFcn = 'zef = zef_KF(zef);';
46+
zef.KF.ApplyButton.ButtonPushedFcn = 'zef_props = properties(zef.KF); for zef_i = 1:length(zef_props); if isprop(zef.KF.(zef_props{zef_i}),''Value''); zef.(zef_props{zef_i}) = str2num(zef.KF.(zef_props{zef_i}).Value); end; end; clear zef_props zef_i;';
47+
zef.KF.CloseButton.ButtonPushedFcn = 'delete(zef.KF);';
48+
zef.KF.filter_type.ValueChangedFcn = 'zef.filter_type = str2num(zef.KF.filter_type.Value);';
49+
50+
%set fonts
51+
set(findobj(zef.KF.UIFigure.Children,'-property','FontSize'),'FontSize',zef.font_size);

0 commit comments

Comments
 (0)
Please sign in to comment.