Update wide vs narrow code

.gitignore

@@ -37,3 +37,5 @@ Manifest.toml
 *.fig
 *.tif
 *.bmp
+*.jpg
+*.jpeg

0 - example/lyh - 相变/捷变频相变/FAR_block_sparsity_phase1.m

@@ -0,0 +1,48 @@
+close all;
+clear all;
+clc;
+M = 4;
+N = 128;
+%block_sparsity = 1;
+tol = 1e-5;
+trial = 20;
+epi = 0.02;
+result = zeros(N,25);
+for col = 4:N
+    for block_sparsity = 10:18
+        success_count = 0;
+        for loop = 1:trial
+            FAR_model = zeros(N,M*N);
+            %Cn = randperm(M)-1
+            for n = 0 : N-1
+                Cn = floor(rand()*M);
+                for q = 0 : N-1
+                    for p = 0:M-1
+                        FAR_model(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi));
+                    end
+                end
+            end
+            col_choose = randperm(N,col);
+            FAR_model = FAR_model(col_choose,:);
+            sparse_signal = zeros(M,N);
+            block = randperm(N,block_sparsity);
+            sparse_signal(:,block) = exp(1i*2*pi*rand(M,block_sparsity));
+            y = FAR_model * sparse_signal(:);
+            cvx_begin
+            variable x(M,N) complex
+            norm21 = 0;
+            for i = 1:N
+                norm21 = norm21 + norm(x(:,i));
+            end
+            minimize(norm21)
+            subject to
+            FAR_model * x(:) == y
+            cvx_end
+            if norm(x(:)-sparse_signal(:))<tol
+                success_count = success_count+1;
+            end
+        end
+        result(col,block_sparsity) = success_count/trial;
+    end
+end
+save('FARblockepsilon2.mat');

0 - example/lyh - 相变/捷变频相变/FAR_sparsity_phase.m

@@ -0,0 +1,47 @@
+close all;
+clear all;
+clc;
+M = 4;
+N = 128;
+%block_sparsity = 1;
+tol = 1e-5;
+trial = 50;
+epi = 0.02;
+result = zeros(N,25);
+for col = 4:4:128
+    for block_sparsity = 1:25
+        success_count = 0;
+        for loop = 1:trial
+            FAR_model = zeros(N,M*N);
+            %Cn = randperm(M)-1
+            for n = 0 : N-1
+                Cn = floor(rand()*M);
+                for q = 0 : N-1
+                    for p = 0:M-1
+                        FAR_model(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi));
+                    end
+                end
+            end
+            
+            col_choose = randperm(N,col);
+            FAR_model = FAR_model(col_choose,:);
+            sparse_signal = zeros(M,N);
+
+            block = randperm(N,block_sparsity);
+            sparse_signal(:,block) = exp(1i*2*pi*rand(M,block_sparsity));
+            
+            y = FAR_model * sparse_signal(:);
+            cvx_begin
+            variable x(M*N) complex
+            minimize(norm(x,1))
+            subject to
+            FAR_model * x == y
+            cvx_end
+            if norm(x-sparse_signal(:))<tol
+                success_count = success_count+1;
+            end
+        end
+        result(col,block_sparsity) = success_count/trial;
+    end
+end
+save('FARepsilon.mat');

0 - example/lyh - 相变/捷变频相变/theoretic.m

@@ -5,4 +5,3 @@ function n = theoretic(m,s,d)
     g = diff(f,t);
     t1 = solve(g);
     n = s*(m+t1^2)+(d-s)*int((u-t1)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t1,inf);
-end

0 - example/lyh - 相变/高斯相变/gauss_phase_real.m

@@ -0,0 +1,27 @@
+N = 100;
+gauss_phase_res = zeros(100,100);
+for col = 1:100
+    %¾ØÕóÉú³É
+    for p =1:100
+        suc = 0;
+        for loop = 1:50
+            x1 = zeros(N,1);
+            q = randperm(N,p);
+            x1(q) = randn(p,1);
+            fai = randn(col,N);
+            b = fai*x1;
+            cvx_begin quiet
+                variable x(N)
+                minimize( norm( x, 1 ) )
+                subject to
+                    fai * x == b
+            cvx_end
+            %disp((norm(x-x1,1)))
+            if (norm(x-x1,1))<10e-5
+                 suc = suc+1;
+             end
+        end
+        gauss_phase_res(col,p)=suc/50;
+    end
+end
+save gauss_phase_real;

1 - experiment/Circulant Matrix/circulant_matrix.m

@@ -4,22 +4,24 @@ clear;
 
 N = 100;
 eps = 1e-10;
+
+% 生成傅里叶矩阵
 F_n = dftmtx(N);
 
+% 输出傅里叶矩阵
 for i = 1: N
     for j = 1: N
         fprintf(string(F_n(i, j)) + ", ");
     end
     fprintf("\b\b\n");
 end
-
 fprintf("\n");
 
-% 检查傅里叶矩阵的正确性
+% 检查傅里叶矩阵
 assert(abs(F_n(2, 2) ^ 2 - F_n(2, 3)) < eps);
 assert(abs(F_n(2, 2) ^ N - 1) < eps);
 
-%% 检查循环卷积
+%% 循环卷积
 x = randn(N, 1);
 z = randn(N, 1);
 p = cconv(x, z, N);
@@ -33,11 +35,11 @@ for i = 1: N
     assert(abs(tmp - p(i)) < eps);
 end
 
