Update Block RIP

Block RIP/Phase transitions/can_recovery.m

@@ -0,0 +1,25 @@
+function flg = can_recovery(Psi, s, eps)
+    [m, n] = size(Psi);
+
+    x = zeros(n, 1);
+    random_indices = randperm(n, s);
+    x(random_indices) = randn(s, 1);
+    x = sign(x);
+    y = Psi * x;
+
+    cvx_begin
+    variable s1(n)
+    minimize(norm(s1, 1))
+    subject to
+    norm(y - Psi * s1) <= eps
+    cvx_end
+
+    p = norm(x - s1, 2);
+
+    if p < eps
+        flg = 1;
+    else
+        flg = 0;
+    end
+
+end

Block RIP/Phase transitions/get_Psi.m

@@ -0,0 +1,11 @@
+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

Block RIP/Phase transitions/main.m

@@ -0,0 +1,34 @@
+clc;
+clear;
+
+N = 32;
+M = 4;
+k = 5;
+epi = 0.02;
+
+Psi_wide = get_Psi(N, M, epi);
+
+max_n = 100;
+max_s = 100;
+trials_time = 50;
+d = 5;
+eps = 1e-5;
+
+prob = zeros(max_n, max_s);
+
+for n = 1:max_n
+    Psi_narrow = get_Psi(n, M, 0);
+    parfor s = 1:max_s
+        if n > s
+            x = 0;
+            for t = 1:trials_time
+                x = x + can_recovery(Psi_narrow, s, eps);
+            end
+            prob(n, s) = x / trials_time;
+        else
+            prob(n, s) = 0;
+        end
+    end
+end
+
+save("Expr1.mat", "prob");

Block RIP/RIP/get_Psi.m

@@ -0,0 +1,11 @@
+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

Block RIP/RIP/get_ric.m

