Initial add.

BlockCadzow.m

@@ -0,0 +1,27 @@
+function Z = BlockCadzow(Z, K, L, reality)
+
+n_iter = 100;
+
+for iter = 1:n_iter
+
+    ind = 0;
+    for m = -(L-1-K):(L-1-K)
+
+        % First average along the diagonals
+
+        Y = Z(ind+1:ind+(L-abs(m)-K), :);
+        [r, n] = size(Y);
+        for d = 0:(r+n-2)
+            M = MaskAdiag(r, n, d);
+            Y(logical(M)) = mean(Y(logical(M)));
+        end
+        Z(ind+1:ind+(L-abs(m)-K), :) = Y;
+        ind = ind + (L-abs(m)-K);
+    end
+
+    % Then do SVD thresholding
+
+    [U, S, V] = svd(Z);
+    S(:, K+1) = 0;
+    Z = U*S*V';
+end    

CmMatrix.m

@@ -0,0 +1,20 @@
+function C_m = CmMatrix(m, L_max, l_min)
+
+if nargin == 2
+    l_min = 0;
+end
+
+l_min = max(abs(m), l_min);
+l_m = l_min:L_max;
+C_m = zeros(length(l_m));
+
+for l = l_m
+    c_lm = PolypartCoeffs(l, abs(m), L_max);
+    N_lm = sqrt((2*l+1)*factorial(l-abs(m))/factorial(l+abs(m))/4/pi);
+
+    if (m < 0)
+        c_lm = (-1)^m * c_lm;
+    end
+
+    C_m(l-l_min+1, :) = N_lm * c_lm((abs(m)+1):(abs(m)+length(l_m)));
+end

CramerRaoBound.m

