Home > src > main > matlab > flaglet_denoising_demo.m

flaglet_denoising_demo

PURPOSE ^

flaglet_denoising_demo

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 flaglet_denoising_demo
 Denoising example (geophysics data)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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