Home > atmlab > gridcreation > optimize_f_grid_amsu.m

optimize_f_grid_amsu

PURPOSE ^

OPTIMIZE_F_GRID_AMSU

SYNOPSIS ^

function fi = optimize_f_grid_amsu(H, y_mono, acc)

DESCRIPTION ^

 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]

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

optimize_f_grid_amsu.m

SOURCE CODE ^

0001 % OPTIMIZE_F_GRID_AMSU
0002 %
0003 % Optimitation of the frequency grid for an ARTS calculation.
0004 %
0005 % This function is for AMSU-type instruments, i.e., instruments
0006 % with double sideband receivers.
0007 %
0008 % This method is very general. The only assumptions made are:
0009 % a) y_mono contains a batch of radiances for differen
0010 % atmospheric cases.
0011 % b) H*y_mono(i) gives a measurement vector y(i).
0012 %
0013 % The task is then to find a subset of frequency indices for
0014 % which the result is similar to the result for all frequencies.
0015 %
0016 % In contrast to optimize_za_grid, the optimization is done over a
0017 % whole batch of atmospheric cases, not just one. Another important
0018 % difference is that the accuray limit should be given in the same
0019 % unit as y_mono (usually Kelvin).
0020 %
0021 % Note the dopairs flag!
0022 %
0023 % FORMAT  fi = optimize_za_grid(H, y_mono, acc)
0024 %
0025 % OUT   fi           a set of frequency indices.
0026 % IN
0027 %       H            H matrix.
0028 %       y_mono  A set of radiance vectors, each vector with a
0029 %                    number of elements fitting the number of
0030 %                    columns of H.
0031 %       acc          Desired accuracy [same unit as y_mono]
0032 
0033 % 2008-09-01 Created by Stefan Buehler
0034 
0035 function fi = optimize_f_grid_amsu(H, y_mono, acc)
0036 
0037 % Important dimensions:
0038 
0039 x = size(H);
0040 nchannels = x(1);
0041 nfmono    = x(2);
0042 
0043 x = size(y_mono);
0044 if (x(1)~=nfmono)
0045   error('Dimensions of H and y_mono do not match!')
0046 end
0047 ncases = x(2);
0048 
0049 clear x;
0050 
0051 
0052 % Calculate the correct reference result:
0053 y_ref = H*y_mono;
0054 
0055 
0056 % The frequency index vector that we want to derive.
0057 fi = [];
0058 
0059 
0060 % The pool of frequencies to choose from.
0061 fpool = 1:nfmono;
0062 
0063 % We optimize the channels separately.
0064 for ichan=1:nchannels
0065 
0066   disp(sprintf('Channel %d.',ichan))
0067 
0068   % The accuracy we have already reached.
0069   this_acc = 1e99;
0070   
0071   % Do the following until accuracy is good enough
0072   while (this_acc>acc && length(fpool)>0 )
0073     
0074     % Calculate all possible results with frequency pairs added.
0075     rms_err = ones(length(fpool), length(fpool)) * 1e99; 
0076     max_err = ones(length(fpool), length(fpool)) * 1e99; 
0077     for i1=1:length(fpool)-1
0078       for i2=i1+1:length(fpool)
0079 
0080         % Add test frequency pairs.
0081         fi_test = [fi, fpool(i1), fpool(i2)];
0082 
0083         % fi_test must be properly sorted.
0084         fi_test = sort(fi_test);
0085         
0086         % Select the new H matrix for testing.
0087         H_test = H(ichan,fi_test);
0088 
0089         % FIXME: Remove?!
0090         %          nn = find(H_test ~= 0);
0091         %          H_test(nn) = H_test(nn) ./ H_test(nn);
0092         
0093         % H_test has to be correctly normalized. (This is important, if
0094         % we just remove frequencies, without re-normalizing, we can
0095         % never get a good result, unless we use all original
0096         % frequencies.)
0097         
0098         % The sum of H_test for each channel.
0099         H_norm = sum(H_test);
0100         
0101         % It is possible that the frequency we are checking here
0102         % is not relevant to this channel. In this case, we
0103         % should simply skip it. We can recognize this case by
0104         % H_norm being zero.
0105         if (H_norm==0)
0106           continue;
0107         end
0108 
0109         % Normalize by dividing by the sum.
0110         H_test = H_test / H_norm;
0111         
0112         % Calculate result with this set of frequencies.
0113         y_test = H_test * y_mono(fi_test,:);
0114         
0115         % Difference to reference case.
0116         diff = y_test - y_ref(ichan,:);
0117         
0118         rms_err(i1,i2) = rms(diff);
0119         max_err(i1,i2) = max(abs(diff(:)));          
0120       end
0121     end
0122 
0123     % Find the i1/i2 for which rms_err is minimal:
0124     [besti1, besti2] = find(rms_err == min(rms_err(:)));
0125 
0126     % If besti1 and besti2 are not scalars, it means that no
0127     % additional pair really improves things. (So we hit here all
0128     % the frequencies that are not relevant to the band.)
0129     if length(besti1)>1
0130       error(['We ran into a dead-end. No other frequency pair ' ...
0131              'improves this channel.'])
0132     end
0133     
0134     
0135     % Set this_acc for accuray test.
0136     this_acc = rms_err(besti1, besti2);
0137     
0138     % Add these indices to fi, and remove them from fpool:
0139     fi = [fi, fpool(besti1), fpool(besti2)];
0140     fpool = [fpool(1:besti1-1),...
0141              fpool(besti1+1:besti2-1),...
0142              fpool(besti2+1:end)];
0143     
0144     % Some output:
0145     disp(sprintf('%dth frequency pair added has indices %d/%d, rms = %g, max = %g.',...
0146                  length(fi)/2,...
0147                  fi(end-1), fi(end),...
0148                  rms_err(besti1, besti2),...
0149                  max_err(besti1, besti2)));
0150     
0151     %      H(ichan,fi)
0152     
0153     %      keyboard
0154   end                                   % While loop
0155 end                                     % Channel loop.
0156 
0157 % fi should be properly sorted.
0158 fi = sort(fi);
0159

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005