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)

- area_weighting AREA_WEIGHTING returns a matrix of weights for geodata based on latitude
- regionize % regionize
- resample_geodata RESAMPLE_GEODATA regrid the data to a new gridsize or fitting to new lat lons
- find_best_freq_set_anneal FIND_BEST_FREQ_SET_ANNEAL
- lse solves A*x = b for X in a least squares sense, given C*x = d

0001 % WEIGHTS Calculate the appropriate weight for each frequency. 0002 % 0003 % We do this by linear regression over all atmospheric cases. 0004 % 0005 % For each channel, we allow only those frequencies for which H is 0006 % nonzero to contribute. (For clear-sky, out of band frequencies 0007 % could probably also be used, but this certainly does not 0008 % generalize well to cloudy cases.) 0009 % 0010 % If there are negative weights, then these weights are set to 0011 % zero, and the regression is repeated without these channels. 0012 % 0013 % At the end, the weights are normalized so that the sum of all 0014 % weights for one channel is exactly 1, as done by Moncet et al., 0015 % 2008. This should improve robustness and generalizability. 0016 % 0017 % FORMAT W = weights(s, y_ref, H, y_mono) 0018 % 0019 % OUT W A matrix of channel weights. (Which can be used 0020 % instead of H, but without normalization) 0021 % 0022 % IN s Frequency set (type logical, dimension must match H) 0023 % y_ref Reference Tbs, result of H*y_mono 0024 % H Original measurement response matrix. We need this 0025 % to decide which frequencies should contribute to 0026 % which channel. 0027 % y_mono A batch of monochromatic Tbs (dimension must match 0028 % y_ref and s) 0029 0030 function W = weights(s, y_ref, H, y_mono) 0031 0032 s_y_ref = size(y_ref); 0033 s_y_mono = size(y_mono); 0034 s_H = size(H); 0035 0036 n_chan = s_y_ref(1); 0037 0038 if n_chan ~= s_H(1) 0039 error('y_ref and H not consistent!'); 0040 end 0041 0042 if s_y_ref(2) ~= s_y_mono(2) 0043 error('Atmospheric case number in *y_ref* and *y_mono* inconsistent!'); 0044 end 0045 0046 % Number of atmospheric cases: 0047 n_cases = s_y_ref(2); 0048 0049 if length(s) ~= s_y_mono(1) 0050 error('Frequency number in *s* and *y_mono* inconsistent!'); 0051 end 0052 0053 % Number of monochromatic frequencies: 0054 n_freqs = length(s); 0055 0056 % Loop channels 0057 W = H*0; 0058 for i=1:n_chan 0059 0060 % We need a loop here to ignore frequencies that would give negative 0061 % weights 0062 ignored = logical(zeros(1,n_freqs)); 0063 0064 all_ok = false; 0065 0066 while ~all_ok 0067 0068 % Find out, which of the active frequencies are relevant for this 0069 % channel. 0070 relevant = H(i,:) > 0; 0071 0072 % Subset of active frequencies for this channel. 0073 % They must be: 0074 % - part of the tested frequency set *s* 0075 % - part of the relevant frequencies for this channel 0076 % *relevant* 0077 % - not part of the frequencies that gave negative weights *ignored* 0078 now_active = s & relevant & ~ignored; 0079 0080 % Indices of the now active frequencies: 0081 i_now_active = find(now_active); 0082 0083 % Number of active frequencies: 0084 n_active_freqs = length(i_now_active); 0085 0086 % Construct design matrix (see multiple regression in Matlab online 0087 % help): 0088 %X = [ones(n_cases,1), y_mono(s,:)']; 0089 X = y_mono(now_active,:)'; 0090 0091 % Do least square fit for the weights, using backslash operator: 0092 w = X\y_ref(i,:)'; 0093 0094 % We do not want negative weights. We identify them, add them 0095 % to the ignore list, and repeat the regression. 0096 neg_weights = find(w<0); 0097 0098 if length(neg_weights) == 0 0099 all_ok = true; 0100 % disp('All freqs ok.') 0101 else 0102 ignored(i_now_active(neg_weights)) = 1; 0103 % disp('Ignoring some freqs.') 0104 end 0105 0106 end 0107 0108 % Construct matrix W, which we can use instead of matrix H. The 0109 % dimensionof W is still the same as that of H, so most elements 0110 % are zero! 0111 W(i,now_active) = w'; 0112 0113 % Set weight for those channels that we ignored to avoid negative 0114 % weights to zero: 0115 W(i,ignored) = 0; 0116 0117 % Normalize weight sum to 1: 0118 weight_sum = sum(W(i,now_active)); 0119 W(i,now_active) = W(i,now_active) / weight_sum; 0120 end 0121

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