Add Radar simulation

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

Phase Transitions in FAR using CS/Expr3.m

@@ -1,4 +1,4 @@
-% Expr2.m to draw Fig 3(a)
+% Expr3.m to draw Fig 3(a)
 clc;
 clear;
 
@@ -8,7 +8,7 @@ epi = 0.02; % \Delta f / f_c = 0.02
 beta = 1;
 
 max_n = 125;
-max_k = 25;
+max_k = 3;
 trials_time = 5;
 eps = 1e-5;
 
@@ -33,7 +33,7 @@ end
 
 for n = 1:max_n
 
-    parfor k = 1:max_k
+    for k = 1:max_k
         s = beta * k * M;
         x = 0;
 

Phase Transitions in FAR using CS/Expr3_can_recovery.m

@@ -13,13 +13,17 @@ function flg = Expr3_can_recovery(Phi_far, N, M, n, s, eps)
     y = Phi * sparse_signal(:);
 
     cvx_begin quiet
-    variable x(M * N) complex
-    minimize(norm(x, 1))
+    variable x(M, N) complex
+        norm21 = 0;
+        for i = 1:N
+            norm21 = norm21 + norm(x(:, i));
+        end
+        minimize(norm21)
     subject to
-    Phi * x == y
+        Phi * x(:) == y
     cvx_end
 
-    p = norm(x - sparse_signal(:), 2);
+    p = norm(x(:) - sparse_signal(:), 2);
 
     if p < eps
         flg = 1;

Phase Transitions in FAR using CS/Lyh/捷变频相变/FAR_block_sparsity_phase1.m

@@ -0,0 +1,48 @@
+close all;
+clear all;
+clc;
+M = 4;
+N = 128;
+%block_sparsity = 1;
+tol = 1e-5;
+trial = 50;
+epi = 0.02;
+result = zeros(N,25);
+for col = 4:4:128
+    for block_sparsity = 10:18
+        success_count = 0;
+        for loop = 1:trial
+            FAR_model = zeros(N,M*N);
+            %Cn = randperm(M)-1
+            for n = 0 : N-1
+                Cn = floor(rand()*M);
+                for q = 0 : N-1
+                    for p = 0:M-1
+                        FAR_model(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi));
+                    end
+                end
+            end
+            col_choose = randperm(N,col);
+            FAR_model = FAR_model(col_choose,:);
+            sparse_signal = zeros(M,N);
+            block = randperm(N,block_sparsity);
+            sparse_signal(:,block) = exp(1i*2*pi*rand(M,block_sparsity));
+            y = FAR_model * sparse_signal(:);
+            cvx_begin
+            variable x(M,N) complex
+            norm21 = 0;
+            for i = 1:N
+                norm21 = norm21 + norm(x(:,i));
+            end
+            minimize(norm21)
+            subject to
+            FAR_model * x(:) == y
+            cvx_end
+            if norm(x(:)-sparse_signal(:))<tol
+                success_count = success_count+1;
+            end
+        end
+        result(col,block_sparsity) = success_count/trial;
+    end
+end
+save('FARblockepsilon2.mat');

Phase Transitions in FAR using CS/Lyh/捷变频相变/FAR_sparsity_phase.m

@@ -0,0 +1,47 @@
+close all;
+clear all;
+clc;
+M = 4;
+N = 128;
+%block_sparsity = 1;
+tol = 1e-5;
+trial = 50;
+epi = 0.02;
+result = zeros(N,25);
+for col = 4:4:128
+    for block_sparsity = 1:25
+        success_count = 0;
+        for loop = 1:trial
+            FAR_model = zeros(N,M*N);
+            %Cn = randperm(M)-1
+            for n = 0 : N-1
+                Cn = floor(rand()*M);
+                for q = 0 : N-1
+                    for p = 0:M-1
+                        FAR_model(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi));
+                    end
+                end
+            end
+            
+            col_choose = randperm(N,col);
+            FAR_model = FAR_model(col_choose,:);
+            sparse_signal = zeros(M,N);
+
+            block = randperm(N,block_sparsity);
+            sparse_signal(:,block) = exp(1i*2*pi*rand(M,block_sparsity));
+            
+            y = FAR_model * sparse_signal(:);
+            cvx_begin
+            variable x(M*N) complex
+            minimize(norm(x,1))
+            subject to
+            FAR_model * x == y
+            cvx_end
+            if norm(x-sparse_signal(:))<tol
+                success_count = success_count+1;
+            end
+        end
+        result(col,block_sparsity) = success_count/trial;
+    end
+end
+save('FARepsilon.mat');

