diff --git a/Phase Transitions in FAR using CS/Expr1.m b/Phase Transitions in FAR using CS/Expr1.m
new file mode 100644
index 0000000..2bd5a61
--- /dev/null
+++ b/Phase Transitions in FAR using CS/Expr1.m	
@@ -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");
diff --git a/Expr 1.mat b/Phase Transitions in FAR using CS/Expr1.mat
similarity index 100%
rename from Expr 1.mat
rename to Phase Transitions in FAR using CS/Expr1.mat
diff --git a/can_recovery.m b/Phase Transitions in FAR using CS/Expr1_can_recovery.m
similarity index 50%
rename from can_recovery.m
rename to Phase Transitions in FAR using CS/Expr1_can_recovery.m
index 84765f2..a0e0fc9 100644
--- a/can_recovery.m
+++ b/Phase Transitions in FAR using CS/Expr1_can_recovery.m	
@@ -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
diff --git a/Phase Transitions in FAR using CS/Expr1_draw.m b/Phase Transitions in FAR using CS/Expr1_draw.m
new file mode 100644
index 0000000..4587798
--- /dev/null
+++ b/Phase Transitions in FAR using CS/Expr1_draw.m	
@@ -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)')
diff --git a/Phase Transitions in FAR using CS/Expr2.m b/Phase Transitions in FAR using CS/Expr2.m
new file mode 100644
index 0000000..4249730
--- /dev/null
+++ b/Phase Transitions in FAR using CS/Expr2.m	
@@ -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");
diff --git a/Phase Transitions in FAR using CS/Expr2_can_recovery.m b/Phase Transitions in FAR using CS/Expr2_can_recovery.m
new file mode 100644
index 0000000..c7df676
--- /dev/null
+++ b/Phase Transitions in FAR using CS/Expr2_can_recovery.m	
@@ -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
diff --git a/Phase Transitions in FAR using CS/Expr3.asv b/Phase Transitions in FAR using CS/Expr3.asv
new file mode 100644
index 0000000..cc6fbaa
--- /dev/null
+++ b/Phase Transitions in FAR using CS/Expr3.asv	
@@ -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");
diff --git a/Phase Transitions in FAR using CS/Expr3.m b/Phase Transitions in FAR using CS/Expr3.m
new file mode 100644
index 0000000..31d44a5
--- /dev/null
+++ b/Phase Transitions in FAR using CS/Expr3.m	
@@ -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");
diff --git a/Phase Transitions in FAR using CS/Expr3.mat b/Phase Transitions in FAR using CS/Expr3.mat
new file mode 100644
index 0000000..b9e7934
Binary files /dev/null and b/Phase Transitions in FAR using CS/Expr3.mat differ
diff --git a/Phase Transitions in FAR using CS/Expr3_can_recovery.asv b/Phase Transitions in FAR using CS/Expr3_can_recovery.asv
new file mode 100644
index 0000000..0e8265e
--- /dev/null
+++ b/Phase Transitions in FAR using CS/Expr3_can_recovery.asv	
@@ -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
diff --git a/Phase Transitions in FAR using CS/Expr3_can_recovery.m b/Phase Transitions in FAR using CS/Expr3_can_recovery.m
new file mode 100644
index 0000000..b67b7e7
--- /dev/null
+++ b/Phase Transitions in FAR using CS/Expr3_can_recovery.m	
@@ -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
diff --git a/Phase Transitions in FAR using CS/I.mat b/Phase Transitions in FAR using CS/I.mat
new file mode 100644
index 0000000..d46c33b
Binary files /dev/null and b/Phase Transitions in FAR using CS/I.mat differ
diff --git a/LICENSE b/Phase Transitions in FAR using CS/LICENSE
similarity index 100%
rename from LICENSE
rename to Phase Transitions in FAR using CS/LICENSE
diff --git a/Phase Transitions in FAR using CS/README.md b/Phase Transitions in FAR using CS/README.md
new file mode 100644
index 0000000..359c146
--- /dev/null
+++ b/Phase Transitions in FAR using CS/README.md	
@@ -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`:
diff --git a/Phase Transitions in FAR using CS/calc_integral.m b/Phase Transitions in FAR using CS/calc_integral.m
new file mode 100644
index 0000000..1d99575
--- /dev/null
+++ b/Phase Transitions in FAR using CS/calc_integral.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
diff --git a/Phase Transitions in FAR using CS/phase_transition.m b/Phase Transitions in FAR using CS/phase_transition.m
new file mode 100644
index 0000000..94c3ce4
--- /dev/null
+++ b/Phase Transitions in FAR using CS/phase_transition.m	
@@ -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
diff --git a/get_recovery_prob.m b/get_recovery_prob.m
deleted file mode 100644
index 73b38ce..0000000
--- a/get_recovery_prob.m
+++ /dev/null
@@ -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
\ No newline at end of file
diff --git a/main.m b/main.m
deleted file mode 100644
index 5a43f1d..0000000
--- a/main.m
+++ /dev/null
@@ -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
\ No newline at end of file