% WEIGHTS Calculate the appropriate weight for each frequency.
%
% We do this by linear regression over all atmospheric cases.
%
% For each channel, we allow only those frequencies for which H is
% nonzero to contribute. (For clear-sky, out of band frequencies
% could probably also be used, but this certainly does not
% generalize well to cloudy cases.)
%
% If there are negative weights, then these weights are set to
% zero, and the regression is repeated without these channels.
%
% At the end, the weights are normalized so that the sum of all
% weights for one channel is exactly 1, as done by Moncet et al.,
% 2008. This should improve robustness and generalizability.
%
% FORMAT W = weights(s, y_ref, H, y_mono)
%
% OUT W A matrix of channel weights. (Which can be used
% instead of H, but without normalization)
%
% IN s Frequency set (type logical, dimension must match H)
% y_ref Reference Tbs, result of H*y_mono
% H Original measurement response matrix. We need this
% to decide which frequencies should contribute to
% which channel.
% y_mono A batch of monochromatic Tbs (dimension must match
% y_ref and s)
function W = weights(s, y_ref, H, y_mono)
s_y_ref = size(y_ref);
s_y_mono = size(y_mono);
s_H = size(H);
n_chan = s_y_ref(1);
if n_chan ~= s_H(1)
error('y_ref and H not consistent!');
end
if s_y_ref(2) ~= s_y_mono(2)
error('Atmospheric case number in *y_ref* and *y_mono* inconsistent!');
end
% Number of atmospheric cases:
n_cases = s_y_ref(2);
if length(s) ~= s_y_mono(1)
error('Frequency number in *s* and *y_mono* inconsistent!');
end
% Number of monochromatic frequencies:
n_freqs = length(s);
% Loop channels
W = H*0;
for i=1:n_chan
% We need a loop here to ignore frequencies that would give negative
% weights
ignored = logical(zeros(1,n_freqs));
all_ok = false;
while ~all_ok
% Find out, which of the active frequencies are relevant for this
% channel.
relevant = H(i,:) > 0;
% Subset of active frequencies for this channel.
% They must be:
% - part of the tested frequency set *s*
% - part of the relevant frequencies for this channel
% *relevant*
% - not part of the frequencies that gave negative weights *ignored*
now_active = s & relevant & ~ignored;
% Indices of the now active frequencies:
i_now_active = find(now_active);
% Number of active frequencies:
n_active_freqs = length(i_now_active);
% Construct design matrix (see multiple regression in Matlab online
% help):
%X = [ones(n_cases,1), y_mono(s,:)'];
X = y_mono(now_active,:)';
% Do least square fit for the weights, using backslash operator:
w = X\y_ref(i,:)';
% We do not want negative weights. We identify them, add them
% to the ignore list, and repeat the regression.
neg_weights = find(w<0);
if length(neg_weights) == 0
all_ok = true;
% disp('All freqs ok.')
else
ignored(i_now_active(neg_weights)) = 1;
% disp('Ignoring some freqs.')
end
end
% Construct matrix W, which we can use instead of matrix H. The
% dimensionof W is still the same as that of H, so most elements
% are zero!
W(i,now_active) = w';
% Set weight for those channels that we ignored to avoid negative
% weights to zero:
W(i,ignored) = 0;
% Normalize weight sum to 1:
weight_sum = sum(W(i,now_active));
W(i,now_active) = W(i,now_active) / weight_sum;
end