%------------------------------------------------------------------------ % NAME: qp_rnd_sensor % % A function to generate the rand fields associated with a Q % structure defining the sensor. The fields needed in Q to % generate the rand fields are stated in /Qpack/README. % So far only measurement noise and baseline of order 0. % The rand fields are transformed as defined by the Hd % set by the Q struc, so if the transformation is into a % space of much smaller dimension, this results in calculations % using a significant smaller amount of memory. Use % qp_rnd_sensor if the fields want to be saved without applying % Hd. % % Use this function in a sequence % qp_rnd_atm_mem_red % qp_rnd_sensor_mem_red % qp_rnd_atmxsensor_mem_red % Information about the files produced in qp_rnd_sensor.m % % NOTE: this function require previous calculation of H % by running qp_H(Q); % % FORMAT: qp_rnd_sensor_mem_red( Q , [message, Qr )] % % OUT: - % IN: Q Setting structure. % OPTIONAL: message variable to print out a random set % name in the screen % Qr If Qr is given the rnd realizations are not done % but copied from Qr. This can be used e.g. to share % rnd realizations between Q and Qr when Q and % Qr have different sensor. % %------------------------------------------------------------------------ % % HISTORY: 2001.08.20 By Carlos Jimenez/ Patrick Eriksson % function qp_rnd_sensor_mem_red( Q , message, Qr) global REPORT_LEVEL r_l = Q.QP_LEVEL; %=== checking inputs if ~exist('message') message=''; end if ~exist('Qr') do_qr = 0;; else do_qr = 1; end %=== loading H matrix do_struc = 0; if isfield(Q,'LRED_ON') if Q.LRED_ON do_struc = 1; end end if ~do_struc load( [Q.OUT,'.hd'], '-mat' ); else eval(['load ',Q.OUT,'.hd ind -mat']); end %=== start printing out(1,1); if ~do_qr %=== creating thermal noise out(1,['Generating ',message,' rnd noise for sensor']); %=== loading sensor variables out(2,'Getting sensor data'); p_abs = read_datafile([Q.CALCGRIDS_DIR,'/',Q.P_ABS], 'VECTOR'); f_mono = read_datafile([Q.CALCGRIDS_DIR,'/',Q.F_MONO], 'VECTOR'); za_p = read_datafile([Q.CALCGRIDS_DIR,'/',Q.ZA_PENCIL],'VECTOR'); REPORT_LEVEL = 0; %== run to get f_sensor,f_y,za_sensor [dum,f_y,dum,f_sensor,za_sensor] = ... qp_H(Q,p_abs,f_mono,za_p,0); REPORT_LEVEL = r_l; %=== Number of realizations n_real = length(Q.NUMBER_DO); %== Setting seed of random generator rand('state',sum(100*clock)); out(2,'Doing measurement noise'); % -- getting covariance matrix % -- NOTE: Doing each za at one time to avoid a covariance % matrix of dimension f_y x f_y, even play with sparse % matrices lis very slow REPORT_LEVEL = 0; S = sFromFile([Q.SENSOR_DIR,'/',Q.MEASNOISE_COVMAT],f_sensor); REPORT_LEVEL = r_l; % -- generating realizations n_f = length(f_sensor); n_za = length(f_y)/n_f; if ~do_struc TN = []; else for g=1:length(ind) for l=1:ind(g) TN{g,l} = []; end end end for l=1:n_real for i= 1:n_za out(2,[' Zenith angle ',num2str(i)]); r = ((chol(S))' * randn(n_f,Q.NUMBER_DO(l))); rm = (mean(r'))'; rm = repmat(rm,1,size(r,2)); r = rm-r; if i==1 TNa = r; else TNa = [TNa;r]; end end % -- Transforming noise if ~do_struc eval(['load ',Q.OUT,'.hd Hd -mat']); TNa = Hd * TNa; TN = [TN TNa]; else for g=1:length(ind) is = num2str(g); for l=1:ind(g) js = num2str(l); eval(['load ',Q.OUT,'.hd Hd',is,js,' -mat']); eval(['H = Hd',is,js,';']) Yaux = H * TNa; TN{g,l} = [TN{g,l} Yaux]; end end end end save( [Q.OUT,'.measnoise'], 'TN' ); clear TN %=== creating constant baseline offset n_real = length(Q.NUMBER_DO); if Q.POLYFIT_DO > 0 if Q.POLYFIT_ORDER == 0 out(1,['Generating ',message,' baseline offset order 0 for sensor']); % -- getting covariance matrix REPORT_LEVEL = 0; S = sFromFile([Q.SENSOR_DIR,'/',sstring_get_i(Q.POLYFIT_COVMATS,1)],za_sensor); REPORT_LEVEL = r_l; n_za = length(za_sensor); if ~do_struc PO = []; else for g=1:length(ind) for l=1:ind(g) PO{g,l} = []; end end end for l=1:n_real % -- random values per za ra = (chol(S))' * randn( n_za, Q.NUMBER_DO(l) ); % -- extending to whole freq vector POa = zeros( n_za * n_f, Q.NUMBER_DO(l) ); for j = 1:n_za out(2,[' Zenith angle ',num2str(j)]); inda = (1 + ( j -1 ) * n_f ) : ( j * n_f); POa(inda,:) = repmat( ra(j,:), n_f, 1); end if 0 for kk = 1:size(POa,2) POa(:,kk) = rem_wave(POa(:,kk),f_sensor)'; end end % -- Transforming baseline if ~do_struc eval(['load ',Q.OUT,'.hd Hd -mat']); POa = Hd * POa; PO = [PO POa]; else for g=1:length(ind) is = num2str(g); for l=1:ind(g) js = num2str(l); eval(['load ',Q.OUT,'.hd Hd',is,js,' -mat']); eval(['H = Hd',is,js,';']) Yaux = H * POa; PO{g,l} = [PO{g,l} Yaux]; end end end end % -- saving polyfit save( [Q.OUT,'.polyfit'], 'PO' ); clear PO else disp('Sorry, baseline offsets of order > 0 not implemented yet.') end end else %do_qr out(1,['Copying ',message,' rnd noise for sensor']); load([Qr.OUT,'.measnoise'], '-mat' ); save( [Q.OUT,'.measnoise'], 'TN' ); clear TN if Q.POLYFIT_DO > 0 & Q.POLYFIT_ORDER == 1 load([Qr.OUT,'.polyfit'], '-mat' ); save( [Q.OUT,'.polyfit'], 'PO' ); clear PO end end %do_qr out(1,-1);