Home > src > main > matlab > noise_prop_check_old.m

noise_prop_check_old

PURPOSE ^

noise_prop_check

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

noise_prop_check

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %noise_prop_check
0002 
0003 niter = 50;
0004 
0005 P = 8 ;
0006 L = 8 ;
0007 R = 420;
0008 
0009 sigma2_noise = 0.1 ;
0010 
0011 B_l = 2 ;
0012 B_n = 2 ;
0013 J_min_l = 0 ;
0014 J_min_n = 0 ;
0015 J_l = ceil(log(L) ./ log(B_l));
0016 J_n = ceil(log(P) ./ log(B_n));
0017 
0018 nodes = slag_sampling(P-1, R);
0019 % Compute all functions K_p(r_i) with p=1:P and i=1:P
0020 K_n_rs = zeros(P, P);
0021 for p = 0:P-1
0022     fp = zeros(1,P);
0023     fp(p+1) = 1.0;
0024     [temp, ~] = slag_synthesis(fp, 'R', R, 'Nodes', nodes) ;
0025     K_n_rs(p+1, :) = temp; 
0026 end
0027 
0028 % Compute coefficients K_{pp'} through gauss-laguerre quadrature
0029 K_pp = zeros(P, P);
0030 [nodes, weights] = slag_gausslaguerre_quadrature(P, 2);
0031 for p1 = 1:P
0032     for p2 = 1:P
0033         K_pp(p1, p2) = sum( K_n_rs(p1, :) .* K_n_rs(p2, :) .* weights ./ nodes.^2 );
0034     end
0035 end
0036 
0037 K_pp
0038 
0039 % Compute tilling and wavelet generating functions
0040 [kappa kappa0] = b3let_axisym_tilling(B_l, B_n, L, N, J_min_l, J_min_n);
0041 
0042 
0043 % Compute all functions K_p(r_i) with p=1:P and i=1:P+1
0044 K_n_rs = zeros(P, P+1);
0045 nodes = slag_sampling(P, R);
0046 for p = 0:P-1
0047     fp = zeros(1,P);
0048     fp(p+1) = 1.0;
0049     [temp, ~] = slag_synthesis(fp, 'R', R) ;
0050     K_n_rs(p+1, :) = temp; 
0051 end
0052 
0053 
0054 % Compute noise model
0055 sigma2_jj_rs = zeros(J_l+1, J_n+1, P+1);
0056 for jl = 0:J_l
0057     for jn = 0:J_n
0058         kappa_jj = kappa{jl+1,jn+1};
0059         for ell = 1:L
0060             for p1 = 1:P
0061                 for p2 = 1:P
0062                     factor = ((2*ell-1)/(4*pi)) * kappa_jj(p1, ell) * conj(kappa_jj(p2, ell)) * K_pp(p1, p2);
0063                     for ri = 1:P+1
0064                          sigma2_jj_rs(jl+1, jn+1, ri) = ...
0065                              sigma2_jj_rs(jl+1, jn+1, ri) + ...
0066                                 sigma2_noise * K_n_rs(p1, ri) * K_n_rs(p2, ri) * factor ; 
0067                     end
0068                 end
0069             end
0070         end
0071     end
0072 end
0073 
0074 % Compute noise realisation
0075 n_std_mean = zeros(J_l+1, J_n+1, P+1);
0076 n_lmp_var = zeros(P, L^2);
0077 
0078 for iter = 1:niter
0079     iter
0080     
0081     noise = zeros(P+1, L, 2*L-1);
0082     for ri = 0:P
0083         noise(ri+1, :, :) = sqrt(sigma2_noise) * randn(L, 2*L-1) / nodes(ri+1) ;
0084     end
0085     n_lmp = flag_analysis(noise, 'R', R, 'Reality', true);
0086     n_lmp_var = n_lmp_var + n_lmp.^2 / niter;
0087     
0088     [n_wav, n_scal] = b3let_axisym_analysis(noise, 'B_l', B_l, 'B_n', B_n, 'J_min_l', J_min_l, 'J_min_n', J_min_n, 'Reality', true);
0089 
0090     % Wavelets
0091     n_std = zeros(1,P+1);
0092     for jl = 0:J_l
0093         for jn = 0:J_n
0094             temp = n_wav{jl+1,jn+1};
0095             for r = 0:P
0096                 layer = temp(r+1,:,:);
0097                 n_std(r+1) = var3d(layer);
0098                 n_std_mean(jl+1, jn+1, r+1) = n_std_mean(jl+1, jn+1, r+1) + n_std(r+1) ./ niter ;
0099             end
0100            
0101         end
0102     end
0103         
0104 end
0105 
0106 sigma2_jj_rs ./ n_std_mean
0107 mean( mean ( mean ( sigma2_jj_rs ./ n_std_mean ) ) )

Generated on Mon 24-Sep-2012 12:26:33 by m2html © 2005