Add Phase Transitions in FAR using CS folder

This commit is contained in:
Ksyer
2023-12-19 17:12:26 +08:00
parent 4881771558
commit 4a7a8ab31b
18 changed files with 390 additions and 28 deletions
+26
View File
@@ -0,0 +1,26 @@
% Expr1.m to draw Fig 2(a)
max_n = 125;
max_s = 35;
trials_time = 50;
m = 4;
d = 32;
eps = 1e-5;
prob = zeros(max_n, max_s);
for n = 1:max_n
parfor s = 1:max_s
x = 0;
for t = 1:trials_time
x = x + Expr1_can_recovery(n, s, m, d, eps);
end
prob(n, s) = x / trials_time;
end
end
save("Expr1.mat", "prob");
@@ -1,8 +1,6 @@
function flg = can_recovery(n, s)
d = 100;
eps = 1e-5;
theta = randn(n, d);
function flg = Expr1_can_recovery(n, s, m, d, eps)
theta = randn(n, m * d);
x = zeros(d, 1);
random_indices = randperm(d, s);
x(random_indices) = randn(s, 1);
@@ -10,16 +8,18 @@ function flg = can_recovery(n, s)
y = theta * x;
cvx_begin
variable s1(d)
minimize(norm(s1, 1))
subject to
norm(y - theta * s1) <= eps
variable s1(d)
minimize(norm(s1, 1))
subject to
norm(y - theta * s1) <= eps
cvx_end
p = norm(x-s1, 2);
p = norm(x - s1, 2);
if p < eps
flg = 1;
else
flg = 0;
end
end
@@ -0,0 +1,11 @@
load 'Expr 1.mat'
p = rot90(prob, 3);
imagesc(p)
set(gca, 'xtick', 1:10:100);
set(gca, 'xticklabel', {'0', '10', '20', '30', '40', '50', '60', '70', '80', '90', '100'})
xlabel('Sparsity of signal (s)')
set(gca, 'ytick', 1:10:100);
set(gca, 'yticklabel', {'0', '10', '20', '30', '40', '50', '60', '70', '80', '90', '100'})
ylabel('Number of Mersurements (n)')
+25
View File
@@ -0,0 +1,25 @@
% Expr2.m to draw Fig 2(b)
max_n = 100;
max_s = 100;
trials_time = 50;
d = 100;
eps = 1e-5;
prob = zeros(max_n, max_s);
for n = 1:max_n
for s = 1:max_s
x = 0;
for t = 1:trials_time
x = x + Expr2_can_recovery(n, s, d, eps);
end
prob(n, s) = x / trials_time;
end
end
save("Expr2.mat", "prob");
@@ -0,0 +1,36 @@
function flg = Expr2_can_recovery(n, s, d, eps)
theta_real = randn(n, d);
theta_imag = randn(n, d);
theta = theta_real + 1j * theta_imag;
x = zeros(d, 1);
random_indices = randperm(d, s);
x(random_indices) = randn(s, 1);
for i = 1:d
if x(i) ~= 0
t = rand(1);
x(i) = cos(t) + 1j * sin(t);
end
end
y = theta * x;
cvx_begin
variable s1(d)
minimize(norm(s1, 1))
subject to
norm(y - theta * s1) <= eps
cvx_end
p = norm(x - s1, 2);
if p < eps
flg = 1;
else
flg = 0;
end
end
@@ -0,0 +1,72 @@
% Expr2.m to draw Fig 3(a)
clc;
clear;
M = 4;
N = 128;
rate = 0.02; % \Delta f / f_c = 0.02
beta = 1;
max_n = 125;
max_k = 25;
trials_time = 50;
prob = zeros(max_n, max_k);
rng(1)
Phi_far = zeros(N, N * M);
C = randi([1, M], 1, N);
xi = 1 + 0.02 * C;
% for i = 1:N
%
% Phi_x = zeros(N, M);
%
% for p = 1:M
%
% for q = 1:N
% range = ((p - 1) * C(i)) / M; % p start from 0 or 1?
% doppler = ((q - 1) * (i - 1) * xi(i)) / N; % q start from 0 or 1?
% Phi_x(q, p) = exp(1j * 2 * pi * (range + doppler)) / sqrt(N);
% end
%
% end
%
% % t == ctranspose(Phi_x(:,1)) * Phi_x(:,1) == 1
% Phi_far(:, (i - 1) * M + 1:i * M) = Phi_x;
%
% end
for i = 1:n
for j = 1:m
n = 60;
k = 5;
s = beta * k * M;
x = 0;
for t = 1:trials_time
x = x + Expr3_can_recovery(Phi_far, N, M, n, s, eps);
end
disp(x);
for n = 1:max_n
parfor k = 1:max_k
s = beta * k * M;
x = 0;
for t = 1:trials_time
x = x + Expr3_can_recovery(Phi_far, N, M, n, s, eps);
end
prob(n, k) = x / trials_time;
end
end
save("Expr3.mat", "prob");
+63
View File
@@ -0,0 +1,63 @@
% Expr2.m to draw Fig 3(a)
clc;
clear;
M = 4;
N = 128;
rate = 0.02; % \Delta f / f_c = 0.02
beta = 1;
max_n = 125;
max_k = 25;
trials_time = 10;
eps = 1e-5;
prob = zeros(max_n, max_k);
rng(1)
Phi_far = zeros(N, N * M);
C = randi([0, M-1], 1, N);
vecn = 0:N-1;
xi = 1 + rate .* C;
for i = 1:N
for j = 1:M
t = (i-1)*M+j;
temp = 1/sqrt(N) .* exp(1j * 2 * pi + (j-1) / M .* C + 1j * 2 * pi * (i-1)/N .* vecn .* xi);
Phi_far(:, t) = temp;
end
disp(ctranspose(Phi_far(:,1)) * Phi_far(:,1))
end
% n = 60; 0
n = 70;
k = 7;
s = beta * k * M;
x = 0;
for t = 1:trials_time
x = x + Expr3_can_recovery(Phi_far, N, M, n, s, eps);
end
disp(x);
for n = 1:max_n
parfor k = 1:max_k
s = beta * k * M;
x = 0;
for t = 1:trials_time
x = x + Expr3_can_recovery(Phi_far, N, M, n, s, eps);
end
prob(n, k) = x / trials_time;
end
end
save("Expr3.mat", "prob");
Binary file not shown.
@@ -0,0 +1,43 @@
function flg = Expr3_can_recovery(Phi_far, N, M, n, s, eps)
Phi = zeros(n, N * M);
indices = randperm(N, n);
for i = 1:n
Phi(i, :) = Phi_far(indices(i), :);
end
d = M*N;
x = zeros(d, 1);
random_indices = randperm(d, s);
x(random_indices) = randn(s, 1);
x = sign(x);
for i = 1:M * N
if x(i) ~= 0
t = rand(1) * 2 * pi;
x(i) = cos(t) + 1j * sin(t);
end
end
y = Phi * x;
cvx_begin quiet
variable s1(d) complex
minimize(norm(s1, 1))
subject to
norm(y - Phi * s1) <= eps
cvx_end
x = sign(real(x));
s1 = real(s1);
p = norm(x - s1, 2);
disp(p);
if p < eps
flg = 1;
else
flg = 0;
end
end
@@ -0,0 +1,44 @@
function flg = Expr3_can_recovery(Phi_far, N, M, n, s, eps)
Phi = zeros(n, N * M);
indices = randperm(N, n);
for i = 1:n
Phi(i, :) = Phi_far(indices(i), :);
end
d = M*N;
x = zeros(d, 1);
random_indices = randperm(d, s);
x(random_indices) = randn(s, 1);
x = sign(x);
for i = 1:M * N
if x(i) ~= 0
t = rand(1) * 2 * pi;
x(i) = cos(t) + 1j * sin(t);
end
end
y = Phi * x;
cvx_begin quiet
variable s1(d) complex
minimize(norm(s1, 1))
subject to
norm(y - Phi * s1) <= eps
cvx_end
x = sign(real(x));
s1 = real(s1);
p = norm(x - s1, 2);
disp(p);
if p < eps
flg = 1;
else
flg = 0;
end
end
Binary file not shown.
@@ -0,0 +1,8 @@
Codes:
- `calc_int.m`: Calc integral $f(x) = (x-tau)^2 * sqrt(2/pi) * exp(-x^2/2)$, range: $[tau, +inf]$.
- `can_recovery.m`: Determine if the signal can be recovered.
- `draw.m`: Draw.
- `get_recovery_prob.m`: Obtain the probability that the signal can be recovered under T experiments.
- `main.m`: Main function entry.
- `phase_transition.m`:
@@ -0,0 +1,15 @@
function I = calc_integral(tau, m)
if (~exist('m', 'var'))
m = 1;
end
if tau < 0
I = 0;
else
syms x f;
f = (x - tau) ^ 2 * exp(-x ^ 2/2) * u ^ (m - 1) / (2 ^ (m / 2 - 1) * gamma(m / 2));
I = double(int(f, [tau, +inf]));
end
end
@@ -0,0 +1,37 @@
% Eq 20
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 = 100;
s_b_range = s_b_min:s_b_max;
result = zeros(32, 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_integral(tau);
end
end
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
result(s_b_idx) = min(f_set);
end
-7
View File
@@ -1,7 +0,0 @@
function x = get_recovery_prob(n, s, N)
x = 0;
for t = 1: N
x = x + can_recovery(n, s);
end
x = x / N;
end
-11
View File
@@ -1,11 +0,0 @@
N = 100;
M = 100;
T = 50;
prob = zeros(N, M);
for n = 1: N
parfor s = 1: M
prob(n, s) = get_recovery_prob(n, s, T);
end
end
c