-% 检查 x * z = A(x)z 的正确性
+% 检查 x * z = A(x)z
 A_x = toeplitz([x(1) fliplr(x(2:end)')], x);
 assert(is_vector_equal(p, A_x' * z, eps));
 
-% 检查卷积定理的正确性
+% 检查卷积定理
 x_hat = fft(x);
 z_hat = fft(z);
 xz_hat = fft(p);

Block RIP/Phase transitions/block_phase_transition.m

@@ -1,39 +1,56 @@
-function pt = block_phase_transition(N, M)
+function [N_b, N_s] = block_phase_transition(N, M)
     tau_min = 0;
     tau_max = 10;
     tau_interval = 0.05;
     tau_range = tau_min:tau_interval:tau_max;
     
     K_min = 1;
-    K_max = 25;
+    K_max = 15;
     K_range = K_min:K_max;
     
-    pt = zeros(25, 1);
+    N_b = zeros(length(K_range), 1);
+    N_s = zeros(length(K_range), 1);
     
-    cache_filename = "I.mat";
+    cache_filename_1 = "I_2.mat";
+    cache_filename_2 = "I_" + string(2 * M) + ".mat";
     
-    if exist(cache_filename, "file")
-        load(cache_filename);
+    if exist(cache_filename_1, "file")
+        load(cache_filename_1, "I_2");
     else
-        I = zeros(length(tau_range), 0);
+        I_2 = zeros(length(tau_range), 0);
     
         for tau_idx = 1:length(tau_range)
             tau = tau_range(tau_idx);
-            I(tau_idx) = calc_block_integral(tau, 4);
+            I_2(tau_idx) = calc_block_integral(tau, 2);
         end
-        save I.mat
+        save(cache_filename_1, "I_2");
+    end
+    
+    if exist(cache_filename_2, "file")
+        load(cache_filename_2, "I_2M");
+    else
+        I_2M = zeros(length(tau_range), 0);
+    
+        parfor tau_idx = 1:length(tau_range)
+            tau = tau_range(tau_idx);
+            I_2M(tau_idx) = calc_block_integral(tau, 2 * M);
+        end
+        save(cache_filename_2, "I_2M");
     end
     
     
     for K_idx = 1:length(K_range)
         K = K_range(K_idx);
-        f_set = zeros(length(tau_range), 1);
+        f_set_1 = zeros(length(tau_range), 1);
+        f_set_2 = zeros(length(tau_range), 1);
     
         for tau_idx = 1:length(tau_range)
             tau = tau_range(tau_idx);
-            f_set(tau_idx) = 1/2 * (K * (2 * M + tau ^ 2) + (N - K) * I(tau_idx));
+            f_set_1(tau_idx) = 1/2 * (K * (2 * M + tau ^ 2) + (N - K) * I_2M(tau_idx));
+            f_set_2(tau_idx) = M/2 * (K * (2 + tau ^ 2) + (N - K) * I_2(tau_idx));
         end
     
-        pt(K_idx) = min(f_set);
+        N_b(K_idx) = min(f_set_1);
+        N_s(K_idx) = min(f_set_2);
     end
 end

Block RIP/Phase transitions/block_phase_transition_0.m

@@ -1,40 +0,0 @@
-function pt = block_phase_transition_0()
-    tau_min = 0;
-    tau_max = 100;
-    tau_interval = 0.1;
-    tau_range = tau_min:tau_interval:tau_max;
-    
-    s_b_min = 1;
-    s_b_max = 25;
-    s_b_range = s_b_min:s_b_max;
-    
-    pt = zeros(25, 1);
-    
-    cache_filename = "I.mat";
-    
-    if exist(cache_filename, "file")
-        load(cache_filename);
-    else
-        I = zeros(length(tau_range), 0);
-    
-        for tau_idx = 1:length(tau_range)
-            tau = tau_range(tau_idx);
-            I(tau_idx) = calc_block_integral(tau, 4);
-        end
-    
-    end
-    
-    save I.mat
-    
-    for s_b_idx = 1:length(s_b_range)
-        s_b = s_b_range(s_b_idx);
-        f_set = zeros(length(tau_range), 1);
-    
-        for tau_idx = 1:length(tau_range)
-            tau = tau_range(tau_idx);
-            f_set(tau_idx) = s_b * (1 + tau ^ 2) + (100 - s_b) * I(tau_idx);
-        end
-    
-        pt(s_b_idx) = min(f_set);
-    end
-end

Block RIP/Phase transitions/block_phase_transition_2.m

@@ -1,20 +0,0 @@
-function pt = block_phase_transition_2(N, M)
-    K_min = 1;
-    K_max = 25;
-    K_range = K_min:K_max;
-    pt = zeros(25, 1);
-    d = N;
-    m = M;
-    
-    for K_idx = 1:length(K_range)
-        K = K_range(K_idx);
-        s = K;
-        syms t;
-        syms u;
-        f = s*(m+t^2)+(d-s)*int((u-t)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t,inf);
-        g = diff(f,t);
-        t1 = solve(g);
-        n = s*(m+t1^2)+(d-s)*int((u-t1)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t1,inf);
-        pt(K_idx) = n;
-    end
-end

Block RIP/Phase transitions/block_phase_transition_formula.m

@@ -0,0 +1,22 @@
+function pt = block_phase_transition_formula(N, M)
+    K_min = 1;
+    K_max = 25;
+    K_range = K_min:K_max;
+    pt = zeros(25, 1);
+    d = N;
+    m = M;
+    
+    for K_idx = 1:length(K_range)
+        s = K_range(K_idx);        
+        pt(K_idx) = theoretic(m, s, d);
+    end
+end
+
+function n = theoretic(m,s,d)
+    syms t;
+    syms u;
+    f = s*(m+t^2)+(d-s)*int((u-t)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t,inf);
+    g = diff(f,t);
+    t1 = solve(g);
+    n = s*(m+t1^2)+(d-s)*int((u-t1)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t1,inf);
+end

Block RIP/Phase transitions/draw.m

@@ -29,7 +29,8 @@ ax = gca;
 grid(ax, 'on');
 set(ax, 'Visible', 'on');
 
-pt = block_phase_transition_2(128, 4);
+% pt = block_phase_transition_2(128, 4);
+pt = block_phase_transition_2(16*16, 16);
 l = line(1:size(pt), pt);
 l.Color = "w";
 l.LineWidth = 5;

Block RIP/Phase transitions/draw_ft_curve.m

@@ -0,0 +1,13 @@
+M = 10;
+N = 320 + 160;
+epi = 0.02;
+[pt_block, pt_sparse] = block_phase_transition(N, M);
+plot(pt_sparse);
+% hold on;
+% plot(pt_sparse);
+% xlim([1, 5]);
+ylim([1, N]);
+xlabel("Sparsity (k)");
+ylabel("Number of measurements (n)")
+legend("Block sparse recovery", "Sparse recovery");
+title("Phase Transition Curve (N = "+string(N)+", M = "+string(M)+")");

FAR vs PD/FAR_simu.m

@@ -1,4 +1,4 @@
-function P_fa = FAR_simu(threshold)
+function threshold = FAR_simu(P_fa_FAR)
     global N M epi trail_times
 
     % Define measurement matrix
@@ -9,64 +9,124 @@ function P_fa = FAR_simu(threshold)
     Lambda_C = setdiff(1:N*M, Lambda);
 
     % Try "trail_times" times
-    NN = N;
-    results = zeros(trail_times, N * M);
+    results_0 = zeros(trail_times, N * M);
+    results_1 = zeros(trail_times, N * M);
+    thresholds_0 = [];
+    thresholds_1 = [];
     
-    if exist("FAR_recovery_1_LASSO_5.mat")
-        load("FAR_recovery_0_debiasedLASSO_500.mat", "results");
+    noise_sigma = 0.01;
+    noise_sigma_2 = noise_sigma ^ 2;
+    
+    if exist("520BP.mat")
+        load("50.mat", "results_0", "results_1");
     else
+        % 0 假设
         for T = 1: trail_times
             % y = Ax + n
-            n = randn(NN, 1) * 0.1;
+            n = randn(N, 1) * noise_sigma;
+            y_noise = n;
+            [x_hat, threshold] = debiased_LASSO(A, y_noise, P_fa_FAR, noise_sigma_2);
+            results_0(T, :) = x_hat;
+            thresholds_0 = [thresholds_0 threshold];
+        end
+        % 1 假设
+        for T = 1: trail_times
+            % y = Ax + n
+            n = randn(N, 1) * noise_sigma;
             y_noise = A * x + n;
-            
-            % x_hat = CS(A, y)
-            x_hat = recovery(A, y_noise);
-            results(T, :) = x_hat;
+            [x_hat, threshold] = debiased_LASSO(A, y_noise, P_fa_FAR, noise_sigma_2);
+            results_1(T, :) = x_hat;
+            thresholds_1 = [thresholds_1 threshold];
         end
     end
 
     % Distribute of H_0 and H_1
-    H_0_distribute = zeros(1, length(Lambda_C));
-    H_1_distribute = zeros(1, length(Lambda));
+    H_00_distribute = zeros(1, length(Lambda_C));
+    H_01_distribute = zeros(1, length(Lambda));
+    H_10_distribute = zeros(1, length(Lambda_C));
+    H_11_distribute = zeros(1, length(Lambda));
 
     i1 = 1;
     i2 = 1;
-    Ts = [];
+    T_00 = [];
+    T_01 = [];
+    T_10 = [];
+    T_11 = [];
+    
+    for threshold = 1: trail_times
+        x_0_hat = results_0(threshold, :);
+        x_1_hat = results_1(threshold, :);
+        t_00 = 0;
+        t_01 = 0;
+        t_10 = 0;
+        t_11 = 0;
 
-    for t = 1: trail_times
-        x_hat = results(t, :);
-        T = 0;
         for i = 1: length(x)
             if ismember(i, Lambda_C)
-                H_0_distribute(i1) = x_hat(i);
+                t_00 = t_00 + abs(x_0_hat(i));
+                t_10 = t_10 + abs(x_1_hat(i));
+                H_00_distribute(i1) = x_0_hat(i);
+                H_10_distribute(i1) = x_1_hat(i);
                 i1 = i1 + 1;
             elseif ismember(i, Lambda)
-                T = T + abs(x_hat(i));
-                H_1_distribute(i2) = x_hat(i);
+                t_01 = t_01 + abs(x_0_hat(i));
+                t_11 = t_11 + abs(x_1_hat(i));
+                H_01_distribute(i2) = x_0_hat(i);
+                H_11_distribute(i2) = x_1_hat(i);
                 i2 = i2 + 1;
             end
         end
-        Ts = [Ts T];
+        T_00 = [T_00 t_00];
+        T_01 = [T_01 t_01];
+        T_10 = [T_10 t_10];
+        T_11 = [T_11 t_11];
     end
 
 
     % Draw
     figure(1)
-    subplot(2, 1, 1);
+    subplot(2, 2, 1);
     title("Freq histogram of H_0");
-    histfit(real(H_0_distribute));
+    histfit(real(H_00_distribute));
     
-    subplot(2, 1, 2);
+    subplot(2, 2, 2);
     title("Freq histogram of H_1");
-    histfit(real(H_1_distribute));
+    histfit(real(H_01_distribute));
 
-    H_0_mean = mean(H_0_distribute);
-    H_0_std = std(H_0_distribute);
-    H_1_mean = mean(H_1_distribute);
-    H_1_std = std(H_1_distribute);
+    subplot(2, 2, 3);
+    title("Freq histogram of H_0");
+    histfit(real(H_10_distribute));
+    
+    subplot(2, 2, 4);
+    title("Freq histogram of H_1");
+    histfit(real(H_11_distribute));
+    
+    H_00_mean = mean(H_00_distribute);
+    H_00_std = std(H_00_distribute);
+    H_01_mean = mean(H_01_distribute);
+    H_01_std = std(H_01_distribute);
+    H_10_mean = mean(H_10_distribute);
+    H_10_std = std(H_10_distribute);
+    H_11_mean = mean(H_11_distribute);
+    H_11_std = std(H_11_distribute);
+
+    fprintf("H_00: mu = %f, std = %f\n", H_00_mean, H_00_std);
+    fprintf("H_01: mu = %f, std = %f\n", H_01_mean, H_01_std);
+    fprintf("H_10: mu = %f, std = %f\n", H_10_mean, H_10_std);
+    fprintf("H_11: mu = %f, std = %f\n", H_11_mean, H_11_std);
+
+    figure(2);
+    subplot(2, 2, 1);
+    histfit(real(T_00));
+    
+    subplot(2, 2, 2);
+    histfit(real(T_01));
+
+    subplot(2, 2, 3);
+    histfit(real(T_10));
+    
+    subplot(2, 2, 4);
+    histfit(real(T_11));
 
-    fprintf("H_0: mu = %f, std = %f\n", H_0_mean, H_0_std);
-    fprintf("H_1: mu = %f, std = %f\n", H_1_mean, H_1_std);
     % P_fa = normcdf(threshold, H_0_mean, H_0_std);
 end

FAR vs PD/Main.m

@@ -7,19 +7,20 @@ N = 64;
 M = 4;
 K = 10;
 lambda = zeros(M * N, 1);
-lambda(:) = 0.3;
-tau = 1e-4;
-iter_max = 500;
-% method = "debiased_LASSO";
+lambda(:) = 0.8;
+tau = 1e-6;
+iter_max = 50;
+method = "debiased_LASSO";
 % method = "LASSO";
-method = "BP";
+% method = "BP";
 % method = "VAMP";
-trail_times = 1e2;
+trail_times = 1000;
 epi = 0;
 extend_target = 1;
 
 %% FAR
-P_fa_FAR = FAR_simu(3);
+P_fa_FAR = 1e-5;
+threshold_FAR = FAR_simu(P_fa_FAR);
 % fprintf("FAR: %f\n", P_fa_FAR);
 
 %% PD

FAR vs PD/check_orth.m

@@ -0,0 +1,25 @@
+clc; clear;
+N = 16;
+M = 4;
+epi = 0;
+
+d_n = floor(rand(N, 1)*M) / M;
+A = get_Psi(N, M, d_n, epi);
+B = get_Psi_2(N, M, d_n, epi);
+C = A - B;
+real_C = real(C);
+imag_C = imag(C);
+
+% subplot(121)
+% heatmap(real_C)
+% 
+% subplot(122)
+% heatmap(imag_C)
+
+r = zeros(N);
+for i = 1: N
+    for j = 1: N
+        r(i, j) = A(i, :) * (A(j, :)');
+    end
+end
+heatmap(abs(r));

FAR vs PD/config.json

@@ -1,11 +0,0 @@
-{
-    "N": 64,
-    "M": 4,
-    "K": 10,
-    "method": "debiased_LASSO",
-    "trail_times": 5e2,
-    "extend_target": true,
-    "debiased_LASSO_params": {
-
-    }
-}

FAR vs PD/get_Psi.m

@@ -1,11 +0,0 @@
-function Psi = get_Psi(N, M, epi)
-    Psi = zeros(N,M*N);
-    for n = 0 : N-1
-        Cn = floor(rand()*M);
-        for q = 0 : N-1
-            for p = 0:M-1
-                Psi(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi)) / sqrt(N);
-            end
-        end
-    end
-end

FAR vs PD/recovery.m

@@ -25,6 +25,14 @@ function x_hat = recovery(A, y_noise)
         x_d_CROD = x_LASSO + A'*(y_noise - A*x_LASSO)/Q_hat;
         x_hat = x_d_CROD;
 
+        RSS = 1/sz(2) * norm(y_noise - A * x_LASSO, 2)^2;
+        noise_sigma_2 = 1e-4;
+        P_fa = normcdf(1, 0, 1);
+        sigma_w_2 = (gamma * (1-gamma)) / ((gamma - Rho)^2) * RSS + noise_sigma_2;
+        k_d = -sigma_w_2 * log(P_fa);
+
+        fprintf("%f %f %f\n", gamma, sigma_w_2, k_d);
+
     elseif method == "LASSO"
         sz = size(A);
         N = sz(2);
@@ -47,7 +55,7 @@ function x_hat = recovery(A, y_noise)
         cvx_end
     else
         global lambda tau iter_max;
-        [x_hat, z_hat_d] = cVAMPro(y_noise, A, lambda, tau, iter_max);
+        [x_hat, x_hat_d] = cVAMPro(y_noise, A, lambda, tau, iter_max);
     end
 end
 

FAR vs PD/untitled.m

@@ -1,12 +0,0 @@
-clc;
-clear;
-N = 100;
-F = dftmtx(N);
-F_inv = inv(F);
-
-result = zeros(N);
-for i = 1: N
-    for j = 1: N
-        result(i, j) = F_inv(i,:) * conj(F_inv(j,:))';
-    end
-end

Random PRI/Main.m

@@ -0,0 +1,41 @@
+clc; clear; close all;
+
+parameters;
+
+global trail_times;
+
+[s_T_narrow, A] = narrow_signal_model(t, 1);
+
+A = A / 10;
+
+AH = A';
+
+H0 = [];
+for T = 1: trail_times
+    x = zeros(10, 1);
+    % x(4) = 0.8;
+    y = A * x;
+    y_noise = awgn(y, 15);
+    
+    x_hat = AH * y_noise;
+
+    H0 = [H0 x_hat'];
+end
+
+H1 = [];
+for T = 1: trail_times
+    x = zeros(10, 1);
+    x(4) = 0.8;
+    y = A * x;
+    y_noise = awgn(y, 15);
+    
+    x_hat = AH * y_noise;
+
+    H1 = [H1 x_hat(4)];
+end
+
+% figure; histfit(real(H0));
+% figure; histfit(real(H1));
+
+figure;
+plot()

Random PRI/narrow_signal_model.m

@@ -0,0 +1,25 @@
+function [s_T_narrow, A_narrow] = narrow_signal_model(t, betas_narrow)
+    global f_c N_narrow N FIGURE N_wide numP f_s K_chirp T_r;
+
+    tt = 0: 1/f_s: numP * T_r - 1/f_s;
+
+    % 使用等效散射系数生成窄带情况信号模型
+    s_T_narrow = zeros(1, N * numP);
+    for j = 1: numP
+        idx = ((j-1) * N + (1:N_narrow));
+        tp = idx / f_s;
+        s_T_narrow(1, idx) = exp(1j * 2 * pi * (f_c * tp + 0.5 * K_chirp * tp .^ 2));
+    end
+
+    % s_T_narrow = exp(1j .* 2 .* pi .* f_c .* t);
+    s_R_narrows = zeros(N_narrow, N * numP);
+    for i = 1: N_narrow
+        p = i * (N_wide / N_narrow);
+        s_R_narrows(i, :) = betas_narrow * [zeros(1, p), s_T_narrow(1: end-p)];
+    end
+    
+    A_narrow = s_R_narrows';
+    if FIGURE
+        subplot(2, 1, 1); plot(tt, real(s_T_narrow));
+        subplot(2, 1, 2); plot(tt, real(s_R_narrows(1, :)));
+    end

Random PRI/parameters.m

@@ -0,0 +1,59 @@
+global bandwidth T_p f_s T_s K_chirp f_c PRF T_r numP t tP c N M R_0...
+    delta_R_wide delta_R_narrow N_wide N_narrow ranges_wide...
+    ranges_narrow betas_wide FIGURE K_wide K_narrow trail_times...
+    method N_high TEST lambda tau iter_max noise_sigma FAR_N FAR_M noise_sigma_2 DEBUG_LEVEL...
+;
+
+f_c = 10e9; % X wave
+bandwidth = 100e6;
+T_p = 1e-5;
+f_s = 2 * bandwidth;
+T_s = 1 / f_s;
+K_chirp = bandwidth / T_p;
+T_r = T_p * 10;
+PRF = 1 / T_r;
+
+numP = 10;
+
+t = 0: 1 / f_s: T_r - 1 / f_s;
+tP = 0: 1 / f_s: T_r * numP - 1 / f_s;
+
+c = 3e8;        % 光速
+
+N = length(t);
+N_high = T_p * f_s;
+M = round(T_p * bandwidth);
+R_0 = 0;
+delta_R_wide = c ./ 2 ./ bandwidth;
+delta_R_narrow = T_p .* c ./ 2;
+N_wide = 100;
+N_narrow = round(N_wide / M);
+
+% alert(N > N_narrow);
+
+ranges_wide = R_0 + (0:N_wide) * delta_R_wide; % [1000, 4000]
+ranges_narrow = R_0 + (0:N_narrow) * delta_R_narrow; % [1000, 4000]
+betas_wide = linspace(1, 0.1, N_wide) + 1j * linspace(0.1, 1, N_wide);
+
+K_wide = round(N_wide / 10);
+K_narrow = round(K_wide / 10);
+
+FIGURE = false;
+trail_times = 5000;
+% method = "debiased_LASSO";
+% method = "LASSO";
+% method = "BP";
+method = "cVAMPro";
+TEST = false;
+
+lambda = 0.01;
+tau = 1e-6;
+iter_max = 1000;
+
+noise_sigma = 0.1;
+noise_sigma_2 = 0.1;
+
+FAR_N = 480;
+FAR_M = 10;
+
+DEBUG_LEVEL = 1;

basic/FISTA.m

@@ -0,0 +1,27 @@
+function [z] = FISTA(y, A, lambda, delta)
+
+x_pre = A'*y;
+t = 1;
+z = x_pre;
+z_pre = z;
+t_pre = t;
+N = size(A, 2);
+diff = 1;
+E = eig(A'*A);
+L = E(end);
+temp1 = A'*y/L;
+temp2 = eye(N) - A'*A/L;
+k = 0;
+while((diff > delta) && (k < 1000))
+    temp = temp1 + temp2 * z_pre;
+    x = sft_thd(temp, lambda/L);
+    t = 0.5*(1 + sqrt(1+4*t_pre*t_pre));
+    z = x + (x - x_pre) * (t_pre-1) / t;
+    diff = mean(abs(z_pre - z)); 
+    x_pre = x;
+    z_pre = z;
+    t_pre = t;
+    k = k + 1;
+end
+
+end

basic/cVAMPro.m

@@ -38,7 +38,6 @@ function [x_hat_wl, x_hat_d] = cVAMPro(y, A, lambda, tau, Kit)
         h_1 = h_1_next;
         Q_1 = Q_1_next;
     end
-
 end
 
 % SoftThreshold function
@@ -47,8 +46,12 @@ function x = ST(h_1, lambda, Q_1)
     x = zeros(N, M);
 
     for i = 1:N
-        sign = h_1(i) ./ abs(h_1(i));
-        diff = abs(h_1(i)) - lambda(i);
+        if h_1(i) == 0
+            sign = 1;
+        else
+            sign = h_1(i) ./ abs(h_1(i));
+        end
+        diff = abs(h_1(i)) - lambda;
         x(i) = sign .* (diff ./ Q_1) .* SF(diff);
     end
 
@@ -73,9 +76,10 @@ function chi_1 = F1(x_1, lambda, Q_1)
     count = 0;
 
     for i = 1:N
-        temp = Q_1 * abs(x_1(i)) + lambda(i);
-        count = count + (2 - lambda(i) / temp) * SF(abs(x_1(i)));
-        %         count = count + (2-lambda(i)/temp) * (abs(x_1(i)) > 1e-4);
+        % temp = Q_1 * abs(x_1(i)) + lambda(i);
+        % count = count + (2 - lambda(i) / temp) * SF(abs(x_1(i)));
+        temp = Q_1 * abs(x_1(i)) + lambda;
+        count = count + (2 - lambda / temp) * SF(abs(x_1(i)));
     end
 
     chi_1 = count / (2 * N * Q_1);

basic/debiased_LASSO.m

@@ -0,0 +1,39 @@
+function [x_hat, sigma_w_2, threshold] = debiased_LASSO(A, y_noise, P_fa, noise_sigma_2, LASSO_lambda)
+
+    if nargin < 5
+        LASSO_lambda = 1e-1;
+    end
+    
+    sz = size(A);
+    N = sz(2);
+
+    gamma = sz(1) / sz(2);
+
+    cvx_begin quiet
+        variable x_LASSO(N) complex
+        minimize(LASSO_lambda * norm(x_LASSO, 1) + norm(y_noise - A * x_LASSO, 2))
+    cvx_end
+    
+    rho_active = sum(abs(x_LASSO) > 1e-3) / N;
+    Q_hat = (gamma - rho_active) / (1 - rho_active);
+    Rho = sum((abs(x_LASSO) > 1e-3) .* (2 - LASSO_lambda ./ (Q_hat * abs(x_LASSO) + LASSO_lambda))) / 2 / N;
+    diff = 1;
+
+    while (diff > 1e-4)
+        Rho_pre = Rho;
+        Rho = sum((abs(x_LASSO) > 1e-3) .* (2 - LASSO_lambda ./ ((gamma - Rho) / (1 - Rho) * abs(x_LASSO) + LASSO_lambda))) / 2 / N;
+        diff = abs(Rho - Rho_pre);
+    end
+
+    Q_hat = (gamma - Rho) / (1 - Rho);
+    x_d_CROD = x_LASSO + A' * (y_noise - A * x_LASSO) / Q_hat;
+    x_hat = x_d_CROD;
+
+    RSS = 1 / sz(1) * norm(y_noise - A * x_LASSO, 2) ^ 2;
+    sigma_w_2 = (gamma * (1 - gamma)) / ((gamma - Rho) ^ 2) * RSS + noise_sigma_2;
+    threshold = -sigma_w_2 * log(P_fa);
+    xx = x_LASSO;
+    FAR_M = 1 / gamma;
+    xx(FAR_M+1: 2*FAR_M) = xx(FAR_M+1: 2*FAR_M) - 1;
+    xx = real(xx);
+end

basic/debiased_LASSO_FISTA.m

@@ -0,0 +1,35 @@
+function [x_hat, sigma_w_2, threshold] = debiased_LASSO_FISTA(A, y_noise, P_fa, noise_sigma_2, LASSO_lambda)
+
+    if nargin < 5
+        LASSO_lambda = 1e-1;
+    end
+    
+    sz = size(A);
+    N = sz(2);
+
+    gamma = sz(1) / sz(2);
+
+    x_LASSO = FISTA(y_noise, A, LASSO_lambda, 1e-5);
+    
+    rho_active = sum(abs(x_LASSO) > 1e-3) / N;
+    Q_hat = (gamma - rho_active) / (1 - rho_active);
+    Rho = sum((abs(x_LASSO) > 1e-3) .* (2 - LASSO_lambda ./ (Q_hat * abs(x_LASSO) + LASSO_lambda))) / 2 / N;
+    diff = 1;
+
+    while (diff > 1e-4)
+        Rho_pre = Rho;
+        Rho = sum((abs(x_LASSO) > 1e-3) .* (2 - LASSO_lambda ./ ((gamma - Rho) / (1 - Rho) * abs(x_LASSO) + LASSO_lambda))) / 2 / N;
+        diff = abs(Rho - Rho_pre);
+    end
+
+    Q_hat = (gamma - Rho) / (1 - Rho);
+    x_d_CROD = x_LASSO + A' * (y_noise - A * x_LASSO) / Q_hat;
+    x_hat = x_d_CROD;
+
+    RSS = 1 / sz(1) * norm(y_noise - A * x_LASSO, 2) ^ 2;
+    sigma_w_2 = (gamma * (1 - gamma)) / ((gamma - Rho) ^ 2) * RSS + noise_sigma_2;
+    threshold = -sigma_w_2 * log(P_fa);
+    xx = x_LASSO;
+    FAR_M = 1 / gamma;
+    mean_xx = mean(xx(FAR_M+1: 2*FAR_M));
+end

basic/get_Psi.m

@@ -0,0 +1,31 @@
+function [C_n, Psi] = get_Psi(N, M, epi)
+    if nargin == 2
+        epi = 0;
+    end
+    k = 1 / sqrt(M * N);
+    C_n = zeros(M, 1);
+
+    % Determined code 1
+    % C_n = repelem(0:M-1, N/M);
+
+    % Determined code 2
+    % C_n = repmat(0:M-1, 1, N/M);
+
+    % Determined code 3
+    C_n = repmat(0:M-1, 1, N/M);
+    C_n = C_n(randperm(length(C_n)));
+
+    Psi = zeros(N,M*N);
+    for n = 0 : N-1
+        col = zeros(1, M*N);
+        % C_n(n+1) = floor(rand()*M);
+        for q = 0 : N-1
+            for p = 0:M-1
+                f1 = p/M*C_n(n+1);
+                f2 = q/N*n*(1+C_n(n+1)*epi);
+                col(q*M+p+1) = exp(1j*2*pi*(f1 + f2));
+            end
+        end
+        Psi(n+1, :) = col * k;
+    end
+end

basic/get_Psi_KR_product.m

@@ -0,0 +1,29 @@
+function Psi = get_Psi_KR_product(N, M, epi)
+    d_n = floor(rand(N, 1)*M) / M;
+    zeta_n = 1 + d_n .* epi; % epi = B / f_c
+
+    R = zeros(N, M);
+    D = zeros(N, N);
+    for n = 0: N-1
+        for m = 0: M-1
+            R(n+1, m+1) = exp(-1j * 2 * pi * m * d_n(n+1));
+        end
+    end
+
+    for n = 0: N-1
+        for l = 0: N-1
+            D(n+1, l+1) = exp(-1j * 2 * pi * l * n * zeta_n(n+1) / N);
+        end
+    end
+
+    Psi = KR_product(R', D')';
+end
+
+function kr = KR_product(F, G)
+    nR_F = size(F, 1);
+    nR_G = size(G, 1);
+    mul = ones(nR_G, 1);
+    FF = kron(F, mul);
+    GG = repmat(G, nR_F, 1);
+    kr = FF .* GG;
+end

basic/get_Psi_sample.m

@@ -0,0 +1,89 @@
+function [s_T, A, freqs] = FAR_signal_model_sample(FAR_N, FAR_M)
+    global c f_c N N_wide B ranges_wide betas_wide FIGURE numP f_s N_high T_r T_p TEST method;
+
+    pulse_number = FAR_N;
+    single_pulse_length = N;
+    all_pulse_length = single_pulse_length * pulse_number;
+    pulse_width = N_high;
+    
+    delta_f = B / FAR_M;
+    delta_r = c / (2 * B);
+    delta_v = c / (2 * f_c * FAR_N * T_r);
+
+    ranges = (0: FAR_M - 1) * delta_r;
+    velos = (0: FAR_N - 1) * delta_v;
+    delays = 2 * ranges / c;
+
+    t_single = (0: 1/f_s: T_p-1/f_s)';
+    t = (0: 1/f_s: FAR_N * T_r-1/f_s)';
+    s_T = zeros(length(t), 1);
+    T_p_N = floor(T_p * f_s);
+    T_r_N = floor(T_r * f_s);
+
+    freqs = zeros(FAR_N, 1);
+    for i = 1: FAR_N
+        C_n = floor(rand() * FAR_M);
+        freqs(i) = f_c + C_n * delta_f;
+        % freqs(i) = f_c + i * delta_f;
+    end
+    
+    for i = 1: FAR_N
+        N_l = (i - 1) * T_r_N + 1;
+        N_r = (i - 1) * T_r_N + T_p_N;
+        s_T(N_l:N_r) = exp(1j * 2 * pi * freqs(i) * t_single);
+    end
+
+    % echo
+    s_R = zeros(length(s_T), FAR_M, FAR_N);
+    s_D = zeros(length(s_T), FAR_M, FAR_N);
+    for range_idx = 1: FAR_M % R
+        r = ranges(range_idx);
+        for doppler_idx = 1: FAR_N % D
+            v = velos(doppler_idx);
+            f_n = freqs(doppler_idx);
+            for pulse_idx = 1: FAR_N % N
+                nT_r = (pulse_idx - 1) * T_r;
+                rr = r + v * nT_r;
+                delay = 2 * rr / c;            
+                t_l = delay + nT_r;
+                N_l = max(1, round(t_l * f_s) + 1);
+                N_r = N_l + T_p_N - 1;
+    
+                s_D(N_l: N_r, range_idx, doppler_idx) = exp(-1j * 2 * pi * f_n * (2 / c) * (r + v * (nT_r + t_single)));
+            end
+        end
+    end
+
+    % Sampling
+    A = zeros(FAR_N, FAR_N * FAR_M);
+    for range_idx = 1: FAR_M % R
+        for doppler_idx = 1: FAR_N % D
+            index = (doppler_idx - 1) * FAR_M + range_idx;
+            for pulse_idx = 1: FAR_N % N
+                t_sample = delays(range_idx) + T_p + (pulse_idx - 1) * T_r;
+                A(pulse_idx, index) = s_D(round(t_sample * f_s), range_idx, doppler_idx);
+            end
+        end
+    end
+    A = get_Psi(FAR_N, FAR_M, 0);
+    AAH = A * A';
+    A = A ./ sqrt(AAH(1, 1));
+
+    if 1 == 0
+        Lambda = 1:10;
+        x_wide = zeros(N_wide, 1);
+        x_wide(Lambda) = 1;
+        n_wide = randn(length(s_T_FAR), 1) * 0.01;
+        y = A_wide * x_wide;
+        y_noise = y + n_wide;
+        x_hat = recovery(A_wide, y_noise, method);
+
+        figure(1);
+        subplot(211); plot(abs(x_wide));
+        subplot(212); plot(abs(x_hat));
+
+        figure(2);
+        subplot(311); plot(real(A_wide * x_wide));
+        subplot(312); plot(real(A_wide * x_hat));
+        subplot(313); plot(real(A_wide * (x_wide - x_hat)));
+    end

basic/get_noise.m

@@ -0,0 +1,11 @@
+function noise = get_noise(noise_sigma, N, M, type)
+    if nargin < 4
+        type = "complex";
+    end
+    
+    if type == "complex"
+        noise = random('Normal', 0, noise_sigma / sqrt(2), N, M) + 1j * random('Normal', 0, noise_sigma / sqrt(2), N, M);
+    else
+        noise = random('Normal', 0, noise_sigma, N, M);
+    end
+end

basic/narrow_signal_model_sample.m

@@ -0,0 +1,87 @@
+function [s_T, A] = narrow_signal_model_sample(FAR_N)
+    FAR_M = 1;
+    global c f_c N N_wide B ranges_wide betas_wide FIGURE numP f_s N_high T_r T_p TEST method;
+
+    pulse_number = FAR_N;
+    single_pulse_length = N;
+    all_pulse_length = single_pulse_length * pulse_number;
+    pulse_width = N_high;
+    
+    delta_f = 0;
+    delta_r = T_p * c / 2;
+    delta_v = c / (2 * f_c * FAR_N * T_r);
+
+    ranges = (0: FAR_M - 1) * delta_r;
+    velos = (0: FAR_N - 1) * delta_v;
+    delays = 2 * ranges / c;
+
+    t_single = (0: 1/f_s: T_p-1/f_s)';
+    t = (0: 1/f_s: FAR_N * T_r-1/f_s)';
+    s_T = zeros(length(t), 1);
+    T_p_N = floor(T_p * f_s);
+    T_r_N = floor(T_r * f_s);
+
+    freqs = zeros(FAR_N, 1);
+    for i = 1: FAR_N
+        C_n = floor(rand() * FAR_M);
+        freqs(i) = f_c + C_n * delta_f;
+        % freqs(i) = f_c + i * delta_f;
+    end
+    
+    for i = 1: FAR_N
+        N_l = (i - 1) * T_r_N + 1;
+        N_r = (i - 1) * T_r_N + T_p_N;
+        s_T(N_l:N_r) = exp(1j * 2 * pi * freqs(i) * t_single);
+    end
+
+    % echo
+    s_R = zeros(length(s_T), FAR_M, FAR_N);
+    s_D = zeros(length(s_T), FAR_M, FAR_N);
+    for range_idx = 1: FAR_M % R
+        r = ranges(range_idx);
+        for doppler_idx = 1: FAR_N % D
+            v = velos(doppler_idx);
+            f_n = freqs(doppler_idx);
+            for pulse_idx = 1: FAR_N % N
+                nT_r = (pulse_idx - 1) * T_r;
+                rr = r + v * nT_r;
+                delay = 2 * rr / c;            
+                t_l = delay + nT_r;
+                N_l = max(1, round(t_l * f_s) + 1);
+                N_r = N_l + T_p_N - 1;
+    
+                s_D(N_l: N_r, range_idx, doppler_idx) = exp(-1j * 2 * pi * f_n * (2 / c) * (r + v * (nT_r + t_single)));
+            end
+        end
+    end
+
+    % Sampling
+    A = zeros(FAR_N, FAR_N * FAR_M);
+    for range_idx = 1: FAR_M % R
+        for doppler_idx = 1: FAR_N % D
+            index = (doppler_idx - 1) * FAR_M + range_idx;
+            for pulse_idx = 1: FAR_N % N
+                t_sample = delays(range_idx) + T_p + (pulse_idx - 1) * T_r;
+                A(pulse_idx, index) = s_D(round(t_sample * f_s), range_idx, doppler_idx);
+            end
+        end
+    end
+
+    if 1 == 0
+        Lambda = 1:10;
+        x_wide = zeros(N_wide, 1);
+        x_wide(Lambda) = 1;
+        n_wide = randn(length(s_T_FAR), 1) * 0.01;
+        y = A_wide * x_wide;
+        y_noise = y + n_wide;
+        x_hat = recovery(A_wide, y_noise, method);
+
+        figure(1);
+        subplot(211); plot(abs(x_wide));
+        subplot(212); plot(abs(x_hat));
+
+        figure(2);
+        subplot(311); plot(real(A_wide * x_wide));
+        subplot(312); plot(real(A_wide * x_hat));
+        subplot(313); plot(real(A_wide * (x_wide - x_hat)));
+    end

basic/sft_thd.m

@@ -0,0 +1,20 @@
+function y = sft_thd(x, thd)
+
+if isequal(size(x), size(thd))
+    
+    tmp = abs(x);
+    y = x;
+    y(tmp <= thd) = 0;
+    y(tmp > thd) = (tmp(tmp > thd) - thd(tmp > thd)) .* x(tmp > thd) ./ tmp(tmp > thd);
+    
+else
+    
+    tmp = abs(x);
+    y = x;
+    y(tmp <= thd) = 0;
+    y(tmp > thd) = (tmp(tmp > thd) - thd) .* x(tmp > thd) ./ tmp(tmp > thd);
+    
+end
+
+
+end

wide vs narrow/Main3.m

@@ -0,0 +1,253 @@
+clc; clear; close all;
+
+parameters;
+
+P_fas = [1e-7, 5e-7, 1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1e0];
+
+P_d_FARs = zeros(length(P_fas), 1);
+P_d_narrows = zeros(length(P_fas), 1);
+
+
+%% FAR
+
+% Signal Model
+% epi = B / f_c;
+epi = 0;
+[C_n, A] = get_Psi(FAR_N, FAR_M, epi);
+fac = A * A';
+A = A / sqrt(abs(fac(1, 1)));
+f_n = f_c + C_n * B / M;
+
+[Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true);
+
+% Define filename
+filename = ...
+    "data2/" + ...
+    "FAR" + "_" + ...
+    string(FAR_N) + "_" + ...
+    string(FAR_M)+ "_" + ...
+    method + "_" + ...
+    string(noise_sigma) + "_" + ...
+    string(length(Lambda)) + "_" + ...
+    string(LASSO_lambda) + "_" + ...
+    trail_times + ...
+    ".mat";
+
+
+for P_fa_idx = 1: length(P_fas)
+    P_fa = P_fas(P_fa_idx);
+    if exist(filename, "file")
+        load( ...
+            filename, ...
+            "recovery_results", "sigma_w2", "thresholds",...
+            "C_n", "A", "Lambda", "Lambda_C", "x" ...
+        );
+    else
+        % Mento Carlo Recovery
+        [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, P_fa, noise_sigma / FAR_M);
+    
+        save( ...
+            filename, ...
+            "recovery_results", "sigma_w2", "thresholds",...
+            "C_n", "A", "Lambda", "Lambda_C", "x" ...
+        );
+    end
+    
+    % Compare sigmas (from paper) and vars (from simulation)
+    % [ss, vs, H0_mean, H0_max, H1_mean, H1_max] = get_sigma_var(recovery_results, sigma_w2, x);
+    
+    T_opt = zeros(2, 1);
+    
+    delta_r = c / (2 * B);
+    tmp_1 = repelem(f_n, FAR_M);
+    tmp_2 = repmat(1: FAR_M, 1, FAR_N)';
+    tmp_3 = exp(-1 * 1j * 2 * pi * 2 * delta_r .* tmp_2 .* tmp_1 / c);
+    
+    idx_1 = 1;
+    idx_2 = 1;
+    i1 = 1; i2 = 1;
+
+    Lambda_distributes = zeros(2, length(Lambda) * trail_times);
+    Lambda_C_distributes = zeros(2, length(Lambda_C) * trail_times);
+
+    for t = 1:trail_times
+        % 0 假设
+        x_0_hat = squeeze(recovery_results(1, t, :));
+        s = x_0_hat .* tmp_3;
+        for i = 1: FAR_N
+            L = (i - 1) * FAR_M + 1;
+            R = L - 1 + FAR_M;
+            T_opt(1, idx_1) = foo(s(L:R));
+            idx_1 = idx_1 + 1;
+        end
+
+        % 1 假设
+        x_1_hat = squeeze(recovery_results(2, t, :));
+        s = x_1_hat .* tmp_3;
+        T_opt(2, idx_2) = foo(s(Lambda));
+        idx_2 = idx_2 + 1;
+
+        for i = 1: length(Lambda) + length(Lambda_C)
+            if ismember(i, Lambda_C)
+                Lambda_C_distributes(1, i1) = x_0_hat(i);
+                Lambda_C_distributes(2, i1) = x_1_hat(i);
+                i1 = i1 + 1;
+            elseif ismember(i, Lambda)
+                Lambda_distributes(1, i2) = x_0_hat(i);
+                Lambda_distributes(2, i2) = x_1_hat(i);
+                i2 = i2 + 1;
+            end
+        end
+    end
+    
+    % figure;
+    % subplot(211); plot(real(sum(squeeze(recovery_results(1, 1:10, :)))));
+    % subplot(212); plot(real(sum(squeeze(recovery_results(2, 1:10, :)))));
+    
+    % figure;
+    % subplot(211); histfit(T_opt(1, 1:idx_1 - 1));
+    % subplot(212); histfit(T_opt(2, 1:idx_2 - 1));
+    
+    figure;
+    subplot(2, 2, 1); histfit(real(Lambda_C_distributes(1, :)));
+    subplot(2, 2, 2); histfit(real(Lambda_distributes(1, :)));
+    subplot(2, 2, 3); histfit(real(Lambda_C_distributes(2, :)));
+    subplot(2, 2, 4); histfit(real(Lambda_distributes(2, :)));
+
+    figure;
+    histfit(real([Lambda_distributes(2, :) Lambda_C_distributes(2, :)]));
+
+    means = [ ...
+                 mean(Lambda_C_distributes(1, :)), ...
+                 mean(Lambda_distributes(1, :)), ...
+                 mean(Lambda_C_distributes(2, :)), ...
+                 mean(Lambda_distributes(2, :)), ...
+                 mean([Lambda_distributes(1, :) Lambda_C_distributes(1, :)]), ...
+                 mean([Lambda_distributes(1, :) Lambda_C_distributes(1, :) Lambda_C_distributes(2, :)]), ...
+                 mean([Lambda_distributes(2, :)]) ...
+             ];
+
+    stds = [ ...
+                std(Lambda_C_distributes(1, :)), ...
+                std(Lambda_distributes(1, :)), ...
+                std(Lambda_C_distributes(2, :)), ...
+                std(Lambda_distributes(2, :)), ...
+                 std([Lambda_distributes(1, :) Lambda_C_distributes(1, :)]), ...
+                 std([Lambda_distributes(1, :) Lambda_C_distributes(1, :) Lambda_C_distributes(2, :)]), ...
+                 std([Lambda_distributes(2, :)]) ...
+            ];
+
+
+    fprintf("H_00: mu = %.15f, std = %.15f\n", means(1), stds(1));
+    fprintf("H_01: mu = %.15f, std = %.15f\n", means(2), stds(2));
+    fprintf("H_10: mu = %.15f, std = %.15f\n", means(3), stds(3));
+    fprintf("H_11: mu = %.15f, std = %.15f\n", means(4), stds(4));
+    fprintf("H_0: mu = %.15f, std = %.15f\n", means(5), stds(5));
+    fprintf("H_0 + Lambda^C: mu = %.15f, std = %.15f\n", means(6), stds(6));
+    fprintf("H_1 Lambda: mu = %.15f, std = %.15f\n\n", means(7), stds(7));
+    
+
+    H_0_mean = mean(T_opt(1, 1:idx_1-1));
+    H_0_std = std(T_opt(1, 1:idx_1-1));
+    H_1_mean = mean(T_opt(2, 1:idx_2-1));
+    H_1_std = std(T_opt(2, 1:idx_2-1));
+
+
+    threshold = norminv(1 - P_fa, H_0_mean, H_0_std);
+    P_d = 1 - normcdf(threshold, H_1_mean, H_1_std);
+    P_d_FARs(P_fa_idx) = P_d;
+end
+
+format long;
+fprintf("H_0 mean: %.15f\n", H_0_mean);
+fprintf("H_0 std: %.15f\n", H_0_std);
+fprintf("H_1 mean: %.15f\n", H_1_mean);
+fprintf("H_1 std: %.15f\n", H_1_std);
+
+%% PD
+[s_T_narrow, A] = narrow_signal_model_2(FAR_N);
+[Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, 1, betas_narrow, false);
+filename = ...
+    "data2/" + ...
+    "PD" + "_" + ...
+    method + "_" + ...
+    string(noise_sigma) + "_" + ...
+    string(length(Lambda)) + "_" + ...
+    trail_times + ...
+    ".mat";
+
+for P_fa_idx = 1: length(P_fas)
+    [s_T, A] = narrow_signal_model_2(FAR_N);
+    P_fa = P_fas(P_fa_idx);
+    if exist(filename, "file")
+        load( ...
+            filename, ...
+            "recovery_results", "sigma_w2", "thresholds",...
+            "f_c", "A", "Lambda", "Lambda_C", "x" ...
+        );
+    else
+        % Mento Carlo Recovery
+        [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, P_fa, noise_sigma);
+    
+        save( ...
+            filename, ...
+            "recovery_results", "sigma_w2", "thresholds",...
+            "f_c", "A", "Lambda", "Lambda_C", "x" ...
+        );
+    end
+
+    T_opt = zeros(2, trail_times);
+        
+    delta_r = T_p * c / 2;
+    tmp_1 = (1: FAR_N)';
+    tmp_3 = exp(-1 * 1j * 2 * pi * 2 * delta_r * f_c .* tmp_1 / c);
+
+    for t = 1:trail_times
+        % 0 假设
+        x_hat = recovery_results(1, t, :);
+        s = squeeze(x_hat) .* tmp_3;
+        for i = 1: FAR_N
+            T_opt(1, idx_1) = foo(s(i));
+            idx_1 = idx_1 + 1;
+        end
+
+        % 1 假设
+        x_hat = recovery_results(2, t, :);
+        s = squeeze(x_hat) .* tmp_3;
+        T_opt(2, idx_2) = foo(s(Lambda));
+        idx_2 = idx_2 + 1;
+    end
+
+    % figure;
+    % subplot(211); plot(real(sum(squeeze(recovery_results(1, 1:10, :)))));
+    % subplot(212); plot(real(sum(squeeze(recovery_results(2, 1:10, :)))));
+
+    H_0_mean = mean(T_opt(1, 1:idx_1-1));
+    H_0_std = std(T_opt(1, 1:idx_1-1));
+    H_1_mean = mean(T_opt(2, 1:idx_2-1));
+    H_1_std = std(T_opt(2, 1:idx_2-1));
+    
+    threshold = norminv(1 - P_fa, H_0_mean, H_0_std);
+    P_d = 1 - normcdf(threshold, H_1_mean, H_1_std);
+    P_d_narrows(P_fa_idx) = P_d;
+end
+
+format long;
+fprintf("H_0 mean: %.15f\n", H_0_mean);
+fprintf("H_0 std: %.15f\n", H_0_std);
+fprintf("H_1 mean: %.15f\n", H_1_mean);
+fprintf("H_1 std: %.15f\n", H_1_std);
+
+
+%% ROC
+
+figure;
+hold on;
+grid on;
+semilogx(P_fas, P_d_FARs);
+semilogx(P_fas, P_d_narrows);
+xlim([8e-8, 1]);
+% ylim([0, 1]);
+xlabel("P_{fa}");
+ylabel("P_{d}");
+title("ROC");

wide vs narrow/PT curve/draw_pt_curve.m

@@ -0,0 +1,11 @@
+clc; clear; close all;
+
+N = 256;
+M = 4;
+
+[K_range, N_b, N_s] = get_phase_transition_curve(N, M);
+
+figure;
+plot(K_range, N_b);
+hold on;
+plot(K_range, N_b);

wide vs narrow/PT curve/get_integral.m

@@ -0,0 +1,20 @@
+function I = get_integral(cache_filename, var_name, tau_range, M)
+    syms x f;
+
+    if exist(cache_filename, "file")
+        data = load(cache_filename, var_name);
+        I = data.(var_name);
+    else
+        I = zeros(length(tau_range), 0);
+
+        for tau_idx = 1:length(tau_range)
+            tau = tau_range(tau_idx);
+            f = (x - tau) ^ 2 * exp(-x ^ 2/2) * x ^ (M - 1) / (2 ^ (M / 2 - 1) * gamma(M / 2));
+            I(tau_idx) = double(int(f, [tau, +inf]));
+        end
+
+        S.(var_name) = I;
+        save(cache_filename, '-struct', 'S', var_name);
+    end
+
+end

wide vs narrow/PT curve/get_phase_transition_curve.m

@@ -0,0 +1,47 @@
+function [K_range, N_b, N_s] = get_phase_transition_curve(N, M, K_max)
+    % This function give the phase transition curve of FAR
+    % input: N, M: FAR parameters
+    % output: N_b, N_s: the phase transition curve of FAR in block case and scalar case
+
+    if nargin < 3
+        K_max = 25;
+    end
+
+    tau_min = 0;
+    tau_max = 15;
+    tau_interval = 0.05;
+    tau_range = tau_min:tau_interval:tau_max;
+
+    K_range = 1:K_max;
+
+    N_b = zeros(length(K_range), 1);
+    N_s = zeros(length(K_range), 1);
+
+    cache_filename_1 = ...
+        "../data/PT_curve_data/" + ...
+        "I_2.mat";
+
+    cache_filename_2 = ...
+        "../data/PT_curve_data/" + ...
+        "I_" + "" + string(2 * M) + ...
+        ".mat";
+
+    I_2 = get_integral(cache_filename_1, "I_2", tau_range, 2);
+    I_2M = get_integral(cache_filename_2, "I_2M", tau_range, 2 * M);
+
+    for K_idx = 1:length(K_range)
+        K = K_range(K_idx);
+        f_set_1 = zeros(length(tau_range), 1);
+        f_set_2 = zeros(length(tau_range), 1);
+
+        for tau_idx = 1:length(tau_range)
+            tau = tau_range(tau_idx);
+            f_set_1(tau_idx) = 1/2 * (K * (2 * M + tau ^ 2) + (N - K) * I_2M(tau_idx));
+            f_set_2(tau_idx) = M / 2 * (K * (2 + tau ^ 2) + (N - K) * I_2(tau_idx));
+        end
+
+        N_b(K_idx) = min(f_set_1);
+        N_s(K_idx) = min(f_set_2);
+    end
+
+end

wide vs narrow/single_point_test/All_Recovery.m

@@ -0,0 +1,104 @@
+function [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, P_fa, noise_sigma)
+
+    global trail_times method LASSO_lambda tau iter_max;
+    sz = size(A);
+    M = sz(1);
+    N = sz(2);
+
+    recovery_results = zeros(2, trail_times, N);
+    sigma_w2 = zeros(2, trail_times, 1);
+    thresholds = zeros(2, trail_times, 1);
+    figure;
+    
+    for hypo = 2: -1: 1
+        % hypo-假设
+        hypo = 3 - hypo;
+        h = waitbar(0, '正在仿真' + string(hypo-1) + '假设情况');
+        for T = 1:trail_times
+            waitbar(T / trail_times, h);
+    
+            noise = get_noise(noise_sigma, M, 1);
+
+            % y = Ax + n
+            if hypo == 1
+                y_noise = noise;
+            else
+                y_noise = A * x + noise;
+            end
+
+            if method == "debiased_LASSO"
+                [x_hat, sigma_w_2, threshold] = debiased_LASSO(A, y_noise, P_fa, noise_sigma^2, LASSO_lambda);
+                sigma_w2(hypo, T) = sigma_w_2;
+                thresholds(hypo, T) = threshold;
+            elseif method == "debiased_LASSO_FISTA"
+                [x_hat, sigma_w_2, threshold] = debiased_LASSO_FISTA(A, y_noise, P_fa, noise_sigma^2, LASSO_lambda);
+                sigma_w2(hypo, T) = sigma_w_2;
+                thresholds(hypo, T) = threshold;
+            elseif method == "cVAMPro"
+                [x_LASSO, x_hat_d] = cVAMPro(y_noise, A, LASSO_lambda, tau, 100);
+                % x_LASSO = FISTA(y_noise, A, LASSO_lambda, 1e-5);
+
+
+                sz = size(A);
+                n = sz(2);
+                gamma = sz(1) / sz(2);
+                lambda = LASSO_lambda;
+
+                % CROD求去偏
+                rho_active = sum(abs(x_LASSO) > 1e-3)/n;
+                Q_hat = (gamma - rho_active)/(1 - rho_active);
+                Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n;
+                diff = 1;
+                while(diff > 1e-4)
+                    Rho_pre = Rho;
+                    Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n;
+                    diff = abs(Rho - Rho_pre);
+                end
+                Q_hat = (gamma-Rho)/(1-Rho);
+                x_d_CROD = x_LASSO + A'*(y_noise - A*x_LASSO)/Q_hat; % x_d_CROD == x_hat_d
+                sigma_n = noise_sigma;
+            
+                % CROD求门限和检验统计量
+                RSS = sum(abs(y_noise - A * x_LASSO).^2)/length(y_noise);
+                chi = Rho*(1 - Rho)/(gamma - Rho);
+                if chi ~= 0
+                    chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi);
+                    z = -(1 - chi + chi_temp) / (2*chi);
+                    z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp);
+                    G_prime = (z + 1/chi);
+                    G_wprime = (z_prime + 1/chi/chi);
+                    chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
+                        + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
+                else
+                    G_prime = gamma;
+                    G_wprime = gamma*(1-gamma);
+                    chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
+                        + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
+                end
+                sigma_CROD = sqrt(2*chi_hat) / Q_hat;
+
+                sigma_w2(hypo, T) = sigma_CROD^2;
+                thresholds(hypo, T) = -sigma_CROD^2 * log(P_fa);
+                % figure; subplot(221); plot(real(x)); subplot(222); plot(real(x_LASSO)); subplot(223); plot(real(x_hat_d)); subplot(224); plot(real(x_d_CROD));
+                if hypo == 2
+                    x_hat = x_hat_d;
+                    % pp = real(x_hat - x);
+                    % [is_not_norm, tmp, tmp] = swtest(pp, 0.1);
+                    % 
+                    % if is_not_norm
+                    %     subplot(211); plot(pp); subplot(212); histfit(pp);
+                    %     is_not_norm
+                    % end
+                else
+                    x_hat = x_d_CROD;
+                end
+            else
+                x_hat = recovery(A, y_noise, method);
+            end
+
+    
+            recovery_results(hypo, T, :) = x_hat;
+        end
+        close(h);
+    end
+end

wide vs narrow/single_point_test/All_Recovery2.m

@@ -0,0 +1,95 @@
+function [recovery_results, sigma_w2, thresholds] = All_Recovery2(A, x, P_fa, noise_sigma, trail_times, LASSO_lambda, tau)
+
+    sz = size(A);
+    M = sz(1);
+    N = sz(2);
+
+    recovery_results = zeros(2, trail_times, N);
+    sigma_w2 = zeros(2, trail_times, 1);
+    thresholds = zeros(2, trail_times, 1);
+    
+    for hypo = 2: -1: 1
+        % hypo-假设
+        hypo = 3 - hypo;
+        h = waitbar(0, '正在仿真' + string(hypo-1) + '假设情况');
+        for T = 1:trail_times
+            waitbar(T / trail_times, h);
+    
+            noise = get_noise(noise_sigma, M, 1);
+
+            % y = Ax + n
+            if hypo == 1
+                y_noise = noise;
+            else
+                y_noise = A * x + noise;
+            end
+
+            [x_LASSO, x_hat_d] = cVAMPro(y_noise, A, LASSO_lambda, tau, 100);
+            % x_LASSO = FISTA(y_noise, A, LASSO_lambda, 1e-5);
+
+
+            sz = size(A);
+            n = sz(2);
+            gamma = sz(1) / sz(2);
+            lambda = LASSO_lambda;
+
+            % CROD求去偏
+            rho_active = sum(abs(x_LASSO) > 1e-3)/n;
+            Q_hat = (gamma - rho_active)/(1 - rho_active);
+            Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n;
+            diff = 1;
+            while(diff > 1e-4)
+                Rho_pre = Rho;
+                Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n;
+                diff = abs(Rho - Rho_pre);
+            end
+            Q_hat = (gamma-Rho)/(1-Rho);
+            x_d_CROD = x_LASSO + A'*(y_noise - A*x_LASSO)/Q_hat; % x_d_CROD == x_hat_d
+            sigma_n = noise_sigma;
+        
+            % CROD求门限和检验统计量
+            RSS = sum(abs(y_noise - A * x_LASSO).^2)/length(y_noise);
+            chi = Rho*(1 - Rho)/(gamma - Rho);
+            if chi ~= 0
+                chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi);
+                z = -(1 - chi + chi_temp) / (2*chi);
+                z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp);
+                G_prime = (z + 1/chi);
+                G_wprime = (z_prime + 1/chi/chi);
+                chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
+                    + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
+            else
+                G_prime = gamma;
+                G_wprime = gamma*(1-gamma);
+                chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)...
+                    + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime);
+            end
+            sigma_CROD = sqrt(2*chi_hat) / Q_hat;
+
+            if sigma_CROD > 1
+                sigma_CROD = sigma_CROD;
+            end
+            sigma_w2(hypo, T) = sigma_CROD^2;
+            thresholds(hypo, T) = -sigma_CROD^2 * log(P_fa);
+            % figure; subplot(221); plot(real(x)); subplot(222); plot(real(x_LASSO)); subplot(223); plot(real(x_hat_d)); subplot(224); plot(real(x_d_CROD));
+
+            % % data4 folder
+            % if hypo == 2
+            %     x_hat = x_d_CROD;
+            %     % x_hat = x_hat_d;
+            % else
+            %     x_hat = x_d_CROD;
+            % end
+
+            % data5 folder
+            if hypo == 2
+                x_hat = x_hat_d;
+            else
+                x_hat = x_d_CROD;
+            end
+
+            recovery_results(hypo, T, :) = x_hat;
+        end
+        close(h);
+    end
+end

wide vs narrow/single_point_test/get_lambda.m

@@ -0,0 +1,19 @@
+function [LASSO_lambda, recovery_results] = get_lambda(A, x, FAR_N, FAR_M, sigma, lambdas)
+
+    noise = get_noise(sigma, FAR_N, 1);
+    y = A * x;
+    y_noise = y + noise;
+
+    tau = 1e-6;
+    iter_max = 1000;
+
+    recovery_results = zeros(length(lambdas), 1);
+    for lambda_idx = 1: length(lambdas)
+        LASSO_lambda = lambdas(lambda_idx);
+        [x_LASSO, x_hat_d] = cVAMPro(y_noise, A, LASSO_lambda, tau, iter_max);
+        recovery_results(lambda_idx) = sum(real(x_LASSO(FAR_M + 1: 2 * FAR_M)));
+        % if abs(recovery_results(lambda_idx) - FAR_M) / FAR_M < 0.01
+        %     return
+        % end
+    end
+end

wide vs narrow/single_point_test/get_sigma_var.m

@@ -0,0 +1,71 @@
+function [sigmas, vars, H0_mean, H0_max, H1_mean, H1_max] = get_sigma_var(recovery_results, sigma_w2, x)
+
+    sz = size(recovery_results);
+    trail_times = sz(2);
+
+    sigmas = zeros(2, trail_times);
+    vars = zeros(2, trail_times);
+    id = 1;
+
+    for T = 1:trail_times
+
+        x_0_hat = squeeze(recovery_results(1, T, :));
+        x_1_hat = squeeze(recovery_results(2, T, :)) - x;
+
+        if anynan(x_1_hat)
+            continue;
+        end
+
+        % if sigma_w2(2, T) > 0.1
+        %     continue;
+        % end
+
+        sigmas(1, id) = sigma_w2(1, T);
+        sigmas(2, id) = sigma_w2(2, T);
+
+        vars(1, id) = var(x_0_hat);
+        vars(2, id) = var(x_1_hat);
+
+        id = id + 1;
+    end
+    sigmas(:, end - trail_times + id: end) = [];
+    vars(:, end - trail_times + id: end) = [];
+
+    REE = abs(vars - sigmas) ./ vars;
+    H0_mean = mean(REE(1, :));
+    H0_max = max(REE(1, :));
+    H1_mean = mean(REE(2, :));
+    H1_max = max(REE(2, :));
+    
+    if 1 == 9
+        figure;
+        subplot(211);
+        plot(sigmas(1, :));
+        hold on;
+        plot(vars(1, :));
+        xlabel("Mento Carlo Times");
+        ylabel("\sigma^2")
+        legend("\sigma^2 (CROD)", "\sigma_w^2 (Expr)");
+        subplot(212);
+        plot(sigmas(2, :));
+        hold on;
+        plot(vars(2, :));
+        xlabel("Mento Carlo Times");
+        ylabel("\sigma^2");
+        legend("\sigma^2 (CROD)", "\sigma_w^2 (Expr)");
+    
+        figure;
+        subplot(211);
+        plot(REE(1, :));
+        yline(H0_mean);
+        xlabel("Mento Carlo Times");
+        ylabel("REE");
+        legend("mean = " + string(H0_mean), "max = " + string(H0_max));
+        subplot(212);
+        plot(REE(2, :));
+        yline(H1_mean);
+        xlabel("Mento Carlo Times");
+        ylabel("REE");
+        legend("mean = " + string(H1_mean), "max = " + string(H1_max));
+    end
+end

wide vs narrow/single_point_test/get_sparse_vector.m

@@ -0,0 +1,28 @@
+function [Lambda, Lambda_C, x] = get_sparse_vector(N, M, Amp, extend_target)
+    x = zeros(N * M, 1);
+
+    if nargin < 4
+        extend_target = true;
+    end
+
+    % Support Set
+    if length(x) > 2 * M && extend_target
+        Lambda = M + 1:2 * M;
+    elseif length(x) < 2 * M && extend_target
+        Lambda = 1:M;
+    else
+        Lambda = 2;
+    end
+
+    if length(Amp) == 1
+        x(Lambda) = Amp;
+    elseif length(Amp) <= max(Lambda)
+        x(Lambda) = Amp(Lambda);
+    else
+        % Warning
+        x(Lambda) = Amp(end);
+    end
+
+    Lambda_C = setdiff(1:N*M, Lambda);
+
+end

wide vs narrow/single_point_test/parameters.m

@@ -0,0 +1,88 @@
+global ...
+    B ...
+    T_p ...
+    f_s ...
+    T_s ...
+    f_c ...
+    PRF ...
+    T_r ...
+    t_single_pulse ...
+    c ...
+    N ...
+    M ...
+    R_0...
+    delta_R_wide ...
+    delta_R_narrow ...
+    N_wide ...
+    N_narrow ...
+    ranges_wide...
+    ranges_narrow ...
+    betas_wide ...
+    betas_narrow ...
+    trail_times...
+    method ...
+    N_high ...
+    LASSO_lambda ...
+    tau ...
+    iter_max ...
+    FAR_N ...
+    FAR_M ...
+    Ns Ms lambdas sigmas
+;
+
+B = 10e5;
+T_p = 1e-5;
+f_s = B;
+T_s = 1 / f_s;
+f_c = 1e9;
+duty_cycle_inv = 10;
+T_r = T_p * duty_cycle_inv;
+PRF = 1 / T_r;
+
+t_single_pulse = 0: 1 / f_s: T_r - 1 / f_s;
+
+c = 3e8;
+
+N = length(t_single_pulse);
+N_high = T_p * f_s;
+M = round(T_p * B);
+R_0 = 0;
+
+delta_R_wide = c ./ 2 ./ B; % 宽带情况下的距离分辨力
+delta_R_narrow = T_p .* c ./ 2; % 窄带情况下的距离分辨力
+
+N_wide = FAR_N; % 宽带情况下，发射 480 个脉冲
+
+N_narrow = round(N_wide / M);
+
+ranges_wide = R_0 + (0:N_wide) * delta_R_wide; % [1000, 4000]
+ranges_narrow = R_0 + (0:N_narrow) * delta_R_narrow; % [1000, 4000]
+
+% betas_wide = (linspace(1, 0.1, FAR_M * FAR_N) + 1j * linspace(0.1, 1, FAR_M * FAR_N))';
+betas_wide = zeros(FAR_M * FAR_N, 1) + 1;
+betas_narrow = zeros(FAR_N, 1);
+for i = 1: FAR_N
+    for j = 1: FAR_M
+        idx = (i-1) * FAR_M + j;
+        betas_narrow(i) = betas_narrow(i) + betas_wide(idx) * exp(1j * 2 * pi * 2 * delta_R_wide * j / c);
+    end
+end
+
+trail_times = 250;
+% method = "debiased_LASSO";
+% method = "LASSO";
+% method = "BP";
+method = "cVAMPro";
+% method = "debiased_LASSO_FISTA";
+
+LASSO_lambda = 0.1;
+tau = 1e-6;
+iter_max = 1000;
+
+Ns = [256];
+Ms = [16];
+sigmas = [0.01];
+lambdas = zeros(length(Ns), 1);
+for i = 1: length(Ns)
+    lambdas(i) = Map(Ns(i) + "_" + Ms(i));
+end

wide vs narrow/single_point_test/query2.m

@@ -0,0 +1,39 @@
+function [recovery_results, sigma_w2, thresholds, C_n, A, Lambda, Lambda_C, x] = query2(FAR_N, FAR_M, sigma, LASSO_lambda)
+    global trail_times method;
+
+    if FAR_N == 2
+        lambda_length = 1;
+    else
+        lambda_length = FAR_M;
+    end
+    
+    filename = ...
+        "./data3/" + ...
+        "FAR" + "_" + ...
+        string(FAR_N) + "_" + ...
+        string(FAR_M)+ "_" + ...
+        method + "_" + ...
+        string(sigma) + "_" + ...
+        string(lambda_length) + "_" + ...
+        string(LASSO_lambda) + "_" + ...
+        trail_times + ...
+        ".mat";
+
+    filename = ...
+        "./data3/" + ...
+        "FAR" + "_" + ...
+        string(FAR_N) + "_" + ...
+        string(FAR_M)+ "_" + ...
+        method + "_" + ...
+        string(sigma) + "_" + ...
+        trail_times + ...
+        ".mat";
+
+    load( ...
+        filename, ...
+        "recovery_results", "sigma_w2", "thresholds",...
+        "C_n", "A", "Lambda", "Lambda_C", "x" ...
+    );
+
+    return; 
+end

wide vs narrow/single_point_test/recovery.m

@@ -0,0 +1,52 @@
+function x_hat = recovery(A, y_noise, method)
+    if method == "debiased_LASSO"
+        sz = size(A);
+        N = sz(2);
+        LASSO_lambda = 0.1;
+
+        gamma = sz(1) / sz(2);
+
+        cvx_begin quiet
+            variable x_LASSO(N) complex
+            minimize(LASSO_lambda * norm(x_LASSO, 1) + norm(y_noise - A * x_LASSO, 2))
+        cvx_end
+
+        rho_active = sum(abs(x_LASSO) > 1e-3)/N;
+        Q_hat = (gamma - rho_active)/(1 - rho_active);
+        Rho = sum((abs(x_LASSO) > 1e-3).* (2 - LASSO_lambda./(Q_hat*abs(x_LASSO) + LASSO_lambda))) / 2 / N;
+        diff = 1;
+        while(diff > 1e-4)
+            Rho_pre = Rho;
+            Rho = sum((abs(x_LASSO) > 1e-3).* (2 - LASSO_lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + LASSO_lambda))) / 2 / N;
+            diff = abs(Rho - Rho_pre);
+        end
+        Q_hat = (gamma-Rho)/(1-Rho);
+        x_d_CROD = x_LASSO + A'*(y_noise - A*x_LASSO)/Q_hat;
+        x_hat = x_d_CROD;
+
+    elseif method == "LASSO"
+        sz = size(A);
+        N = sz(2);
+        LASSO_lambda = 0.1;
+        cvx_begin quiet
+            variable x_LASSO(N) complex
+            minimize(LASSO_lambda * norm(x_LASSO, 1) + norm(y_noise - A * x_LASSO, 2))
+        cvx_end
+
+        x_hat = x_LASSO;
+
+    elseif method == "BP"
+        sz = size(A);
+        N = sz(2);
+        cvx_begin quiet
+            variable x_hat(N) complex
+            minimize(norm(x_hat, 1))
+        subject to
+            A * x_hat == y_noise
+        cvx_end
+    else
+        global lambda tau iter_max;
+        [x_hat_wl, x_hat] = cVAMPro(y_noise, A, lambda, tau, iter_max);
+    end
+end
+

