Home > src > main > matlab > flaglet_denoising_demo2.m

flaglet_denoising_demo2

PURPOSE ^

flaglet_denoising_demo2

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 flaglet_denoising_demo2
 Denoising example (large-scale structure data)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % flaglet_denoising_demo2
0002 % Denoising example (large-scale structure data)
0003 
0004 % Resolution parameters
0005 L = 128 ;
0006 P = 128 ;
0007 
0008 % Load LSS data
0009 load('lss_flmp_312')
0010 R = 420;
0011 
0012 f_lmp = zeros(P, L^2);
0013 f_lmp(:,:) = f_lmp_full(1:P, 1:L^2);
0014 
0015 disp('Flag_synthesis...')
0016 f = flag_synthesis(f_lmp, 'R', R, 'reality', true);
0017 
0018 % Wavelet decomposition parameters
0019 Downsample = true;
0020 B_l = 2 ;
0021 B_p = 2 ;
0022 J_min_l = 0 ;
0023 J_min_p = 0 ;
0024 J_l = ceil(log(L) ./ log(B_l))
0025 J_p = ceil(log(P) ./ log(B_p))
0026 
0027 disp('Compute tilling of harmonic space...')
0028 [kappa kappa0] = flaglet_axisym_tilling(B_l, B_p, L, P, J_min_l, J_min_p);
0029 
0030 disp('Find noise covariance from SNR')
0031 SNR_ini = 8.0;
0032 f_power = sum(sum(abs(f_lmp).^2));
0033 noise_power = f_power * ( 10.0^(-SNR_ini/10.0) );
0034 factor = sum( ((1:P)./P).^2 ) * L * L ;
0035 sigma_noise = sqrt( noise_power / factor )
0036 
0037 disp('Generate random noise n_lmp of std dev sigma')
0038 n_lmp = sigma_noise * randn(P, L^2);
0039     for en = 1:P
0040       for el = 0:L-1
0041         ind = ssht_elm2ind(el, 0);
0042         n_lmp(en,ind) = sigma_noise .* randn * (en/P) ;
0043         for m = 1:el
0044           ind = ssht_elm2ind(el, m);
0045           n_lmp(en,ind) = ...
0046             sqrt(sigma_noise^2./2) .* randn * (en/P)  ...
0047             + sqrt(-1) * sqrt(sigma_noise^2./2) .* randn * (en/P)  ;
0048         end
0049       end
0050     end
0051     
0052 noise_power = sum(sum(abs(n_lmp).^2));
0053 SNR_before = 10*log10( f_power / noise_power )
0054 
0055 disp('Construct the noise in real space...')
0056 initial_noise = flag_synthesis(n_lmp, 'R', R, 'Reality', true);
0057 
0058 disp('Add the noise to the initial dataset')
0059 %g_lmp = f_lmp + n_lmp ;
0060 g = f + initial_noise ;
0061 
0062 disp('Perform wavelet decomposition...')
0063 [g_wav, g_scal] = flaglet_axisym_analysis(g, 'B_l', B_l, 'B_p', B_p, 'L', L, 'P', P, 'J_min_l', J_min_l, 'J_min_p', J_min_p, 'Reality', true, 'Downsample', Downsample);
0064 
0065 disp('Construct noise model...')
0066 noisemodel = cell(J_l+1, J_p+1);
0067 for jl = J_min_l:J_l
0068     if Downsample == true
0069         bandlimit_l = min([ ceil(B_l^(jl+1)) L ]);
0070     else
0071         bandlimit_l = L;
0072     end
0073     for jp = J_min_p:J_p
0074         if Downsample == true
0075             bandlimit_p = min([ ceil(B_p^(jp+1)) P ]);
0076         else
0077             bandlimit_p = P;
0078         end
0079         [jl, jp]
0080         temp = kappa{jl+1,jp+1};  
0081         contrib = zeros(1, bandlimit_p);
0082         nodes = slag_sampling(bandlimit_p, R);
0083         for r = 0:bandlimit_p-1
0084             birs = [nodes(r+1), R];
0085             for n = 0:bandlimit_p-1
0086                 fn = zeros(1,P);
0087                 fn(n+1) = 1.0;
0088                 [K_n_s, ~] = slag_synthesis(fn, 'Nodes', birs) ;
0089                 K_n_s = K_n_s(1);
0090                 for l = 0:bandlimit_l-1
0091                     contrib(r+1) = ...
0092                         contrib(r+1) + ((2*l+1)/(4*pi)) * (n/P) * (  K_n_s  * temp(n+1,l+1)).^2 ;
0093                 end
0094                   
0095             end
0096         end
0097         noisemodel{jl+1,jp+1} = sigma_noise * sqrt(contrib);
0098     end  
0099 end
0100 
0101 g_scal_rec = g_scal ;
0102 
0103 disp('Threshold wavelet coefficients...')
0104 g_wav_rec = cell(J_l+1-J_min_l, J_p+1-J_min_p);
0105 nb = 0;
0106 for jl = J_min_l:J_l
0107     if Downsample == true
0108         bandlimit_l = min([ ceil(B_l^(jl+1)) L ]);
0109     else
0110         bandlimit_l = L;
0111     end
0112     for jp = J_min_p:J_p
0113         if Downsample == true
0114             bandlimit_p = min([ ceil(B_p^(jp+1)) P ]);
0115         else
0116             bandlimit_p = P;
0117         end
0118         temp = g_wav{jl+1-J_min_l,jp+1-J_min_p};
0119         treshold = 3 * noisemodel{jl+1,jp+1};
0120         [jl, jp]
0121         for p = 1:2*bandlimit_l-1
0122             for t = 1:bandlimit_l
0123                 for r = 1:bandlimit_p
0124                     if abs(temp(r, t, p)) < treshold(r) % Hard thresholding
0125                         temp(r, t, p) = 0;
0126                         nb = nb + 1;
0127                     end
0128                 end
0129             end
0130         end
0131         g_wav_rec{jl+1-J_min_l,jp+1-J_min_p} = temp;
0132     end
0133 end
0134 
0135 disp('Reconstruct the denoised field from the wavelets...')
0136 g_denoised = flaglet_axisym_synthesis(g_wav_rec, g_scal_rec, 'B_l', B_l, 'B_p', B_p, 'L', L, 'P', P, 'J_min_l', J_min_l, 'J_min_p', J_min_p, 'Reality', true, 'Downsample', Downsample);
0137 g_lmp_denoised = flag_analysis(g_denoised, 'Reality', true);
0138 
0139 remaining_noise = g_denoised - f;
0140 f_power = sum(sum(abs(f_lmp).^2))
0141 
0142 noise_power = sum(sum(abs(n_lmp).^2))
0143 n_lmp_denoised = flag_analysis(remaining_noise, 'Reality', true);
0144 remaining_noise_power = sum(sum(abs(n_lmp_denoised).^2))
0145 
0146 SNR_before = 10*log10( f_power / noise_power )
0147 SNR_after = 10*log10( f_power / remaining_noise_power )
0148 
0149 
0150 % Shell to plot
0151 layer = 230;
0152 Lplot = 256;
0153 Pplot = 256;
0154 
0155 % Layout parameters
0156 minbound = 5; 
0157 maxbound = 100;
0158 zoomfactor = 2.2;
0159 
0160 disp('Plot...')
0161 figure('Position',[100 100 1000 850]) 
0162 
0163 subplot(2,2,1); 
0164 flaglet_plot_f( f, 'L', Lplot, 'P', Pplot, 'layer', layer   )
0165 zoom(zoomfactor)
0166 caxis([minbound maxbound])
0167 colormap(flipud(gray))
0168 
0169 subplot(2,2,2);
0170 flaglet_plot_f( initial_noise, 'L', Lplot, 'P', Pplot, 'layer', layer  )
0171 zoom(zoomfactor)
0172 caxis([minbound 0.6*maxbound])
0173 colormap(flipud(gray))
0174 
0175 subplot(2,2,3); 
0176 flaglet_plot_f( g, 'L', Lplot, 'P', Pplot, 'layer', layer   )
0177 zoom(zoomfactor)
0178 caxis([minbound 0.6*maxbound])
0179 colormap(flipud(gray))
0180 
0181 subplot(2,2,4); 
0182 flaglet_plot_f( g_denoised, 'L', Lplot, 'P', Pplot, 'layer', layer   )
0183 zoom(zoomfactor)
0184 caxis([minbound maxbound])
0185 colormap(flipud(gray))

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