Description of s2let_demo_curvelet

0001 % s2let_demo_curvelet_covariance
0002 % Run curvelet covariance demo.
0003 %
0004 % Demo to compare theoretical covariance of curvelet coefficients with
0005 % empirical data from using our transform functions.
0006 % The empirical covariance is computed for several sets of
0007 % harmonic coefficients,
0008 % and the theoretical covariance is compared to the mean of those
0009 % calculations in units of its standard deviation.
0010 %
0011 % Default usage is given by
0012 %
0013 %   s2let_demo_curvelet_covariance
0014 %
0015 % ---------------------------------------------------------
0016 % S2LET package to perform Wavelet Transform on the Sphere.
0017 % Copyright (C) 2012-2016  Boris Leistedt, Jennifer Chan & Jason McEwen
0018 % See LICENSE.txt for license details
0019 % -----------------------------------------------------------
0020 clear all;
0021 
0022 % Define parameters
0023 B = 2;
0024 L = 16;
0025 Spin = 2;
0026 J_min = 2;
0027 J = s2let_jmax(L, B);
0028 
0029 var_flm = 1; 
0030 
0031 % Compute theoretical covariance of curvelet coefficients.
0032 % The covariance <Wj(rho)Wj*(rho)> is given by the following expression:
0033 %
0034 % sigma^2 Sum(l,n) |Psij_ln|^2
0035 %
0036 % Where sigma^2 is the variance of the harmonic coefficients and Psij_lm
0037 % are the harmonic coefficients of the j-th curvelet.
0038 %
0039 % A similar expression applies for the scaling function coefficients.
0040 
0041 
0042 % ---------------
0043 % Tile curvelets:
0044 % ---------------
0045 [cur_lm scal_l] = s2let_curvelet_tiling(B, L, J_min, ...
0046                                              'Spin', Spin, 'SpinLowered', false,...
0047                                              'SpinLoweredFrom',0);
0048 % reshape scaling functions
0049 kappa0_cur = zeros(L,1);
0050 for el = abs(Spin):L-1
0051  kappa0_cur(el+1,1) = scal_l(el^2+el+1,1);
0052 end
0053 % reshape curvelets
0054 kappa_cur = zeros(J+1,L^2);
0055 for j = J_min:J
0056  ind = Spin*Spin + 1;
0057  for el = abs(Spin):L-1
0058   for m= -el:el
0059    kappa_cur(j+1,ind) = cur_lm{j-J_min+1}(1,ind);
0060    ind = ind +1; 
0061   end
0062  end
0063 end 
0064 
0065 
0066 covar_W_phi_theory = 0;
0067 for el = 0:L-1
0068     covar_W_phi_theory = covar_W_phi_theory + kappa0_cur(el+1) * kappa0_cur(el+1)';
0069 end
0070 
0071 covar_W_psi_theory = zeros(J-J_min+1,1);
0072 for j = J_min:J
0073     ind = 1;
0074     for el = 0:L-1
0075         for n = -el:el
0076             covar_W_psi_theory(j-J_min+1) = covar_W_psi_theory(j-J_min+1) + kappa_cur(j+1,ind) * kappa_cur(j+1,ind)';
0077             ind = ind + 1;
0078         end
0079     end
0080 end
0081 
0082 covar_W_phi_theory = covar_W_phi_theory .* var_flm
0083 covar_W_psi_theory = covar_W_psi_theory .* var_flm
0084 
0085 runs = 100;
0086 if J_min == 0
0087     % Ergodicity fails for J_min = 0, because the scaling function will
0088     % only cover f00. Hence <flm flm*> will be 0 in that case and the
0089     % scaling coefficients will all be the same. So, if we do have
0090     % J_min = 0, we take the variance over all realisations instead (of
0091     % course, we then won't have a standard deviation to compare it to the
0092     % theoretical variance).
0093     f_scals = zeros(runs,1);
0094 else
0095     covar_W_phi_data = zeros(1, runs);
0096 end
0097 
0098 covar_W_psi_data = zeros(J-J_min+1, runs);
0099 
0100 
0101 for i = 1:runs
0102     % Generate normally distributed random flmn of complex signal
0103     % with mean 0 and variance 1
0104     flm = zeros(L^2 - Spin^2, 1);
0105     flm = [zeros(Spin^2,1); (randn(size(flm)) + 1i*randn(size(flm)))/sqrt(2)*sqrt(var_flm)];
0106 
0107     % Compute inverse then forward transform.
0108     f = ssht_inverse(flm, L, 'Spin', Spin);
0109 
0110     
0111     % ---------------
0112     % Curvelet transform:
0113     % ---------------
0114     [f_cur, f_scal] = s2let_transform_curvelet_analysis_px2cur(complex(real(f), imag(f)), ...
0115                                                                'B', B, 'L', L,  'J_min', J_min,  ...
0116                                                                'Spin', Spin,  ...
0117                                                                'Reality', false, ...
0118                                                                'Upsample', true, ...
0119                                                                'SpinLowered', false,...
0120                                                                'SpinLoweredFrom',0, ...
0121                                                                'Sampling','MW');
0122     if J_min == 0
0123         f_scals(i) = f_scal(1);
0124     else
0125         covar_W_phi_data(i) = var(f_scal(:));
0126     end
0127 
0128     for j = J_min:J
0129         f_cur_j = f_cur{j-J_min+1};
0130         covar_W_psi_data(j-J_min+1,i) = var(f_cur_j(:));
0131     end
0132 end
0133 
0134 if J_min == 0
0135     covar_W_phi_data = var(f_scals(:))
0136     W_phi_error_absolute = abs(covar_W_phi_data - covar_W_phi_theory)
0137 else
0138     mean_covar_W_phi_data = mean(covar_W_phi_data)
0139     std_covar_W_phi_data = std(covar_W_phi_data)
0140     W_phi_error_in_std = abs(mean_covar_W_phi_data - covar_W_phi_theory)/std_covar_W_phi_data
0141 end
0142 
0143 mean_covar_W_psi_data = mean(covar_W_psi_data,2)
0144 std_covar_W_psi_data = std(covar_W_psi_data,0,2)
0145 W_psi_error_in_std = abs(mean_covar_W_psi_data - covar_W_psi_theory)./std_covar_W_psi_data
s2let_demo_curvelet_covariance

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
