% OPTIMIZE_F_GRID_AMSU
%
% Optimitation of the frequency grid for an ARTS calculation.
%
% This function is for AMSU-type instruments, i.e., instruments
% with double sideband receivers.
%
% This method is very general. The only assumptions made are:
% a) y_mono contains a batch of radiances for differen
% atmospheric cases.
% b) H*y_mono(i) gives a measurement vector y(i).
%
% The task is then to find a subset of frequency indices for
% which the result is similar to the result for all frequencies.
%
% In contrast to optimize_za_grid, the optimization is done over a
% whole batch of atmospheric cases, not just one. Another important
% difference is that the accuray limit should be given in the same
% unit as y_mono (usually Kelvin).
%
% Note the dopairs flag!
%
% FORMAT fi = optimize_za_grid(H, y_mono, acc)
%
% OUT fi a set of frequency indices.
% IN
% H H matrix.
% y_mono A set of radiance vectors, each vector with a
% number of elements fitting the number of
% columns of H.
% acc Desired accuracy [same unit as y_mono]
% 2008-09-01 Created by Stefan Buehler
function fi = optimize_f_grid_amsu(H, y_mono, acc)
% Important dimensions:
x = size(H);
nchannels = x(1);
nfmono = x(2);
x = size(y_mono);
if (x(1)~=nfmono)
error('Dimensions of H and y_mono do not match!')
end
ncases = x(2);
clear x;
% Calculate the correct reference result:
y_ref = H*y_mono;
% The frequency index vector that we want to derive.
fi = [];
% The pool of frequencies to choose from.
fpool = 1:nfmono;
% We optimize the channels separately.
for ichan=1:nchannels
disp(sprintf('Channel %d.',ichan))
% The accuracy we have already reached.
this_acc = 1e99;
% Do the following until accuracy is good enough
while (this_acc>acc && length(fpool)>0 )
% Calculate all possible results with frequency pairs added.
rms_err = ones(length(fpool), length(fpool)) * 1e99;
max_err = ones(length(fpool), length(fpool)) * 1e99;
for i1=1:length(fpool)-1
for i2=i1+1:length(fpool)
% Add test frequency pairs.
fi_test = [fi, fpool(i1), fpool(i2)];
% fi_test must be properly sorted.
fi_test = sort(fi_test);
% Select the new H matrix for testing.
H_test = H(ichan,fi_test);
% FIXME: Remove?!
% nn = find(H_test ~= 0);
% H_test(nn) = H_test(nn) ./ H_test(nn);
% H_test has to be correctly normalized. (This is important, if
% we just remove frequencies, without re-normalizing, we can
% never get a good result, unless we use all original
% frequencies.)
% The sum of H_test for each channel.
H_norm = sum(H_test);
% It is possible that the frequency we are checking here
% is not relevant to this channel. In this case, we
% should simply skip it. We can recognize this case by
% H_norm being zero.
if (H_norm==0)
continue;
end
% Normalize by dividing by the sum.
H_test = H_test / H_norm;
% Calculate result with this set of frequencies.
y_test = H_test * y_mono(fi_test,:);
% Difference to reference case.
diff = y_test - y_ref(ichan,:);
rms_err(i1,i2) = rms(diff);
max_err(i1,i2) = max(abs(diff(:)));
end
end
% Find the i1/i2 for which rms_err is minimal:
[besti1, besti2] = find(rms_err == min(rms_err(:)));
% If besti1 and besti2 are not scalars, it means that no
% additional pair really improves things. (So we hit here all
% the frequencies that are not relevant to the band.)
if length(besti1)>1
error(['We ran into a dead-end. No other frequency pair ' ...
'improves this channel.'])
end
% Set this_acc for accuray test.
this_acc = rms_err(besti1, besti2);
% Add these indices to fi, and remove them from fpool:
fi = [fi, fpool(besti1), fpool(besti2)];
fpool = [fpool(1:besti1-1),...
fpool(besti1+1:besti2-1),...
fpool(besti2+1:end)];
% Some output:
disp(sprintf('%dth frequency pair added has indices %d/%d, rms = %g, max = %g.',...
length(fi)/2,...
fi(end-1), fi(end),...
rms_err(besti1, besti2),...
max_err(besti1, besti2)));
% H(ichan,fi)
% keyboard
end % While loop
end % Channel loop.
% fi should be properly sorted.
fi = sort(fi);