Phase Transitions in FAR using CS/Lyh/捷变频相变/theoretic.m

@@ -0,0 +1,7 @@
+function n = theoretic(m,s,d)
+    syms t;
+    syms u;
+    f = s*(m+t^2)+(d-s)*int((u-t)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t,inf);
+    g = diff(f,t);
+    t1 = solve(g);
+    n = s*(m+t1^2)+(d-s)*int((u-t1)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t1,inf);

Phase Transitions in FAR using CS/Lyh/高斯相变/gauss_phase_real.m

@@ -0,0 +1,27 @@
+N = 100;
+gauss_phase_res = zeros(100,100);
+for col = 1:100
+    %¾ØÕóÉú³É
+    for p =1:100
+        suc = 0;
+        for loop = 1:50
+            x1 = zeros(N,1);
+            q = randperm(N,p);
+            x1(q) = randn(p,1);
+            fai = randn(col,N);
+            b = fai*x1;
+            cvx_begin quiet
+                variable x(N)
+                minimize( norm( x, 1 ) )
+                subject to
+                    fai * x == b
+            cvx_end
+            %disp((norm(x-x1,1)))
+            if (norm(x-x1,1))<10e-5
+                 suc = suc+1;
+             end
+        end
+        gauss_phase_res(col,p)=suc/50;
+    end
+end
+save gauss_phase_real;

Radar simulation/R_d_func.m

@@ -0,0 +1,13 @@
+function y = R_d_func(t)
+    global c f_n T_r;
+
+    ti = t - (2/c) * r(t);
+    p = rem(ti, T_r); % p = t - nT_r
+    n = round((ti - p) / T_r);
+    
+    if (n < 0)
+        n = 0;
+    end
+
+    y = R_x_func(t) * exp(1j * -2 * pi * f_n(n + 1) * (t - n * T_r));
+end

Radar simulation/R_x_func.m

@@ -0,0 +1,5 @@
+function y = R_x_func(t)
+    global scatter_coef c;
+    ti = t - (2/c) * r(t);
+    y = scatter_coef * T_x_func(ti);
+end

Radar simulation/T_x_func.m

@@ -0,0 +1,14 @@
+function y = T_x_func(t)
+    global T_r T_p f_n;
+    p = rem(t, T_r); % p = t - nT_r
+    n = round((t - p) / T_r);
+
+    if (p > T_p)
+        y = 0;
+    elseif (p <= 0)
+        y = 0;
+    else
+        y = exp(1j * 2 * pi * f_n(n + 1) * p);
+    end
+
+end

Radar simulation/configure_parameters.m

@@ -0,0 +1,40 @@
+global M N K d epsilon f_c Delta_f c scatter_coef B f_s T_p T_r r_0 velocity lambda delta_t max_t range_t len freqs slow_len
+
+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 = 100 / f_s; % 单载频脉冲下的采样周期 / 脉冲宽度
+T_r = T_p * 10;
+
+r_0 = 4; % 初始距离 r(0)
+velocity = 3e4; % 目标速度（假设目标做匀速直线运动）
+
+lambda = c / f_c; % 雷达工作波长
+
+% 仿真时间
+delta_t = 1e-2 * T_p;
+max_t = 100 * T_r;
+range_t = 0:delta_t:max_t-delta_t;
+len = round(max_t / delta_t);
+
+freqs = ((0:len-1) * f_s) / len;
+
+% 绘图
+figure_flag_1 = false;
+figure_flag_2 = false;
+figure_flag_3 = false;
+
+slow_len = 128;

Radar simulation/get_doppler.m

@@ -0,0 +1,51 @@
+function doppler = get_doppler(s_T, s_R)
+    global T_r delta_t len f_s c f_c velocity max_t slow_len;
+    
+    freqs = ((0:len-1) * f_s) / len;
+
+    slow_freqs = ((0:slow_len-1) * (1 / T_r)) / slow_len;
+    range_t_slow = 1:slow_len;
+
+    num_slow = max_t / T_r;
+    s_R_slow = zeros(1, slow_len);
+    k = floor(T_r / delta_t);
+
+    init_idx = 1;
+    while abs(s_R(init_idx)) == 0
+        init_idx = init_idx + 1;
+    end
+    init_idx
+
+    for i = 0: num_slow - 1
+        while abs(s_R(init_idx + i * k)) == 0
+            init_idx = init_idx + 1;
+        end
+        s_R_slow(i + 1) = s_R(init_idx + i * k);
+    end
+
+    s_R_fft = fft(s_R_slow);
+    [~, max_index_s_R] = max(s_R_fft);
+    freq_s_R = slow_freqs(max_index_s_R);
+
+    figure(2);
+    subplot(2,1,1);
+    plot(range_t_slow, abs(s_R_slow));
+    title(sprintf('s_R'));
+    subplot(2,1,2);
+    plot(slow_freqs, abs(s_R_fft));
+    title(sprintf('s_R_fft, freq = %E', freq_s_R));
+    xlabel('频率 (Hz)');
+
+    f_d = 1/T_r - freq_s_R;
+
+    fprintf("s_R_freq = %d\n", max_index_s_R);
+    fprintf("doppler_freq = %E\n\n", f_d);
+
+    fprintf("%E\n", 2 * f_c * velocity / c);
+
+    doppler_v = (c * f_d) / (2 * f_c);
+    
+    fprintf("doppler_v: %E\nvelocity: %E\n", doppler_v, velocity);
+
+    doppler = doppler_v;
+end