@@ -0,0 +1,30 @@
+function CRLB = CramerRaoBound(theta, phi, sigma2, t0, p0, a0, L, posonly)
+
+DfDt0 = zeros(size(theta(:)));
+DfDp0 = zeros(size(theta(:)));
+DfDa0 = zeros(size(theta(:)));
+
+for l = 0:(L-1)
+    for m = -l:l
+        Y_lm_n  = SphericalHarmonic(l, m, theta, phi);
+        Y_lm1_0 = SphericalHarmonic(l, m+1, t0, p0);
+        Y_lm_0  = SphericalHarmonic(l, m, t0, p0);
+
+        DfDa0 = DfDa0 + conj(Y_lm_0) * Y_lm_n;
+        DfDt0 = DfDt0 + a0 * (m*cot(t0)*conj(Y_lm_0) ...
+              + sqrt((l-m)*(l+m+1)) * exp(1i*p0) ...
+             .* conj(Y_lm1_0)) * Y_lm_n;
+        DfDp0 = DfDp0 + a0 * (-1i*m) * conj(Y_lm_0) * Y_lm_n;
+    end
+end
+
+if posonly
+    Df = [DfDt0 DfDp0].';
+else
+    Df = [DfDt0 DfDp0 DfDa0].';
+end
+
+
+% Df   = [DfDt0 DfDp0].';
+I    = 1/sigma2 * (Df * Df');
+CRLB = inv(I);

DiracSpectrum.m

@@ -0,0 +1,11 @@
+function f = DiracSpectrum(S, L)
+
+i = 0;
+f = zeros(L^2, 1);
+for l = 0:L-1
+    for m = -l:l
+        i = i + 1;
+        Y    = conj(SphericalHarmonic(l, m, S(:, 2), S(:, 3)));
+        f(i) = sum(S(:, 1) .* Y);
+    end
+end

Example10_CRLB.m

@@ -0,0 +1,120 @@
+%--------------------------------------------------------------------------
+% CRLB example for 1 Dirac
+%--------------------------------------------------------------------------
+
+% Dirac parameters
+t0_list = [pi/2 pi/4];
+p0_list = [1 1];
+a0_list  = [1 0.1];
+
+% Sampling grid
+L = 2;
+[T, P] = ssht_sampling(L, 'Grid', true);
+
+% Simulation parameters
+SNR_range = -10:1:50;
+n_iter    = 100;
+
+
+% Compute the samples of the bandlimited Dirac
+
+% Simulate and plot
+for ifig = 1:2
+    h = figure(ifig);
+    col = summer;
+    col = col(1:32:end, :);
+    set(h,'NextPlot','replacechildren')
+    set(h,'DefaultAxesColorOrder', col);
+    clf; hold all;
+end
+
+matlabpool(3);
+
+var_out = zeros(length(t0_list), length(SNR_range));
+var_theta = length(t0_list), length(SNR_range);
+var_noise = length(t0_list), length(SNR_range);
+CRLB = length(t0_list), length(SNR_range);
+CRLB_theta = length(t0_list), length(SNR_range);
+
+for i_pos = 1:length(t0_list)
+    t0 = t0_list(i_pos);
+    p0 = p0_list(i_pos);
+    a0 = a0_list(i_pos);
+
+    flm = DiracSpectrum([a0, t0, p0], L);
+    f   = real(ssht_inverse(flm, L));
+
+    % Signal power
+    P_signal = var(f(:));
+
+
+    parfor i_SNR = 1:length(SNR_range)
+        SNR = SNR_range(i_SNR)
+        estimates = zeros(n_iter, 3);
+        for iter = 1:n_iter
+            % Get the samples of the bandlimited Dirac
+            P_noise = P_signal / 10^(SNR/10);
+            noise   = randn(size(f));
+            noise   = sqrt(P_noise / var(noise(:))) * noise;
+
+            f_noise = f + noise;
+            flm_noise = ssht_forward(f_noise, L, 'Reality', true);
+            [t_est, p_est, a_est] = SphereFRI(flm_noise, 1, false, true);
+            estimates(iter, :) = [t_est, p_est, a_est];
+        end
+        var_out(i_pos, i_SNR) = sum(sum(bsxfun(@minus, real(estimates(:, 1:2)), [t0 p0]).^2))/n_iter; 
+        var_theta(i_pos, i_SNR) = sum(real(estimates(:, 1) - t0).^2)/n_iter; 
+        var_noise(i_pos, i_SNR) = P_noise;
+
+        F_inv = CramerRaoBound(T(:), P(:), P_noise, t0, p0, a0, L, 1);
+        CRLB(i_pos, i_SNR) = (F_inv(1,1) + F_inv(2,2));
+        CRLB_theta(i_pos, i_SNR) = F_inv(1,1);
+    end
+
+    % Here you do the variance of the estimator... You should instead do the
+    % MSE. They interestingly coincide... so unbiased...
+
+    % Figure 1
+    figure(1);
+    plot(SNR_range, 10*log10(var_out(i_pos, :)), ... 
+        ':', ...
+        'LineWidth', 1, ...
+        'MarkerSize', 3, ...
+        'MarkerFaceColor', [0.5,0.5,0.5], ...
+        'MarkerEdgeColor', 'k');
+    plot(SNR_range, 10*log10(CRLB(i_pos, :)), ...
+            '-', ...
+            'LineWidth', 1, ...
+            'MarkerSize', 3, ...
+            'MarkerFaceColor', [0.5,0.5,0.5], ...
+            'MarkerEdgeColor', 'k');
+
+    % Figure 2 (only theta)
+    figure(2);
+    plot(SNR_range, 10*log10(var_theta(i_pos, :)), ... 
+        ':', ...
+        'LineWidth', 1, ...
+        'MarkerSize', 3, ...
+        'MarkerFaceColor', [0.5,0.5,0.5], ...
+        'MarkerEdgeColor', 'k');
+    plot(SNR_range, 10*log10(CRLB_theta(i_pos, :)), ...
+            '-', ...
+            'LineWidth', 1, ...
+            'MarkerSize', 3, ...
+            'MarkerFaceColor', [0.5,0.5,0.5], ...
+            'MarkerEdgeColor', 'k');
+    filename = sprintf('CRLB_data_%d.mat', i_pos);
+    save(filename);
+end
+
+xlabel('Input SNR [dB]');
+ylabel('10log10(Output MSE)');
+legend('MSE, \pi/2', 'CRLB, \pi/2', 'MSE, \pi/4', 'CRLB, \pi/4');
+axis tight;
+
+figWidth = 3.5;
+ExportifyFigure(h, [figWidth, 1/2*figWidth]);
+
+filename = sprintf('../papers/TSP/Figures/CRLB_plots.pdf');
+% export_fig('format','eps', filename);
+

Example1_Diracs.m

@@ -0,0 +1,19 @@
+% Setup some parameteres
+K = 5;
+L = ceil(K + sqrt(K));
+
+% Generate a random stream of Diracs, and compute its spectrum
+S = RandomDiracs(K, 1);
+
+f = DiracSpectrum(S, L);
+
+% Reconstruct the signal using SphereFRI
+[t, p, a] = SphereFRI(f, K);
+
+% Error
+[~, idx1] = sort(t);
+[~, idx2] = sort(S(:, 2));
+
+e = [(t(idx1)-S(idx2,2)) (p(idx1)-S(idx2,3)) (a(idx1)-S(idx2,1))]
+
+sum(abs(e).^2)

Example2_Green.m

@@ -0,0 +1,106 @@
+%==========================================================================
+%  This example plots Helmholtz Green's functions on the sphere
+%==========================================================================
+
+% NOTE: It would be more efficient to directly compute the spectral
+% coefficients, but this way we are trying to emulate what's happpening in
+% the real world.
+
+% Some parameters
+L = 100;            % Cutoff bandwidth for spectral computation (0...L-1)
+Lcutoff = 30;       % Cutoff bandwidth for Green's function computation (0...L-1)
+r = 0.25;           % Rigid sphere radius
+f = 500;            % Frequency in Hertz
+k = 2*pi*f/343;     % Wavenumber (== 2*pi*f/c)
+r_s = [0 0 1;
+       0 0 2; 
+       0 0 3; 
+       0 0 4; 
+       0 0 5]';
+
+% Sampling grid
+[t, p] = ssht_sampling(L, 'Grid', true);
+[x, y, z] = SphCart(t, p, r*ones(size(t)));
+
+g = zeros([size(t, 1) size(t, 2) size(r_s, 2)]);
+g_lm = zeros(L^2, size(r_s, 2));
+
+for i = 1:size(r_s, 2)
+    % Compute the Green's function
+    disp('Computing Green''s Function...');
+    g_i = GreenRigidSphere([x(:) y(:) z(:)]', r_s(:, i), k, Lcutoff);
+    g(:, :, i) = reshape(g_i, size(t));P
+
+    % And its Spherical Harmonic Transform
+    disp('Computing SHT...');
+    g_lm(:, i) = ssht_forward(g(:, :, i), L);
+end
+
+
+%% Do the plots
+
+% Green's function on the sphere
+figure(1);
+subplot(2,2,1);
+ssht_plot_sphere(abs(g(:, :, 1)), L, 'Type', 'colour', ...
+    'ColourBar', true, 'Lighting', true);
+
+subplot(2,2,2);
+ssht_plot_sphere(abs(g(:, :, 2)), L, 'Type', 'colour', ...
+    'ColourBar', true, 'Lighting', true);
+
+subplot(2,2,3);
+ssht_plot_sphere(abs(g(:, :, 3)), L, 'Type', 'colour', ...
+    'ColourBar', true, 'Lighting', true);
+
+subplot(2,2,4);
+ssht_plot_sphere(abs(g(:, :, 4)), L, 'Type', 'colour', ...
+    'ColourBar', true, 'Lighting', true);
+
+
+% Green's function spectrum
+Lplot = 20;
+l = 0:Lplot;
+l0 = l.^2 + l + 1;
+
+figure(2);
+clf;
+hold all;
+
+stem(0:length(l0)-1, abs(g_lm(l0, 1)), 'bv', 'filled', 'MarkerSize', 6, 'LineWidth', 1.2);
+stem(0:length(l0)-1, abs(g_lm(l0, 2)), 'rs', 'filled', 'MarkerSize', 6, 'LineWidth', 1.2);
+stem(0:length(l0)-1, abs(g_lm(l0, 3)), 'ko', 'filled', 'MarkerSize', 6, 'LineWidth', 1.2);
+stem(0:length(l0)-1, abs(g_lm(l0, 4)), 'mp', 'filled', 'MarkerSize', 6, 'LineWidth', 1.2);
+ax = axis;
+grid;
+xlabel('l');
+ylabel('|G(l, 0)|');
+
+g_aux = g(:, p==0, 1:4);
+g_aux = squeeze(g_aux(:, 1, :));
+pbaspect([1 1 1]);
+
+% Green's function along a meridian
+s = figure(3);
+clf;
+hold all;
+
+t0 = t(p==0);
+plot(t0, abs(g_aux(:, 1)), 'b--', 'LineWidth', 1.2);
+plot(t0, abs(g_aux(:, 2)), 'r-.', 'LineWidth', 1.2);
+plot(t0, abs(g_aux(:, 3)), 'k:',  'LineWidth', 1.2);
+plot(t0, abs(g_aux(:, 4)), 'm-',  'LineWidth', 1.2);
+
+sp = 10;
+plot(t0(1:sp:end), abs(g_aux(1:sp:end, 1)), 'bv', 'LineWidth', 1.5);
+plot(t0(1:sp:end), abs(g_aux(1:sp:end, 2)), 'rs', 'LineWidth', 1.5);
+plot(t0(1:sp:end), abs(g_aux(1:sp:end, 3)), 'ko',  'LineWidth', 1.5);
+plot(t0(1:sp:end), abs(g_aux(1:sp:end, 4)), 'mp',  'LineWidth', 1.5);
+
+legend('1 m', '2 m', '3 m', '4 m');
+
+axis tight;
+pbaspect([1 .7 1]);
+grid;
+xlabel('\theta');
+ylabel('g');