diff --git a/Block RIP/Phase transitions/Expr1.mat b/Block RIP/Phase transitions/Expr1.mat new file mode 100644 index 0000000..1e88895 Binary files /dev/null and b/Block RIP/Phase transitions/Expr1.mat differ diff --git a/Block RIP/Phase transitions/can_recovery.m b/Block RIP/Phase transitions/can_recovery.m new file mode 100644 index 0000000..9bf1f99 --- /dev/null +++ b/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 diff --git a/Block RIP/Phase transitions/get_Psi.m b/Block RIP/Phase transitions/get_Psi.m new file mode 100644 index 0000000..3042b5d --- /dev/null +++ b/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 diff --git a/Block RIP/Phase transitions/main.m b/Block RIP/Phase transitions/main.m new file mode 100644 index 0000000..fb13276 --- /dev/null +++ b/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"); diff --git a/Block RIP/RIP/FARblockepsilon2.mat b/Block RIP/RIP/FARblockepsilon2.mat new file mode 100644 index 0000000..35dcfd6 Binary files /dev/null and b/Block RIP/RIP/FARblockepsilon2.mat differ diff --git a/Block RIP/RIP/get_Psi.m b/Block RIP/RIP/get_Psi.m new file mode 100644 index 0000000..3042b5d --- /dev/null +++ b/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 diff --git a/Block RIP/RIP/get_ric.m b/Block RIP/RIP/get_ric.m new file mode 100644 index 0000000..8b2f454 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Block RIP/RIP/main.m b/Block RIP/RIP/main.m new file mode 100644 index 0000000..616f6b4 --- /dev/null +++ b/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"); \ No newline at end of file diff --git a/Block RIP/bound/get_our_ub.m b/Block RIP/bound/get_our_ub.m new file mode 100644 index 0000000..27105eb --- /dev/null +++ b/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 diff --git a/Block RIP/bound/get_wl_ub.m b/Block RIP/bound/get_wl_ub.m new file mode 100644 index 0000000..7810f54 --- /dev/null +++ b/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 diff --git a/Block RIP/bound/main.m b/Block RIP/bound/main.m new file mode 100644 index 0000000..6de5259 --- /dev/null +++ b/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}") \ No newline at end of file diff --git a/Block RIP/config_parameters.m b/Block RIP/config_parameters.m deleted file mode 100644 index 8dfd365..0000000 --- a/Block RIP/config_parameters.m +++ /dev/null @@ -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; diff --git a/Block RIP/draw.m b/Block RIP/draw.m deleted file mode 100644 index ffb6c28..0000000 --- a/Block RIP/draw.m +++ /dev/null @@ -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); diff --git a/Block RIP/main.m b/Block RIP/main.m deleted file mode 100644 index a2a4097..0000000 --- a/Block RIP/main.m +++ /dev/null @@ -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); diff --git a/Block RIP/main2.m b/Block RIP/main2.m deleted file mode 100644 index b0d51ec..0000000 --- a/Block RIP/main2.m +++ /dev/null @@ -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); diff --git a/Block RIP/solve_test.m b/Block RIP/solve_test.m deleted file mode 100644 index bdf23e6..0000000 --- a/Block RIP/solve_test.m +++ /dev/null @@ -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