wide vs narrow/single_point_test/signle_point_mento_carlo_test.m

@@ -0,0 +1,269 @@
+% 单目标检测点测试
+
+clc; clear; close all;
+
+parameters;
+
+% Ns = [4, 8, 16, 32, 64, 128];
+
+TEST_1 = true;
+TEST_2 = false;
+TEST_3 = false;
+TEST_4 = false;
+
+if TEST_1
+    for N_idx = 1: length(Ns)
+        for M_idx = 1: length(Ms)
+            for sigma_idx = 1: length(sigmas)
+                FAR_N = Ns(N_idx);
+                FAR_M = Ms(M_idx);
+                sigma = sigmas(sigma_idx);
+                for lambda_idx = 1: length(lambdas)
+                    LASSO_lambda = lambdas(lambda_idx);
+    
+                    % epi = 0;
+                    % [C_n, A] = get_Psi(FAR_N, FAR_M, epi);
+                    % fac = A * A';
+                    % A = A ./ sqrt(abs(fac(1, 1)));
+                    % f_n = f_c + C_n * B / M;
+                    % 
+                    % [Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true);
+                    % [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, LASSO_lambda, sigma);
+            
+                    [recovery_results, sigma_w2, thresholds, C_n, A, Lambda, Lambda_C, x] = query2(FAR_N, FAR_M, sigma, LASSO_lambda);
+            
+                    i1 = 1; i2 = 1;
+                    Lambda_distributes = zeros(2, length(Lambda) * trail_times);
+                    Lambda_C_distributes = zeros(2, length(Lambda_C) * trail_times);
+                    x_1_add = length(length(Lambda) + length(Lambda_C));
+            
+                    for t = 1:trail_times
+                        x_0_hat = squeeze(recovery_results(1, t, :));
+                        x_1_hat = squeeze(recovery_results(2, t, :)) - x;
+                        if anynan(x_1_hat)
+                            continue;
+                        end
+            
+                        for i = 1: length(Lambda) + length(Lambda_C)
+                            if ismember(i, Lambda_C)
+                                Lambda_C_distributes(1, i1) = x_0_hat(i);
+                                Lambda_C_distributes(2, i1) = x_1_hat(i);
+                                i1 = i1 + 1;
+                            elseif ismember(i, Lambda)
+                                Lambda_distributes(1, i2) = x_0_hat(i);
+                                Lambda_distributes(2, i2) = x_1_hat(i);
+                                i2 = i2 + 1;
+                            end
+                        end
+                    end
+                    real_H00 = real(Lambda_C_distributes(1, 1:i1-1));
+                    real_H01 = real(Lambda_distributes(1, 1:i2-1));
+                    real_H10 = real(Lambda_C_distributes(2, 1:i1-1));
+                    real_H11 = real(Lambda_distributes(2, 1:i2-1));
+
+                    figure;
+                    subplot(2, 2, 1); histfit(real_H00); xlabel("x\_hat"); ylabel("times");
+                    subplot(2, 2, 2); histfit(real_H01); xlabel("x\_hat"); ylabel("times");
+                    subplot(2, 2, 3); histfit(real_H10); xlabel("x\_hat"); ylabel("times");
+                    subplot(2, 2, 4); histfit(real_H11); xlabel("x\_hat"); ylabel("times");
+                   
+                    % figure;
+                    % data = real(Lambda_C_distributes(2, :));
+                    % xx = linspace(min(data), max(data), 1e3);
+                    % yy = normpdf(x, mean(data), std(data));
+                    % hold on;
+                    % histfit(data);
+            
+                    xlabel("x\_hat");
+                    ylabel("times");
+                    legend("lambda = " + string(LASSO_lambda) + ", sigma = " + string(sigma));
+    
+                    means = [ ...
+                                 mean(Lambda_C_distributes(1, 1:i1-1)), ...
+                                 mean(Lambda_distributes(1, i2-1)), ...
+                                 mean(Lambda_C_distributes(2, 1:i1-1)), ...
+                                 mean(Lambda_distributes(2, i2-1)), ...
+                                 mean([Lambda_distributes(1, i2-1) Lambda_C_distributes(1, 1:i1-1)]) ...
+                             ];
+                
+                    stds = [ ...
+                                std(Lambda_C_distributes(1, 1:i1-1)), ...
+                                std(Lambda_distributes(1, i2-1)), ...
+                                std(Lambda_C_distributes(2, 1:i1-1)), ...
+                                std(Lambda_distributes(2, i2-1)), ...
+                                std([Lambda_distributes(1, i2-1) Lambda_C_distributes(1, 1:i1-1)]) ...
+                            ];
+                
+                
+                    fprintf("H_00: mu = %.15f, std = %.15f\n", means(1), stds(1));
+                    fprintf("H_01: mu = %.15f, std = %.15f\n", means(2), stds(2));
+                    fprintf("H_10: mu = %.15f, std = %.15f\n", means(3), stds(3));
+                    fprintf("H_11: mu = %.15f, std = %.15f\n", means(4), stds(4));
+                    fprintf("H_0: mu = %.15f, std = %.15f\n", means(5), stds(5));
+                    
+                    [ss, vs, H0_mean, H0_max, H1_mean, H1_max] = get_sigma_var(recovery_results, sigma_w2, x);
+                end
+            end
+        end
+    end
+end
+
+if TEST_2
+    lg_lambdas = -5: 0.01: -0.75;
+    % lg_lambdas = -3: 0.01: -2.5;
+    lambdas = 10 .^ lg_lambdas;
+    LASSO_x_hat_maps = zeros(length(Ns), length(Ms), length(sigmas), length(lambdas));
+    LASSO_lambdas = zeros(length(Ns), length(Ms), length(sigmas));
+    for sigma_idx = 1: length(c)
+        sigma = sigmas(sigma_idx);
+        figure;
+        for N_idx = 1: length(Ns)
+            FAR_N = Ns(N_idx);
+            for M_idx = 1: length(Ms)
+                FAR_M = Ms(M_idx);
+                subplot(length(Ns), length(Ms), (N_idx-1)*length(Ms)+M_idx);
+                epi = 0;
+                [C_n, A] = get_Psi(FAR_N, FAR_M, epi);
+                fac = A * A';
+                A = A ./ sqrt(abs(fac(1, 1)));
+                f_n = f_c + C_n * B / M;
+                
+                [Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true);
+                [LASSO_lambda, LASSO_x_hat_map] = get_lambda(A, x, FAR_N, FAR_M, sigma, lambdas);
+                LASSO_x_hat_maps(N_idx, M_idx, sigma_idx, :) = LASSO_x_hat_map;
+                LASSO_lambdas(N_idx, M_idx, sigma_idx) = LASSO_lambda;
+
+                semilogx(lambdas, LASSO_x_hat_map / FAR_M);
+                xlabel("\lambda in LASSO");
+                ylabel("mean(x\_hat(\Lambda))");
+                yline(1);
+                ylim([-0.2, 1.2]);
+                legend("N = " + string(FAR_N), "M = " + string(FAR_M));
+            end
+        end
+    end
+    save("tmp.mat", lambdas, LASSO_x_hat_maps, LASSO_lambdas);
+end
+
+if TEST_3
+    for N_idx = 1: length(Ns)
+        for M_idx = 1: length(Ms)
+            for sigma_idx = 1: length(sigmas)
+                FAR_N = Ns(N_idx);
+                FAR_M = Ms(M_idx);
+                sigma = sigmas(sigma_idx);
+                for lambda_idx = 1: length(lambdas)
+                    LASSO_lambda = lambdas(lambda_idx);
+    
+                    % epi = 0;
+                    % [C_n, A] = get_Psi(FAR_N, FAR_M, epi);
+                    % fac = A * A';
+                    % A = A ./ sqrt(abs(fac(1, 1)));
+                    % f_n = f_c + C_n * B / M;
+                    % 
+                    % [Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true);
+                    % [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, LASSO_lambda, sigma);
+            
+                    [recovery_results, sigma_w2, thresholds, C_n, A, Lambda, Lambda_C, x] = query2(FAR_N, FAR_M, sigma, LASSO_lambda);
+            
+                    i1 = 1; i2 = 1;
+                    Lambda_distributes = zeros(2, length(Lambda) * trail_times);
+                    Lambda_C_distributes = zeros(2, length(Lambda_C) * trail_times);
+                    x_1_add = length(length(Lambda) + length(Lambda_C));
+            
+                    for t = 1:trail_times
+                        x_0_hat = squeeze(recovery_results(1, t, :));
+                        x_1_hat = squeeze(recovery_results(2, t, :)) - x;
+            
+                        for i = 1: length(Lambda) + length(Lambda_C)
+                            if ismember(i, Lambda_C)
+                                Lambda_C_distributes(1, i1) = x_0_hat(i);
+                                Lambda_C_distributes(2, i1) = x_1_hat(i);
+                                i1 = i1 + 1;
+                            elseif ismember(i, Lambda)
+                                Lambda_distributes(1, i2) = x_0_hat(i);
+                                Lambda_distributes(2, i2) = x_1_hat(i);
+                                i2 = i2 + 1;
+                            end
+                        end
+                    end
+                    real_H00 = real(Lambda_C_distributes(1, 1:i1-1));
+                    real_H01 = real(Lambda_distributes(1, 1:i2-1));
+                    real_H10 = real(Lambda_C_distributes(2, 1:i1-1));
+                    real_H11 = real(Lambda_distributes(2, 1:i2-1));
+    
+                    
+                    figure;
+                    subplot(2, 2, 1); histfit(real_H00); xlabel("x\_hat"); ylabel("times");
+                    subplot(2, 2, 2); histfit(real_H01); xlabel("x\_hat"); ylabel("times");
+                    subplot(2, 2, 3); histfit(real_H10); xlabel("x\_hat"); ylabel("times");
+                    subplot(2, 2, 4); histfit(real_H11); xlabel("x\_hat"); ylabel("times");
+                   
+                    figure;
+                    data = real(Lambda_C_distributes(2, :));
+                    xx = linspace(min(data), max(data), 1e3);
+                    yy = normpdf(x, mean(data), std(data));
+                    hold on;
+                    histfit(data);
+            
+                    xlabel("x\_hat");
+                    ylabel("times");
+                    legend("lambda = " + string(LASSO_lambda) + ", sigma = " + string(sigma));
+    
+                     means = [ ...
+                                 mean(Lambda_C_distributes(1, :)), ...
+                                 mean(Lambda_distributes(1, :)), ...
+                                 mean(Lambda_C_distributes(2, :)), ...
+                                 mean(Lambda_distributes(2, :)), ...
+                                 mean([Lambda_distributes(1, :) Lambda_C_distributes(1, :)]), ...
+                                 mean([Lambda_distributes(1, :) Lambda_C_distributes(1, :) Lambda_C_distributes(2, :)]), ...
+                                 mean([Lambda_distributes(2, :)]) ...
+                             ];
+                
+                    stds = [ ...
+                                std(Lambda_C_distributes(1, :)), ...
+                                std(Lambda_distributes(1, :)), ...
+                                std(Lambda_C_distributes(2, :)), ...
+                                std(Lambda_distributes(2, :)), ...
+                                 std([Lambda_distributes(1, :) Lambda_C_distributes(1, :)]), ...
+                                 std([Lambda_distributes(1, :) Lambda_C_distributes(1, :) Lambda_C_distributes(2, :)]), ...
+                                 std([Lambda_distributes(2, :)]) ...
+                            ];
+                
+                
+                    fprintf("H_00: mu = %.15f, std = %.15f\n", means(1), stds(1));
+                    fprintf("H_01: mu = %.15f, std = %.15f\n", means(2), stds(2));
+                    fprintf("H_10: mu = %.15f, std = %.15f\n", means(3), stds(3));
+                    fprintf("H_11: mu = %.15f, std = %.15f\n", means(4), stds(4));
+                    fprintf("H_0: mu = %.15f, std = %.15f\n", means(5), stds(5));
+                    fprintf("H_0 + Lambda^C: mu = %.15f, std = %.15f\n", means(6), stds(6));
+                    fprintf("H_1 Lambda: mu = %.15f, std = %.15f\n\n", means(7), stds(7));
+                    
+                    [ss, vs, H0_mean, H0_max, H1_mean, H1_max] = get_sigma_var(recovery_results, sigma_w2, x);
+                end
+            end
+        end
+    end
+end
+
+
+% for N_idx = 1: length(Ns)
+%     for M_idx = 1: length(Ms)
+%         for lambda_idx = 1: length(lambdas)
+%             for sigma_idx = 1: length(sigmas)
+%                 FAR_N = Ns(N_idx);
+%                 FAR_M = Ms(M_idx);
+%                 LASSO_lambda = lambdas(lambda_idx);
+%                 sigma = sigmas(sigma_idx);
+%                 [recovery_results, sigma_w2, thresholds, C_n, A, Lambda, Lambda_C, x] = query(FAR_N, FAR_M, sigma, LASSO_lambda);
+%             end
+%         end
+%     end
+% end
+
+
+% figure;
+% subplot(2, 2, 1); histfit(real(Lambda_C_distributes(1, :))); xlabel("x\_hat"); ylabel("times");
+% subplot(2, 2, 2); histfit(real(Lambda_distributes(1, :))); xlabel("x\_hat"); ylabel("times");
+% subplot(2, 2, 3); histfit(real(Lambda_C_distributes(2, :))); xlabel("x\_hat"); ylabel("times");
+% subplot(2, 2, 4); histfit(real(Lambda_distributes(2, :))); xlabel("x\_hat"); ylabel("times");

