Update BlockRIP experiment

Block RIP/config_parameters.m

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

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

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

@@ -0,0 +1,87 @@
+% 仿真入口，确定 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

@@ -0,0 +1,20 @@
+% 由于函数单调，因此可以用二分法求 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