s2let_demo_covariance

 s2let_demo_covariance - Run covariance demo.

 Demo to compare theoretical covariance of wavelet coefficients with
 empirical data from using our transform functions. The empirical
 covariance is computed for several sets of harmonic coefficients,
 and the theoretical covariance is compared to the mean of those
 calculations in units of its standard deviation.

 Default usage is given by

   s2let_demo_covariance

 Authors: Martin Büttner (m.buettner.d@gmail.com)
          Jason McEwen (www.jasonmcewen.org)

0001 % s2let_demo_covariance - Run covariance demo.
0002 %
0003 % Demo to compare theoretical covariance of wavelet coefficients with
0004 % empirical data from using our transform functions. The empirical
0005 % covariance is computed for several sets of harmonic coefficients,
0006 % and the theoretical covariance is compared to the mean of those
0007 % calculations in units of its standard deviation.
0008 %
0009 % Default usage is given by
0010 %
0011 %   s2let_demo_covariance
0012 %
0013 % Authors: Martin Büttner (m.buettner.d@gmail.com)
0014 %          Jason McEwen (www.jasonmcewen.org)
0015 
0016 % S2LET package to perform Wavelet transforms on the sphere.
0017 % Copyright (C) 2012-2015-2014 Boris Leistedt and Jason McEwen
0018 % See LICENSE.txt for license details
0019 
0020 clear all;
0021 
0022 % Define parameters
0023 B = 2;
0024 L = 16;
0025 N = L;
0026 Spin = 0;
0027 J_min = 2;
0028 J = s2let_jmax(L, B);
0029 
0030 var_flm = 1; % Should we use the actual variance var(flmn) of each
0031              % realization here instead?
0032 
0033 % Compute theoretical covariance of wavelet coefficients.
0034 % The covariance <Wj(rho)Wj*(rho)> is given by the following expression:
0035 %
0036 % sigma^2 Sum(l,n) |Psij_ln|^2
0037 %
0038 % Where sigma^2 is the variance of the harmonic coefficients and Psij_lm
0039 % are the harmonic coefficients of the j-th wavelet.
0040 %
0041 % A similar expression applies for the scaling function coefficients.
0042 
0043 [psi_lm, phi_l] = s2let_wavelet_tiling(B, L, N, Spin, J_min);
0044 
0045 covar_W_phi_theory = 0;
0046 for el = 0:L-1
0047     covar_W_phi_theory = covar_W_phi_theory + phi_l(el+1) * phi_l(el+1)';
0048 end
0049 
0050 covar_W_psi_theory = zeros(J-J_min+1,1);
0051 for j = J_min:J
0052     ind = 1;
0053     for el = 0:L-1
0054         for n = -el:el
0055             covar_W_psi_theory(j-J_min+1) = covar_W_psi_theory(j-J_min+1) + psi_lm(ind,j+1) * psi_lm(ind,j+1)';
0056             ind = ind + 1;
0057         end
0058     end
0059 end
0060 
0061 covar_W_phi_theory = covar_W_phi_theory .* var_flm
0062 covar_W_psi_theory = covar_W_psi_theory .* var_flm
0063 
0064 runs = 100;
0065 if J_min == 0
0066     % Ergodicity fails for J_min = 0, because the scaling function will
0067     % only cover f00. Hence <flm flm*> will be 0 in that case and the
0068     % scaling coefficients will all be the same. So, if we do have
0069     % J_min = 0, we take the variance over all realisations instead (of
0070     % course, we then won't have a standard deviation to compare it to the
0071     % theoretical variance).
0072     f_scals = zeros(runs,1);
0073 else
0074     covar_W_phi_data = zeros(1, runs);
0075 end
0076 
0077 covar_W_psi_data = zeros(J-J_min+1, runs);
0078 
0079 
0080 for i = 1:runs
0081     % Generate normally distributed random flmn of complex signal
0082     % with mean 0 and variance 1
0083     flm = zeros(L^2 - Spin^2, 1);
0084     flm = [zeros(Spin^2,1); (randn(size(flm)) + 1i*randn(size(flm)))/sqrt(2)*sqrt(var_flm)];
0085 
0086     % Compute inverse then forward transform.
0087     f = ssht_inverse(flm, L, 'Spin', Spin);
0088 
0089     [f_wav, f_scal] = s2let_transform_analysis_mw(complex(real(f), imag(f)), 'B', B, 'L', L, 'N', N,    ...
0090                                                   'Spin', Spin, 'J_min', J_min, ...
0091                                                   'Upsample', true);
0092 
0093     if J_min == 0
0094         f_scals(i) = f_scal(1);
0095     else
0096         covar_W_phi_data(i) = var(f_scal(:));
0097     end
0098 
0099     for j = J_min:J
0100         f_wav_j = f_wav{j-J_min+1};
0101         covar_W_psi_data(j-J_min+1,i) = var(f_wav_j(:));
0102     end
0103 end
0104 
0105 if J_min == 0
0106     covar_W_phi_data = var(f_scals(:))
0107     W_phi_error_absolute = abs(covar_W_phi_data - covar_W_phi_theory)
0108 else
0109     mean_covar_W_phi_data = mean(covar_W_phi_data)
0110     std_covar_W_phi_data = std(covar_W_phi_data)
0111     W_phi_error_in_std = abs(mean_covar_W_phi_data - covar_W_phi_theory)/std_covar_W_phi_data
0112 end
0113 
0114 mean_covar_W_psi_data = mean(covar_W_psi_data,2)
0115 std_covar_W_psi_data = std(covar_W_psi_data,0,2)
0116 W_psi_error_in_std = abs(mean_covar_W_psi_data - covar_W_psi_theory)./std_covar_W_psi_data