Radar simulation/get_downconversion_pulse.m

@@ -0,0 +1,41 @@
+function s_R_tlide = get_downconversion_pulse(s_R, f_n)
+
+    config_parameters;
+
+    s_R_tlide = complex(zeros(1, len));
+    real_s_R_tlide = zeros(1, len);
+    imag_s_R_tlide = zeros(1, len);
+
+    % down conversion
+    for t_idx = 1:len
+        t = range_t(t_idx);
+    
+        p = rem(t, T_r);
+    
+        % if (p > T_p)
+        %     continue
+        % end
+    
+        n = (t - p) / T_r;
+        s_R_tlide(t_idx) = s_R(t_idx) * exp(-1 * 1j * 2 * pi * f_n(t_idx) * (t - n * T_r));
+        real_s_R_tlide(t_idx) = real(s_R_tlide(t_idx));
+        imag_s_R_tlide(t_idx) = imag(s_R_tlide(t_idx));
+    end
+
+    if figure_flag_3
+        figure(3)
+        title("Received signal (after down conversion)")
+        subplot(2, 1, 1);
+        plot(range_t, imag_s_R_tlide);
+        xlabel("Time (s)");
+        ylabel("sin(\omega t)");
+        ylim([-1, 1]);
+    
+        subplot(2, 1, 2);
+        plot(range_t, real_s_R_tlide);
+        xlabel("Time (s)");
+        ylabel("cos(\omega t)");
+        ylim([-1, 1]);
+    end
+
+end

Radar simulation/get_range.m

@@ -0,0 +1,29 @@
+function [range, range_idx] = get_range(s_T, s_R)
+    global T_r delta_t len c r_0 range_t freqs;
+
+    range_N = T_r / delta_t;
+    s_T_first = [s_T(1:range_N), zeros(1, len - range_N)];
+    s_R_first = [s_R(1:range_N), zeros(1, len - range_N)];
+    
+    figure(1);
+    plot(range_t, s_T_first, color='red');
+    hold on;
+    plot(range_t, s_R_first, color='blue');
+
+    s_T_fft = fft(s_T_first, len);
+    s_R_fft = fft(s_R_first, len);
+
+    figure(2);
+    plot(freqs, s_T_fft, color='red');
+    hold on;
+    plot(freqs, s_R_fft, color='blue');
+    
+    t = conj(s_T_fft);
+    p = ifft(s_R_fft .* t);
+    norm_p = real(p).^2 + imag(p).^2;
+    % plot(range_t, norm_p);
+    
+    [~, range_idx] = max(norm_p);
+    range = range_t(range_idx) * c / 2;
+    fprintf("range = %f, r_0 = %f\n", range, r_0);
+end

Radar simulation/get_reflect_pulse.m

@@ -0,0 +1,49 @@
+function s_R = get_reflect_pulse(s_T)
+
+    config_parameters;
+
+    s_R = complex(zeros(1, len));
+    real_s_R = zeros(1, len);
+    imag_s_R = zeros(1, len);
+
+    for t_idx = 1:len
+        t = range_t(t_idx);
+
+        p = rem(t, T_r);
+
+        if (p > T_p)
+            continue
+        end
+
+        n = (t - p) / T_r;
+        t_diff = 2 * (r_0 + velocity * n * T_r) / c;
+
+        new_t_idx = floor((t + t_diff) / delta_t);
+
+        if new_t_idx > len
+            continue;
+        end
+
+        s_R(new_t_idx) = scatter_coef * s_T(t_idx);
+    end
+
+    real_s_R = real(s_R);
+    imag_s_R = imag(s_R);
+
+    if figure_flag_2
+        figure(2)
+        title("Received signal")
+        subplot(2, 1, 1);
+        plot(range_t, imag_s_R);
+        xlabel("Time (s)");
+        ylabel("sin(\omega t)");
+        ylim([-1, 1]);
+
+        subplot(2, 1, 2);
+        plot(range_t, real_s_R);
+        xlabel("Time (s)");
+        ylabel("cos(\omega t)");
+        ylim([-1, 1]);
+    end
+
+end