wide vs narrow/single_point_test/signle_point_mento_carlo_train.m

@@ -0,0 +1,53 @@
+% 单目标检测点测试
+
+clc; clear; close all;
+
+parameters;
+
+for N_idx = 1: length(Ns)
+    for M_idx = 1: length(Ms)
+        for sigma_idx = 1: length(sigmas)
+            FAR_N = Ns(N_idx);
+            FAR_M = Ms(M_idx);
+            % LASSO_lambda = lambdas(lambda_idx);
+            sigma = sigmas(sigma_idx);
+            for lambda_idx = 1: length(lambdas)
+                LASSO_lambda = lambdas(lambda_idx);
+            
+                epi = 0;
+                Psi_filename = "./results/Signal_Model_" + string(FAR_N) + "_" + string(FAR_M) + ".mat";
+                load(Psi_filename, "C_n", "A", "ref_lambda");
+                LASSO_lambda = ref_lambda;
+                % [C_n, A] = get_Psi(FAR_N, FAR_M, epi);
+                f_n = f_c + C_n * B / M;
+                
+                [Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true);
+                filename = ...
+                    "./data3/" + ...
+                    "FAR" + "_" + ...
+                    string(FAR_N) + "_" + ...
+                    string(FAR_M)+ "_" + ...
+                    method + "_" + ...
+                    string(sigma) + "_" + ...
+                    trail_times + ...
+                    ".mat";
+                if exist(filename, "file")
+                    load( ...
+                        filename, ...
+                        "recovery_results", "sigma_w2", "thresholds",...
+                        "C_n", "A", "Lambda", "Lambda_C", "x" ...
+                    );
+                else
+                    % Mento Carlo Recovery
+                    [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, LASSO_lambda, sigma);
+                
+                    save( ...
+                        filename, ...
+                        "recovery_results", "sigma_w2", "thresholds",...
+                        "C_n", "A", "Lambda", "Lambda_C", "x" ...
+                    );
+                end
+            end
+        end
+    end
+end

