Home > src > main > matlab > denoising_demo3.m

denoising_demo3

PURPOSE ^

denoising_demo3

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 denoising_demo3
 Denoising example (geophysics data) and plot all wavelets

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % denoising_demo3
0002 % Denoising example (geophysics data) and plot all wavelets
0003 
0004 % Resolution parameters
0005 L = 92 ;
0006 P = 92 ;
0007 R = 3000;
0008 
0009 % Layout parameters
0010 minbound = -2; 
0011 maxbound = 2;
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 = 5 ;
0043 B_p = 5 ;
0044 J_min_l = 2 ;
0045 J_min_p = 2 ;
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 of f...')
0085 [f_wav, f_scal] = flaglet_axisym_analysis(f, '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('Perform wavelet decomposition of g...')
0088 [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);
0089 
0090 
0091 % Layout parameters
0092 minbound = -1.0; 
0093 maxbound = 1.0;
0094 Lplot = 192;
0095 Pplot = 192;
0096 layer = 173;
0097 zoomfactor = 2.0;
0098 
0099 nx = 2;
0100 ny = 2;
0101 
0102 figure('Position',[100 100 1000 850]) 
0103 ind = 1;
0104 subplot(ny, nx, ind);
0105 flaglet_plot_f( f, 'L', Lplot, 'P', Pplot, 'layer', layer )
0106 zoom(zoomfactor)
0107 caxis([minbound maxbound])
0108 title('Data')
0109 
0110 ind = 2;
0111 subplot(ny, nx, ind);
0112 flaglet_plot_f( g, 'L', Lplot, 'P', Pplot, 'layer', layer )
0113 zoom(zoomfactor)
0114 caxis([minbound maxbound])
0115 title('Noisy data')
0116 
0117 ind = 3;
0118 subplot(ny, nx, ind);
0119 flaglet_plot_f( f_scal, 'L', Lplot, 'P', Pplot, 'layer', layer )
0120 zoom(zoomfactor)
0121 caxis([minbound maxbound])
0122 title('Data')
0123 
0124 ind = 4;
0125 subplot(ny, nx, ind);
0126 flaglet_plot_f( g_scal, 'L', Lplot, 'P', Pplot, 'layer', layer )
0127 zoom(zoomfactor)
0128 caxis([minbound maxbound])
0129 title('Noisy data')
0130 
0131 
0132 nx = 2;
0133 ny = 4;
0134 
0135 J_l = [ 1 2 ];
0136 J_n = [ 1 2 ];
0137 
0138 figure('Position',[100 100 1000 1200])
0139 ind = 0;
0140 for jl = J_l
0141     for jn = J_n
0142         ind = ind + 1; 
0143         subplot(ny, nx, ind);
0144         flaglet_plot_f( f_wav{jl, jn}, 'L', Lplot, 'P', Pplot, 'layer', layer )
0145         zoom(zoomfactor)
0146         title(['F : jl=',int2str(jl),' jn=',int2str(jn)])
0147         caxis([minbound maxbound])
0148         
0149         ind = ind + 1; 
0150         subplot(ny, nx, ind);
0151         flaglet_plot_f( g_wav{jl, jn}, 'L', Lplot, 'P', Pplot, 'layer', layer )
0152         zoom(zoomfactor)
0153         title(['G : jl=',int2str(jl),' jn=',int2str(jn)])
0154         caxis([minbound maxbound])
0155     end
0156 end

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