@@ -0,0 +1,16 @@
+function [RICs, RIC] = get_ric(A, K)
+    sample_times = 500;
+    [m, n] = size(A);
+    RICs = zeros(2 * sample_times, 1);
+
+    for i = 1: sample_times
+        indices = randperm(n, K);
+        A_Lambda = A(:, indices(:));
+        B = conj(A_Lambda') * A_Lambda;
+        eig_B = real(eig(B));
+        RICs(2 * i - 1) = min(eig_B) - 1;
+        RICs(2 * i) = max(eig_B) - 1;
+    end
+
+    RIC = max(RICs);
+end

Block RIP/RIP/main.m

@@ -0,0 +1,26 @@
+clc;
+clear;
+
+N = 32;
+M = 4;
+k = 5;
+epi = 0.02;
+
+Psi_narrow = get_Psi(N, M, 0);
+Psi_wide = get_Psi(N, M, epi);
+
+fprintf("RIC_narrow: ");
+[RICs_narrow, RIC_narrow] = get_ric(Psi_narrow, k);
+fprintf("\n\n");
+
+fprintf("RIC_wide: ");
+[RICs_wide, RIC_wide] = get_ric(Psi_wide, k);
+fprintf("\n\n");
+
+RICs_narrow = sort(RICs_narrow);
+RICs_wide = sort(RICs_wide);
+
+plot(1:length(RICs_narrow), RICs_narrow, 1:length(RICs_wide), RICs_wide);
+xlabel("Experiment times");
+ylabel("Block RIC");
+title("RIC");

Block RIP/bound/get_our_ub.m

@@ -0,0 +1,3 @@
+function result = get_our_ub(delta, N, C, M, eps)
+    result = delta^2 * N / (C * M * (log(M) - log(eps)));
+end

Block RIP/bound/get_wl_ub.m

@@ -0,0 +1,13 @@
+function result = get_wl_ub(N, M, eps)
+    k1 = sqrt((M - 1) / N);
+    k2 = log(M * N);
+    k3 = sqrt(k2 - log(eps));
+    k4 = sqrt(log(2 * M) - log(eps));
+    
+    delta_1 = 24 * k1 * k2 * (2 * k3 + 1);
+    delta_2 = 1.5 * k1 * (2 * k4 + 1);
+
+    numerator = N * (1/8 - delta_1 - delta_2) ^ 2;
+    denominator = 81 * M * k2 * (1 + 2 * delta_2 / 3);
+    result = numerator / denominator;
+end

Block RIP/bound/main.m

@@ -0,0 +1,25 @@
+clc;
+clear;
+
+M = 50;
+eps = 1e-5;
+delta = 0.3;
+C = 10.66;
+
+nums = 10:0.1:50;
+num = length(nums);
+N = zeros(1, num);
+upper_bound_C1 = zeros(1, num);
+upper_bound_C2 = zeros(1, num);
+
+for i = 1: length(nums)
+    N(i) = 10^nums(i);
+    upper_bound_C1(i) = get_our_ub(1, N(i), C, M, eps);
+    upper_bound_C2(i) = get_wl_ub(N(i), M, eps);
+    N(i) = nums(i);
+end
+
+semilogy(N, upper_bound_C1, 'r-', N, upper_bound_C2, 'b');
+xlabel("log_{10}(N)");
+ylabel("Upper bound of K");
+title("Prob(\delta_K < 0.3) >= 1 - 10^{-5}")

Block RIP/config_parameters.m

@@ -1,34 +0,0 @@
-N = 128; % 脉冲个数
-M = 1; % 频点个数
-K = 10; % 目标个数
-d = 32;
-epsilon = 1e-5; % 误差
-
-f_c = 10e9; % 初始载频 10GHz
-Delta_f = 8e6; % 载频步进间隔 8MHz
-% c = 299792458; % 光速
-c = 3e8;
-scatter_coef = 0.3; % 目标散射强度
-
-B = 64e6; % 带宽 64MHz
-% B_0 = 1e9;
-
-f_s = 3 * f_c; % 快时间采样率
-T_p = 10 / f_s; % 单载频脉冲下的采样周期 / 脉冲宽度
-T_r = T_p * 10;
-
-r_0 = 1; % 初始距离 r(0)
-velocity = 5e5; % 目标速度（假设目标做匀速直线运动）
-
-lambda = c / f_c; % 雷达工作波长
-
-% 仿真时间
-delta_t = 1e-3 * T_p;
-max_t = 30 * T_r;
-range_t = 0:delta_t:max_t-delta_t;
-len = length(range_t);
-
-% 绘图
-figure_flag_1 = false;
-figure_flag_2 = false;
-figure_flag_3 = false;

Block RIP/draw.m

@@ -1,7 +0,0 @@
-filename = "/Users/ksyer/CLionProjects/BlockRIP/cmake-build-debug/1.txt";
-
-df = dlmread(filename);
-x = df(:, 1);
-y = df(:, 2);
-
-scatter(x, y);

Block RIP/main.m

@@ -1,21 +0,0 @@
-% 
-M = 200;
-N = 1e12;
-P = N / M;
-s = 10;
-epsilon = 1e-5;
-
-x = 22:0.2:30;
-N = zeros(size(x));
-P = zeros(size(x));
-sigma = zeros(size(x));
-
-for i = 1: length(x)
-    N(i) = 10^x(i);
-    P(i) = N(i) / M;
-    ita_1_lb = sqrt(172.24 * 32.0 * s * (log(4 * s) ^ 2) * log(8 * N(i)) * log(9 * P(i)) / P(i));
-    ita_2_lb = sqrt(32.0 / 3.0 * s * (-log(epsilon)) / P(i));
-    sigma(i) = ita_1_lb * (1 + ita_1_lb) + ita_2_lb;
-end
-
-semilogx(P, sigma);

Block RIP/main2.m

@@ -1,87 +0,0 @@
-% 仿真入口，确定 N 和 s 后，通过遍历 eta_1 和 eta_2 来计算出 P
-
-% 设置用于遍历 eta 的参数
-eta_1_step = 1e-2;
-eta_2_step = 1e-8;
-
-eta_1_start = eta_1_step;
-eta_2_start = eta_2_step;
-
-eta_1_end = (sqrt(5) - 1) / 2; % 大于这个值时，eta_1 ^2 + eta_1 必定会大于 1
-eta_2_end = 1e-3;
-
-eta_1 = eta_1_start:eta_1_step:eta_1_end;
-eta_2 = eta_2_start:eta_2_step:eta_2_end;
-
-eta_1_N = length(eta_1);
-eta_2_N = length(eta_2);
-
-N = 1e14;
-s = 10;
-epsilon = 1e-5;
-C_1 = 5576;
-C_2 = 10.66;
-
-result = zeros(eta_1_N, eta_2_N);
-metric_1 = zeros(eta_1_N, 1);
-metric_2 = zeros(eta_2_N, 1);
-
-for i = 1:eta_1_N
-    metric_1(i) = (C_1 * s * (log(4*s))^2 * log(8*N)) / (eta_1(i) ^ 2);
-end
-
-for j = 1:eta_2_N
-    metric_2(j) = (C_2 * s * log(epsilon^-1)) / (eta_2(j) ^ 2);
-end
-
-
-sigmas = zeros(eta_1_N, 1 * eta_2_N);
-results = zeros(eta_1_N, 1 * eta_2_N);
-k = 0;
-d = 1e2;
-for i = 1:eta_1_N
-    for j = 1:eta_2_N
-        sigma = eta_1(i) * (1 + eta_1(i)) + eta_2(j);
-        if sigma > 1
-            continue
-        end
-        sigmas(i, j) = sigma;
-        k = k + 1;
-        m1 = metric_1(i);
-        m2 = metric_2(j);
-
-        m = solve_test(m1);
-        if m1 < m2
-            m = max(m, m2);
-        end
-
-        if m > N
-            results(i, j) = 0;
-        else
-            results(i, j) = N / m;
-        end
-    end
-end
-
-figure;
-% semilogy(sigmas, results);
-h = heatmap(results);
-h.GridVisible = false;
-ax = gca;
-xn = length(ax.XDisplayLabels);
-yn = length(ax.YDisplayLabels);
-
-for i = 1:length(ax.XDisplayLabels)
-    % if rem(i, rem(xn, 10)) ~= 0
-    %     ax.XDisplayLabels(i) = {nan};
-    % end
-end
-
-for i = 1:length(ax.YDisplayLabels)
-    % if rem(i, rem(yn, 10)) ~= 0
-    %     ax.YDisplayLabels(i) = {nan};
-    % end
-end
-
-
-% heatmap(result);

Block RIP/solve_test.m

@@ -1,20 +0,0 @@
-% 由于函数单调，因此可以用二分法求 x/log(9x) = k 的解
-
-function result = solve_test(k)
-    eps = 1e-3;
-    l = 1;
-    r = k;
-    while r - l > eps
-        mid = (l+r) / 2;
-        if (foo(mid) < foo(r))
-            l = mid;
-        else
-            r = mid;
-        end
-    end
-    result = l;
-end
-
-function f = foo(x)
-    f = x / log(9 * x);
-end