Home > src > main > matlab > noise_prop_check.m

noise_prop_check

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 = 20;
0004 
0005 N = 16 ;
0006 L = 16 ;
0007 R = 430.4;
0008 
0009 sigma_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(N) ./ log(B_n));
0017 
0018 [kappa kappa0] = b3let_axisym_tilling(B_l, B_n, L, N, J_min_l, J_min_n);
0019 
0020 nodes = slag_sampling(N, R);
0021 
0022 noisemodel0 = zeros(1,N+1);
0023 for r = 0:N
0024     birs = [nodes(r+1), R];
0025     for n = 0:N-1
0026         fn = zeros(1,N);
0027         fn(n+1) = 1.0;
0028         [K_n_s, ~] = slag_synthesis(fn, 'Nodes', birs) ;
0029         K_n_s = K_n_s(1);
0030         for l = 0:L-1
0031             noisemodel0(r+1) = ...
0032                 noisemodel0(r+1) + ((2*l+1)/(4*pi)) * (  K_n_s * kappa0(n+1,l+1)).^2 ;
0033         end
0034     end
0035 end
0036 noisemodel0 = sigma_noise * sqrt(noisemodel0);
0037 
0038 noisemodel = zeros(J_l+1, J_n+1, N+1);
0039 for jl = 0:J_l
0040     for jn = 0:J_n
0041         temp = kappa{jl+1,jn+1};  
0042         for r = 0:N
0043             birs = [nodes(r+1), R];
0044             for n = 0:N-1
0045                 fn = zeros(1,N);
0046                 fn(n+1) = 1.0;
0047                 [K_n_s, ~] = slag_synthesis(fn, 'Nodes', birs) ;
0048                 K_n_s = K_n_s(1);
0049                 for l = 0:L-1
0050                     noisemodel(jl+1,jn+1,r+1) = ...
0051                         noisemodel(jl+1,jn+1,r+1) + ((2*l+1)/(4*pi)) * (  K_n_s * temp(n+1,l+1)).^2 ;
0052                 end
0053                   
0054             end
0055         end
0056     end  
0057 end
0058 noisemodel = sigma_noise * sqrt(noisemodel);
0059 
0060 n_scal_std_mean = zeros(1,N+1); 
0061 n_std_mean = zeros(J_l+1, J_n+1, N+1);
0062 
0063 for iter = 1:niter
0064     iter
0065     
0066     % Generate random noise n_lmn of std dev sigma
0067     %n_lmn = sqrt(sigma_noise_ini) * randn(N,L^2);%randreal(N, L);
0068 
0069     n_lmn = sigma_noise * randn(N, L^2);
0070     for en = 1:N
0071       for el = 0:L-1
0072         ind = ssht_elm2ind(el, 0);
0073         n_lmn(en,ind) = sigma_noise .* randn ;
0074         for m = 1:el
0075           ind = ssht_elm2ind(el, m);
0076           n_lmn(en,ind) = ...
0077             sqrt(sigma_noise^2./2) .* randn  ...
0078             + sqrt(-1) * sqrt(sigma_noise^2./2) .* randn ;
0079         end
0080       end
0081     end
0082     sqrt(var2d(n_lmn));
0083     %n_lmn = n_lmn * sqrt(sigma_noise_ini/sigma_noise);
0084     %sigma_noise = var2d(n_lmn);
0085     noise = flag_synthesis(n_lmn, 'R', R, 'Reality', true);
0086     [n_wav, n_scal] = b3let_axisym_analysis(noise, 'B_l', B_l, 'B_p', B_n, 'L', L, 'P', N, 'J_min_l', J_min_l, 'J_min_p', J_min_n, 'Reality', true, 'downsample', false);
0087 
0088     % Scaling function
0089     n_scal_std = zeros(1,N); 
0090     for r = 0:N
0091         layer = n_scal(r+1,:,:);
0092         n_scal_std(r+1) = sqrt(var3d(layer));
0093     end
0094     n_scal_std_mean(:) = n_scal_std_mean(:) + n_scal_std(:) ./ niter ;
0095     
0096     % Wavelets
0097     n_std = zeros(N);
0098     for jl = 0:J_l
0099         for jn = 0:J_n
0100             temp = n_wav{jl+1,jn+1};
0101             for r = 0:N
0102                 layer = temp(r+1,:,:);
0103                 n_std(r+1) = sqrt(var3d(layer));
0104                 n_std_mean(jl+1, jn+1, r+1) = n_std_mean(jl+1, jn+1, r+1) + n_std(r+1) ./ niter ;
0105             end
0106            
0107         end
0108     end
0109         
0110 end
0111 
0112 noisemodel0 ./ n_scal_std_mean
0113 noisemodel ./ n_std_mean
0114 mean( mean ( mean ( noisemodel ./ n_std_mean ) ) )
0115

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