wide vs narrow/single_point_test/simulator.m

@@ -0,0 +1,196 @@
+%% Initial
+clc; clear; close all;
+
+trail_times = 250;
+method = "cVAMPro";
+Ns = [64, 128, 256, 512, 1024];
+Ms = [4, 8, 16, 32];
+sigma = 0.01;
+tau = 1e-6;
+iter_max = 100;
+epi = 0;
+
+lg_lambdas = -5: 0.05: -1;
+lambdas = 10 .^ lg_lambdas;
+
+H0_REE_means = zeros(length(Ns), length(Ms));
+H0_REE_maxs = zeros(length(Ns), length(Ms));
+H1_REE_means = zeros(length(Ns), length(Ms));
+H1_REE_maxs = zeros(length(Ns), length(Ms));
+
+figure;
+
+for N_idx = 1: length(Ns)
+    for M_idx = 1: length(Ms)
+        subplot(length(Ns), length(Ms), (N_idx - 1) * length(Ms) + M_idx);
+        FAR_N = Ns(N_idx);
+        FAR_M = Ms(M_idx);
+
+        fprintf("Simluating: N = %d, M = %d\n\n", FAR_N, FAR_M);
+
+        MN = FAR_N * FAR_M;
+        betas_wide = zeros(FAR_M * FAR_N, 1) + 1;
+        [Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true);
+
+        % get lambda
+        signal_model_filename = "./Signal_Model/Signal_Model_" + string(FAR_N) + "_" + string(FAR_M) + ".mat";
+        if exist(signal_model_filename, "file")
+            load(signal_model_filename, "A", "C_n", "ref_lambda");
+        else
+            curr_metric = FAR_N;
+            for TT = 1: 20
+                [C_n, A] = get_Psi(FAR_N, FAR_M, epi);
+                
+                noise = get_noise(sigma, FAR_N, 1);
+                y_noise = A * x + noise;
+                
+                MSEs = zeros(length(lambdas), 1);
+                all_x_hat = zeros(length(x), length(lambdas));
+                
+                for lambda_idx = 1: length(lambdas)
+                    lambda = lambdas(lambda_idx);
+                     [x_LASSO, x_hat_d] = cVAMPro(y_noise, A, lambda, tau, iter_max);
+                    MSEs(lambda_idx) = sum(real(x_LASSO - x));
+                    all_x_hat(:, lambda_idx) = x_LASSO;
+                end
+    
+                [metric, idx] = min(abs(MSEs));
+                if curr_metric > metric
+                    ref_lambda = lambdas(idx);
+                    save(signal_model_filename, "A", "C_n", "ref_lambda");
+                    curr_metric = metric;
+                end
+            end
+        end
+        fprintf("get lambda: lambda = %.15f. \n", ref_lambda);
+        fprintf("Signal Model filename: " + signal_model_filename + "\n\n");
+
+        % train
+        LASSO_lambda = ref_lambda;
+        simulate_results_filename = ...
+            "./data5/" + ...
+            "FAR" + "_" + ...
+            string(FAR_N) + "_" + ...
+            string(FAR_M)+ "_" + ...
+            method + "_" + ...
+            string(sigma) + "_" + ...
+            trail_times + ...
+            ".mat";
+        if exist(simulate_results_filename, "file")
+            load( ...
+                simulate_results_filename, ...
+                "recovery_results", "sigma_w2", "thresholds",...
+                "C_n", "A", "Lambda", "Lambda_C", "x" ...
+            );
+        else
+            % Mento Carlo Recovery
+            [recovery_results, sigma_w2, thresholds] = All_Recovery2(A, x, LASSO_lambda, sigma, trail_times, LASSO_lambda, tau);
+        
+            save( ...
+                simulate_results_filename, ...
+                "recovery_results", "sigma_w2", "thresholds",...
+                "C_n", "A", "Lambda", "Lambda_C", "x" ...
+            );
+        end
+        fprintf("Simulate results filename: " + simulate_results_filename + "\n\n");
+
+        % test
+        i1 = 1; i2 = 1;
+        Lambda_distributes = zeros(2, length(Lambda) * trail_times);
+        Lambda_C_distributes = zeros(2, length(Lambda_C) * trail_times);
+        
+        for t = 1:trail_times
+            x_0_hat = squeeze(recovery_results(1, t, :));
+            x_1_hat = squeeze(recovery_results(2, t, :)) - x;
+            if anynan(x_1_hat)
+                continue;
+            end
+
+            for i = 1: length(Lambda) + length(Lambda_C)
+                if ismember(i, Lambda_C)
+                    Lambda_C_distributes(1, i1) = x_0_hat(i);
+                    Lambda_C_distributes(2, i1) = x_1_hat(i);
+                    i1 = i1 + 1;
+                elseif ismember(i, Lambda)
+                    Lambda_distributes(1, i2) = x_0_hat(i);
+                    Lambda_distributes(2, i2) = x_1_hat(i);
+                    i2 = i2 + 1;
+                end
+            end
+        end
+        real_H00 = real(Lambda_C_distributes(1, 1:i1-1));
+        real_H01 = real(Lambda_distributes(1, 1:i2-1));
+        real_H10 = real(Lambda_C_distributes(2, 1:i1-1));
+        real_H11 = real(Lambda_distributes(2, 1:i2-1));
+
+        % figure;
+        % subplot(2, 2, 1); histfit(real_H00); xlabel("x\_hat"); ylabel("times"); title("H_0 (not in support set)")
+        % subplot(2, 2, 2); histfit(real_H01); xlabel("x\_hat"); ylabel("times"); title("H_0 (in support set)")
+        % subplot(2, 2, 3); histfit(real_H10); xlabel("x\_hat"); ylabel("times"); title("H_1 (not in support set)")
+        % subplot(2, 2, 4); histfit(real_H11); xlabel("x\_hat"); ylabel("times"); title("H_1 (in support set)")
+        % sgtitle("N = " + string(FAR_N) + ", M = " + string(FAR_M));
+
+        histfit(real_H11); xlabel("x\_hat"); ylabel("times"); title("N = " + string(FAR_N) + ", M = " + string(FAR_M));
+
+        means = [ ...
+                     mean(Lambda_C_distributes(1, 1:i1-1)), ...
+                     mean(Lambda_distributes(1, i2-1)), ...
+                     mean(Lambda_C_distributes(2, 1:i1-1)), ...
+                     mean(Lambda_distributes(2, i2-1)), ...
+                     mean([Lambda_distributes(1, i2-1) Lambda_C_distributes(1, 1:i1-1)]) ...
+                 ];
+    
+        stds = [ ...
+                    std(Lambda_C_distributes(1, 1:i1-1)), ...
+                    std(Lambda_distributes(1, i2-1)), ...
+                    std(Lambda_C_distributes(2, 1:i1-1)), ...
+                    std(Lambda_distributes(2, i2-1)), ...
+                    std([Lambda_distributes(1, i2-1) Lambda_C_distributes(1, 1:i1-1)]) ...
+                ];
+    
+    
+        fprintf("H_00: mu = %.15f, std = %.15f\n", means(1), stds(1));
+        fprintf("H_01: mu = %.15f, std = %.15f\n", means(2), stds(2));
+        fprintf("H_10: mu = %.15f, std = %.15f\n", means(3), stds(3));
+        fprintf("H_11: mu = %.15f, std = %.15f\n", means(4), stds(4));
+        fprintf("H_0: mu = %.15f, std = %.15f\n", means(5), stds(5));
+        
+        [ss, vs, H0_mean, H0_max, H1_mean, H1_max] = get_sigma_var(recovery_results, sigma_w2, x);
+        H0_REE_means(N_idx, M_idx) = H0_mean;
+        H0_REE_maxs(N_idx, M_idx) = H0_max;
+        H1_REE_means(N_idx, M_idx) = H1_mean;
+        H1_REE_maxs(N_idx, M_idx) = H1_max;
+        fprintf("Test Complete \n\n");
+    end
+end
+sgtitle("H_1 (in support set)");
+
+figure;
+
+subplot(221);
+h = heatmap(Ns, Ms, H0_REE_means');
+h.XLabel = "N";
+h.YLabel = "M";
+h.Title = "H_0, mean(REE)";
+
+subplot(222);
+h = heatmap(Ns, Ms, H0_REE_maxs');
+h.XLabel = "N";
+h.YLabel = "M";
+h.Title = "H_0, max(REE)";
+
+subplot(223);
+h = heatmap(Ns, Ms, H1_REE_means');
+h.XLabel = "N";
+h.YLabel = "M";
+h.Title = "H_1, mean(REE)";
+
+subplot(224);
+h = heatmap(Ns, Ms, H1_REE_maxs');
+h.XLabel = "N";
+h.YLabel = "M";
+h.Title = "H_1, max(REE)";
+
+
+
+

wide vs narrow/single_point_test/swtest.m

@@ -0,0 +1,273 @@
+function [H, pValue, W] = swtest(x, alpha)
+%SWTEST Shapiro-Wilk parametric hypothesis test of composite normality.
+%   [H, pValue, SWstatistic] = SWTEST(X, ALPHA) performs the
+%   Shapiro-Wilk test to determine if the null hypothesis of
+%   composite normality is a reasonable assumption regarding the
+%   population distribution of a random sample X. The desired significance 
+%   level, ALPHA, is an optional scalar input (default = 0.05).
+%
+%   The Shapiro-Wilk and Shapiro-Francia null hypothesis is: 
+%   "X is normal with unspecified mean and variance."
+%
+%   This is an omnibus test, and is generally considered relatively
+%   powerful against a variety of alternatives.
+%   Shapiro-Wilk test is better than the Shapiro-Francia test for
+%   Platykurtic sample. Conversely, Shapiro-Francia test is better than the
+%   Shapiro-Wilk test for Leptokurtic samples.
+%
+%   When the series 'X' is Leptokurtic, SWTEST performs the Shapiro-Francia
+%   test, else (series 'X' is Platykurtic) SWTEST performs the
+%   Shapiro-Wilk test.
+% 
+%    [H, pValue, SWstatistic] = SWTEST(X, ALPHA)
+%
+% Inputs:
+%   X - a vector of deviates from an unknown distribution. The observation
+%     number must exceed 3 and less than 5000.
+%
+% Optional inputs:
+%   ALPHA - The significance level for the test (default = 0.05).
+%  
+% Outputs:
+%  SWstatistic - The test statistic (non normalized).
+%
+%   pValue - is the p-value, or the probability of observing the given
+%     result by chance given that the null hypothesis is true. Small values
+%     of pValue cast doubt on the validity of the null hypothesis.
+%
+%     H = 0 => Do not reject the null hypothesis at significance level ALPHA.
+%     H = 1 => Reject the null hypothesis at significance level ALPHA.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                Copyright (c) 17 March 2009 by Ahmed Ben Sada          %
+%                 Department of Finance, IHEC Sousse - Tunisia           %
+%                       Email: ahmedbensaida@yahoo.com                   %
+%                    $ Revision 3.0 $ Date: 18 Juin 2014 $               %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+%
+% References:
+%
+% - Royston P. "Remark AS R94", Applied Statistics (1995), Vol. 44,
+%   No. 4, pp. 547-551.
+%   AS R94 -- calculates Shapiro-Wilk normality test and P-value
+%   for sample sizes 3 <= n <= 5000. Handles censored or uncensored data.
+%   Corrects AS 181, which was found to be inaccurate for n > 50.
+%   Subroutine can be found at: http://lib.stat.cmu.edu/apstat/R94
+%
+% - Royston P. "A pocket-calculator algorithm for the Shapiro-Francia test
+%   for non-normality: An application to medicine", Statistics in Medecine
+%   (1993a), Vol. 12, pp. 181-184.
+%
+% - Royston P. "A Toolkit for Testing Non-Normality in Complete and
+%   Censored Samples", Journal of the Royal Statistical Society Series D
+%   (1993b), Vol. 42, No. 1, pp. 37-43.
+%
+% - Royston P. "Approximating the Shapiro-Wilk W-test for non-normality",
+%   Statistics and Computing (1992), Vol. 2, pp. 117-119.
+%
+% - Royston P. "An Extension of Shapiro and Wilk's W Test for Normality
+%   to Large Samples", Journal of the Royal Statistical Society Series C
+%   (1982a), Vol. 31, No. 2, pp. 115-124.
+%
+
+%
+% Ensure the sample data is a VECTOR.
+%
+
+if numel(x) == length(x)
+    x  =  x(:);               % Ensure a column vector.
+else
+    error(' Input sample ''X'' must be a vector.');
+end
+
+%
+% Remove missing observations indicated by NaN's and check sample size.
+%
+
+x  =  x(~isnan(x));
+
+if length(x) < 3
+   error(' Sample vector ''X'' must have at least 3 valid observations.');
+end
+
+if length(x) > 5000
+    warning('Shapiro-Wilk test might be inaccurate due to large sample size ( > 5000).');
+end
+
+if max(x) == min(x)
+    H = 0;
+    pValue = 0;
+    W = 0;
+    return;
+end
+
+%
+% Ensure the significance level, ALPHA, is a 
+% scalar, and set default if necessary.
+%
+
+if (nargin >= 2) && ~isempty(alpha)
+   if ~isscalar(alpha)
+      error(' Significance level ''Alpha'' must be a scalar.');
+   end
+   if (alpha <= 0 || alpha >= 1)
+      error(' Significance level ''Alpha'' must be between 0 and 1.'); 
+   end
+else
+   alpha  =  0.05;
+end
+
+% First, calculate the a's for weights as a function of the m's
+% See Royston (1992, p. 117) and Royston (1993b, p. 38) for details
+% in the approximation.
+
+x       =   sort(x); % Sort the vector X in ascending order.
+n       =   length(x);
+mtilde  =   norminv(((1:n)' - 3/8) / (n + 1/4));
+weights =   zeros(n,1); % Preallocate the weights.
+
+if kurtosis(x) > 3
+    
+    % The Shapiro-Francia test is better for leptokurtic samples.
+    
+    weights =   1/sqrt(mtilde'*mtilde) * mtilde;
+
+    %
+    % The Shapiro-Francia statistic W' is calculated to avoid excessive
+    % rounding errors for W' close to 1 (a potential problem in very
+    % large samples).
+    %
+
+    W   =   (weights' * x)^2 / ((x - mean(x))' * (x - mean(x)));
+
+    % Royston (1993a, p. 183):
+    nu      =   log(n);
+    u1      =   log(nu) - nu;
+    u2      =   log(nu) + 2/nu;
+    mu      =   -1.2725 + (1.0521 * u1);
+    sigma   =   1.0308 - (0.26758 * u2);
+
+    newSFstatistic  =   log(1 - W);
+
+    %
+    % Compute the normalized Shapiro-Francia statistic and its p-value.
+    %
+
+    NormalSFstatistic =   (newSFstatistic - mu) / sigma;
+    
+    % Computes the p-value, Royston (1993a, p. 183).
+    pValue   =   1 - normcdf(real(NormalSFstatistic), 0, 1);
+    
+else
+    
+    % The Shapiro-Wilk test is better for platykurtic samples.
+
+    c    =   1/sqrt(mtilde'*mtilde) * mtilde;
+    u    =   1/sqrt(n);
+
+    % Royston (1992, p. 117) and Royston (1993b, p. 38):
+    PolyCoef_1   =   [-2.706056 , 4.434685 , -2.071190 , -0.147981 , 0.221157 , c(n)];
+    PolyCoef_2   =   [-3.582633 , 5.682633 , -1.752461 , -0.293762 , 0.042981 , c(n-1)];
+
+    % Royston (1992, p. 118) and Royston (1993b, p. 40, Table 1)
+    PolyCoef_3   =   [-0.0006714 , 0.0250540 , -0.39978 , 0.54400];
+    PolyCoef_4   =   [-0.0020322 , 0.0627670 , -0.77857 , 1.38220];
+    PolyCoef_5   =   [0.00389150 , -0.083751 , -0.31082 , -1.5861];
+    PolyCoef_6   =   [0.00303020 , -0.082676 , -0.48030];
+
+    PolyCoef_7   =   [0.459 , -2.273];
+
+    weights(n)   =   polyval(PolyCoef_1 , u);
+    weights(1)   =   -weights(n);
+    
+    if n > 5
+        weights(n-1) =   polyval(PolyCoef_2 , u);
+        weights(2)   =   -weights(n-1);
+    
+        count  =   3;
+        phi    =   (mtilde'*mtilde - 2 * mtilde(n)^2 - 2 * mtilde(n-1)^2) / ...
+                (1 - 2 * weights(n)^2 - 2 * weights(n-1)^2);
+    else
+        count  =   2;
+        phi    =   (mtilde'*mtilde - 2 * mtilde(n)^2) / ...
+                (1 - 2 * weights(n)^2);
+    end
+        
+    % Special attention when n = 3 (this is a special case).
+    if n == 3
+        % Royston (1992, p. 117)
+        weights(1)  =   1/sqrt(2);
+        weights(n)  =   -weights(1);
+        phi = 1;
+    end
+
+    %
+    % The vector 'WEIGHTS' obtained next corresponds to the same coefficients
+    % listed by Shapiro-Wilk in their original test for small samples.
+    %
+
+    weights(count : n-count+1)  =  mtilde(count : n-count+1) / sqrt(phi);
+
+    %
+    % The Shapiro-Wilk statistic W is calculated to avoid excessive rounding
+    % errors for W close to 1 (a potential problem in very large samples).
+    %
+
+    W   =   (weights' * x) ^2 / ((x - mean(x))' * (x - mean(x)));
+
+    %
+    % Calculate the normalized W and its significance level (exact for
+    % n = 3). Royston (1992, p. 118) and Royston (1993b, p. 40, Table 1).
+    %
+
+    newn    =   log(n);
+
+    if (n >= 4) && (n <= 11)
+    
+        mu      =   polyval(PolyCoef_3 , n);
+        sigma   =   exp(polyval(PolyCoef_4 , n));    
+        gam     =   polyval(PolyCoef_7 , n);
+    
+        newSWstatistic  =   -log(gam-log(1-W));
+    
+    elseif n > 11
+    
+        mu      =   polyval(PolyCoef_5 , newn);
+        sigma   =   exp(polyval(PolyCoef_6 , newn));
+    
+        newSWstatistic  =   log(1 - W);
+    
+    elseif n == 3
+        mu      =   0;
+        sigma   =   1;
+        newSWstatistic  =   0;
+    end
+
+    %
+    % Compute the normalized Shapiro-Wilk statistic and its p-value.
+    %
+
+    NormalSWstatistic   =   (newSWstatistic - mu) / sigma;
+    
+    % NormalSWstatistic is referred to the upper tail of N(0,1),
+    % Royston (1992, p. 119).
+    pValue       =   1 - normcdf(NormalSWstatistic, 0, 1);
+    
+    % Special attention when n = 3 (this is a special case).
+    if n == 3
+        pValue  =   6/pi * (asin(sqrt(W)) - asin(sqrt(3/4)));
+        % Royston (1982a, p. 121)
+    end
+    
+end
+
+%
+% To maintain consistency with existing Statistics Toolbox hypothesis
+% tests, returning 'H = 0' implies that we 'Do not reject the null 
+% hypothesis at the significance level of alpha' and 'H = 1' implies 
+% that we 'Reject the null hypothesis at significance level of alpha.'
+%
+
+H  = (alpha >= pValue);

wide vs narrow/single_point_test/test_lambda.m

@@ -0,0 +1,102 @@
+clc; clear; close all;
+
+N = 256;
+M = 16;
+MN = N * M;
+sigmas = [0.01];
+% noise_sigma = 0.01;
+
+filename = "./results/Signal_Model_" + string(N) + "_" + string(M) + ".mat";
+curr_metric = N;
+for TT = 1: 100
+    TT
+tau = 1e-6;
+iter_max = 25;
+epi = 0;
+[C_n, A] = get_Psi(N, M, epi);
+betas_wide = zeros(M * N, 1) + 1;
+[Lambda, Lambda_C, x] = get_sparse_vector(N, M, betas_wide, true);
+
+lg_lambdas = -3: 0.05: -1;
+lambdas = 10 .^ lg_lambdas;
+method = "cVAMPro";
+
+
+for sigma_idx = 1: length(sigmas)
+    sigma = sigmas(sigma_idx);
+    noise = get_noise(sigma, N, 1);
+    y_noise = A * x + noise;
+    
+    MSEs = zeros(length(lambdas), 1);
+    all_x_hat = zeros(length(x), length(lambdas));
+    for lambda_idx = 1: length(lambdas)
+        lambda = lambdas(lambda_idx);
+        if method == "cVAMPro"
+            [x_LASSO, x_hat_d] = cVAMPro(y_noise, A, lambda, tau, iter_max);
+        elseif method == "debiased\_LASSO"
+            [x_LASSO, sigma_w_2, threshold] = debiased_LASSO(A, y_noise, 0, sigma^2, lambda);
+        elseif method == "LASSO"
+            cvx_begin quiet
+                variable x_LASSO(MN) complex
+                minimize(lambda * norm(x_LASSO, 1) + norm(y_noise - A * x_LASSO, 2))
+            cvx_end
+        elseif method == "FISTA"
+            x_LASSO = FISTA(y_noise, A, LASSO_lambda, 1e-5);
+        elseif method == "debiased\_LASSO\_FISTA"
+            [x_LASSO, sigma_w_2, threshold] = debiased_LASSO_FISTA(A, y_noise, 0, noise_sigma^2, lambda);
+        elseif method == "BP"
+            sz = size(A);
+            N = sz(2);
+            cvx_begin quiet
+                variable x_LASSO(N) complex
+                minimize(norm(x_LASSO, 1))
+            subject to
+                A * x_LASSO == y_noise
+            cvx_end
+        end
+        MSEs(lambda_idx) = sum(real(x_LASSO - x));
+        all_x_hat(:, lambda_idx) = x_LASSO;
+    end
+    abs_MSEs = abs(MSEs);
+    [metric, idx] = min(abs_MSEs);
+    if curr_metric > metric
+        ref_lambda = lambdas(idx);
+        save(filename, "A", "C_n", "ref_lambda");
+        curr_metric = metric;
+    end
+    if curr_metric < 0.1
+        figure;
+        subplot(211);
+        semilogx(lambdas, MSEs);
+        yline(0);
+        ylim([-length(Lambda) - 0.1, 0.1]);
+        xlabel("\lambda");
+        ylabel("MSE (x\_hat - x)");
+        title("N = " + string(N) + ", M = " + string(M) + ", \sigma = " + string(sigma) + ", " + method);
+        subplot(212);
+        plot(C_n);
+        xlabel("Frequency code");
+        ylabel("Times");
+        break
+    end
+    % for i = 1: length(C_n)
+    %     fprintf("%d, ", C_n(i))
+    % end
+    % fprintf("%.15f\n", sum(MSEs));
+
+    % filename = ...
+    %     string(N) + "_" + ...
+    %     string(M) + "_" + ...
+    %     string(sigma) + "_" + ...
+    %     method + ".mat";
+    % save(filename, "lambdas", "MSEs", "x", "all_x_hat");
+    
+    % figure;
+    % semilogx(lambdas, xx);
+    % yline(0);
+    % ylim([-length(Lambda) - 0.1, 0.1]);
+    % xlabel("\lambda");
+    % ylabel("MSE (x\_hat - x)");
+    % title("N = " + string(FAR_N) + ", M = " + string(FAR_M) + ", \sigma = " + string(sigma) + ", " + method);
+end
+end

wide vs narrow/test_equal_RCS.m

@@ -0,0 +1,36 @@
+clc; clear; close all;
+
+run("../parameters.m");
+
+global f_c N_narrow N N_wide numP f_s N_high betas_wide delta_R_wide c;
+f_c = f_c * 1e3;
+
+numP = 1;
+s_T_wide = zeros(1, N * numP);
+
+for j = 1:numP
+    idx = ((j - 1) * N + (1:N_wide));
+    tp = idx / f_s;
+    s_T_wide(1, idx) = exp(1j * 2 * pi * f_c * tp);
+end
+
+s_R_wides = zeros(N_wide, N_wide * numP);
+
+for k = 1:N_wide
+    s_R_wides(k, :) = [zeros(1, k), s_T_wide(1:end - k)];
+end
+
+echo_1 = zeros(1, N * numP);
+
+for i = 1:10
+    echo_1 = echo_1 + betas_wide(i) * s_R_wides(i, :);
+end
+
+equal_beta = 0 + 0j;
+
+for i = 1:10
+    equal_beta = equal_beta + betas_wide(i) * exp(1j * 2 * pi * 2 * delta_R_wide * i / c);
+end
+
+fprintf("理论等效RCS: " + string(equal_beta) + "\n");
+fprintf("验证等效RCS: " + string(echo_1(50)) + "\n");