Radar simulation/get_transmitted_pulse.m

@@ -0,0 +1,40 @@
+function [C_n, f_n, s_T] = get_transmitted_pulse()
+
+    config_parameters;
+    
+    s_T = complex(zeros(1, len));
+    C_n = zeros(1, len);
+    f_n = zeros(1, len);
+    
+    for t_idx = 1:len
+        t = range_t(t_idx);
+        p = rem(t, T_r); % p = t - nT_r
+    
+        if (p > T_p)
+            continue
+        end
+    
+        C_n(t_idx) = floor(rand * (M - 1));
+        f_n(t_idx) = f_c + C_n(t_idx) * Delta_f;
+        s_T(t_idx) = exp(1j * 2 * pi * f_n(t_idx) * p);
+    end
+
+    real_s_T = real(s_T);
+    imag_s_T = imag(s_T);
+
+    if figure_flag_1
+        figure(1)
+        title("Transmitted signal")
+        subplot(2, 1, 1);
+        plot(range_t, imag_s_T);
+        xlabel("Time (s)");
+        ylabel("sin(\omega t)");
+        ylim([-1, 1]);
+    
+        subplot(2, 1, 2);
+        plot(range_t, real_s_T);
+        xlabel("Time (s)");
+        ylabel("cos(\omega t)");
+        ylim([-1, 1]);
+    end
+end

Radar simulation/main.m

@@ -0,0 +1,13 @@
+% This file to show the transmitted pulse for single signal
+
+clc;
+clear;
+
+config_parameters;
+
+[C_n, f_n, s_T] = get_transmitted_pulse();
+s_R = get_reflect_pulse(C_n, f_n);
+s_R_tlide = get_downconversion_pulse(s_R, f_n);
+
+r = get_range(s_T, s_R);
+dopp = get_doppler(s_T, s_R_tlide);

Radar simulation/main2.m

@@ -0,0 +1,37 @@
+%%
+clc;
+clear;
+
+configure_parameters;
+
+global C_n f_n;
+C_n = zeros(1, len);
+f_n = zeros(1, len);
+
+for t_idx = 1:len
+    C_n(t_idx) = floor(rand * (M - 1));
+    f_n(t_idx) = f_c + C_n(t_idx) * Delta_f;
+end
+
+T_x = zeros(1, len);
+R_x = zeros(1, len);
+R_d = zeros(1, len);
+
+for i = 1:len
+    t = range_t(i);
+    T_x(i) = T_x_func(t);
+    R_x(i) = R_x_func(t);
+    R_d(i) = R_d_func(t);
+end
+
+%%
+if true
+    figure(5)
+    plot(range_t, T_x, color="red");
+    hold on;
+    plot(range_t, R_d, color="blue");
+    xlim([0, 10 * T_r])
+end
+
+[r, range_idx] = get_range(T_x, R_x);
+doppler = get_doppler(T_x, R_d);

Radar simulation/r.m

@@ -0,0 +1,4 @@
+function x = r(t)
+    global r_0 velocity;
+    x = r_0 + velocity * t;
+end

Randomized Stepped FR Exploiting BS of ET/test.m

@@ -0,0 +1,17 @@
+N = 32e9;
+M = 3;
+eps = 1e-4;
+
+delta_1 = 24 * sqrt((M-1)/N) * log(M*N) * (2*sqrt(log(M * N) - log(eps)) + 1);
+delta_2 = 3/2 * sqrt((M-1)/N) * (2*sqrt(log(M * 2) - log(eps)) + 1);
+
+K = N * (1/8 - delta_1 - delta_2)^2 / (81 * M * log(M * N) * (1 + 2/3 * delta_2));
+
+x = 0:25;
+Ns = N * (1+randn(size(x)));
+delta_1 = 24 * sqrt((M-1)./N) * log(M.*N) * (2*sqrt(log(M .* N) - log(eps)) + 1);
+delta_2 = 3/2 * sqrt((M-1)./N) * (2*sqrt(log(M * 2) - log(eps)) + 1);
+Ks = Ns * (1/8 - delta_1 - delta_2)^2 / (81 * M * log(M .* Ns) * (1 + 2/3 * delta_2));
+
+rate = K .* M .* log(M .* Ns) ./ Ns;
+plot(x, rate)