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

Macroscopic Analysis of VAMP in a MMS/cVAMP.m

@@ -1,77 +1,77 @@
-% Input:y,A,lambda,tau,Kit
+% Input:y,A,lambda,tau,Kit
-% Output:x_hat_wl,x_hat_d
+% Output:x_hat_wl,x_hat_d
-
+
-function [x_hat_wl, x_hat_d] = cVAMP(y, A, lambda, tau, Kit)
+function [x_hat_wl, x_hat_d] = cVAMP(y, A, lambda, tau, Kit)
-    % Initialization
+    % Initialization
-    gamma = 768 ./ 1024;
+    gamma = 768 ./ 1024;
-    k = 0;
+    k = 0;
-    p = ctranspose(A) * y;
+    p = ctranspose(A) * y;
-    h_1 = p;
+    h_1 = p;
-    Q_1 = gamma;
+    Q_1 = gamma;
-    tau_d = 1;
+    tau_d = 1;
-
+
-    % while
+    % while
-    while (k < Kit) && (tau_d > tau)
+    while (k < Kit) && (tau_d > tau)
-        % Factorized Part
+        % Factorized Part
-        x_1 = ST(h_1, lambda, Q_1); % \hat{x}_1^{(k)}
+        x_1 = ST(h_1, lambda, Q_1); % \hat{x}_1^{(k)}
-        chi_1 = F1(x_1, lambda, Q_1); % \chi_1^{(k)}
+        chi_1 = F1(x_1, lambda, Q_1); % \chi_1^{(k)}
-        % Message Passing
+        % Message Passing
-        h_2 = x_1 ./ chi_1 - h_1; % h_2^{(k)}
+        h_2 = x_1 ./ chi_1 - h_1; % h_2^{(k)}
-        Q_2 = 1 ./ chi_1 - Q_1; % \hat{Q}_2^{(k)}
+        Q_2 = 1 ./ chi_1 - Q_1; % \hat{Q}_2^{(k)}
-        % Gaussian Part
+        % Gaussian Part
-        t1 = (p + h_2) ./ Q_2;
+        t1 = (p + h_2) ./ Q_2;
-        t2 = ctranspose(A) * (A * (p + h_2)) / ((Q_2 + 1) * Q_2);
+        t2 = ctranspose(A) * (A * (p + h_2)) / ((Q_2 + 1) * Q_2);
-        x_2 = t1 + t2; % \hat{x}_2^{(k)}
+        x_2 = t1 + t2; % \hat{x}_2^{(k)}
-        chi_2 = gamma ./ (Q_2 + 1) + (1 - gamma) ./ Q_2;
+        chi_2 = gamma ./ (Q_2 + 1) + (1 - gamma) ./ Q_2;
-        % Message Passing
+        % Message Passing
-        h_1_next = x_2 ./ chi_2 - h_2;
+        h_1_next = x_2 ./ chi_2 - h_2;
-        Q_1_next = 1 ./ chi_2 - Q_2;
+        Q_1_next = 1 ./ chi_2 - Q_2;
-        tau_d = norm(h_1_next - h_1) ./ norm(h_1_next);
+        tau_d = norm(h_1_next - h_1) ./ norm(h_1_next);
-        k = k + 1;
+        k = k + 1;
-        % output
+        % output
-        x_hat_wl = x_1;
+        x_hat_wl = x_1;
-        x_hat_d = h_1_next ./ Q_1_next;
+        x_hat_d = h_1_next ./ Q_1_next;
-        % next
+        % next
-        h_1 = h_1_next;
+        h_1 = h_1_next;
-        Q_1 = Q_1_next;
+        Q_1 = Q_1_next;
-    end
+    end
-
+
-end
+end
-
+
-% SoftThreshold function
+% SoftThreshold function
-function x = ST(h_1, lambda, Q_1)
+function x = ST(h_1, lambda, Q_1)
-    [N, M] = size(h_1);
+    [N, M] = size(h_1);
-    x = zeros(N, M);
+    x = zeros(N, M);
-
+
-    for i = 1:N
+    for i = 1:N
-        % sign = h_1(i) ./ abs(h_1(i));
+        % sign = h_1(i) ./ abs(h_1(i));
-        diff = abs(h_1(i)) - lambda(i);
+        diff = abs(h_1(i)) - lambda(i);
-        x(i) = sign(h_1(i)) .* (diff ./ Q_1) .* SF(diff);
+        x(i) = sign(h_1(i)) .* (diff ./ Q_1) .* SF(diff);
-    end
+    end
-
+
-end
+end
-
+
-% Heaviside's step function
+% Heaviside's step function
-function v = SF(a)
+function v = SF(a)
-    % if a > 0
+    % if a > 0
-    %     v = 1;
+    %     v = 1;
-    % elseif a == 0
+    % elseif a == 0
-    %     v = 0;  % at zero points
+    %     v = 0;  % at zero points
-    % else
+    % else
-    %     v = 0;
+    %     v = 0;
-    % end
+    % end
-    v = heaviside(a);
+    v = heaviside(a);
-end
+end
-
+
-% SoftThreshold function
+% SoftThreshold function
-function v = F1(x_1, lambda, Q_1)
+function v = F1(x_1, lambda, Q_1)
-    [N, M] = size(x_1);
+    [N, M] = size(x_1);
-    count = 0;
+    count = 0;
-
+
-    for i = 1:N
+    for i = 1:N
-        temp = Q_1 .* abs(x_1(i)) + lambda(i);
+        temp = Q_1 .* abs(x_1(i)) + lambda(i);
-        count = count + (2 - lambda(i) ./ temp) .* SF(abs(x_1(i)));
+        count = count + (2 - lambda(i) ./ temp) .* SF(abs(x_1(i)));
-    end
+    end
-
+
-    v = count ./ (2 .* N .* Q_1);
+    v = count ./ (2 .* N .* Q_1);
-end
+end

Macroscopic Analysis of VAMP in a MMS/demo1.m

@@ -1,46 +1,46 @@
-% Input:y,A,lambda,tau,Kit
+% Input:y,A,lambda,tau,Kit
-% Output:x_hat_wl,x_hat_d
+% Output:x_hat_wl,x_hat_d
-clear;
+clear;
-clc;
+clc;
-% rng(1); % 随机种子
+% rng(1); % 随机种子
-
+
-%% test_稀疏向量
+%% test_稀疏向量
-% 设定稀疏度
+% 设定稀疏度
-k = 100; % 设定稀疏度
+k = 100; % 设定稀疏度
-
+
-% 构造感知矩阵D
+% 构造感知矩阵D
-m = 768; % 感知矩阵行数
+m = 768; % 感知矩阵行数
-n = 1024; % 感知矩阵列数 （n>>m）
+n = 1024; % 感知矩阵列数 （n>>m）
-
+
-% D = randn(m,n);     % 生成满足高斯分布的感知矩阵 64*256
+% D = randn(m,n);     % 生成满足高斯分布的感知矩阵 64*256
-
+
-F = dftmtx(n);
+F = dftmtx(n);
-row_indices = randperm(n, m);
+row_indices = randperm(n, m);
-D = F(row_indices, :);
+D = F(row_indices, :);
-
+
-% 构造稀疏信号X——共n个元素，其中k个元素不为0
+% 构造稀疏信号X——共n个元素，其中k个元素不为0
-X = zeros(n, 1);
+X = zeros(n, 1);
-index = randperm(n, k);
+index = randperm(n, k);
-val = randn(1, k);
+val = randn(1, k);
-X(index) = val;
+X(index) = val;
-
+
-% 得到观测矩阵（压缩后）
+% 得到观测矩阵（压缩后）
-A = D * X;
+A = D * X;
-
+
-%%
+%%
-% % 通过cVMAP算法完成恢复X，得到恢复后信号
+% % 通过cVMAP算法完成恢复X，得到恢复后信号
-lambda = ones(n, 1) ./ 10;
+lambda = ones(n, 1) ./ 10;
-[x_hat_wl, x_hat_d] = cVAMP(A, D, lambda, 1e-4, 200);
+[x_hat_wl, x_hat_d] = cVAMP(A, D, lambda, 1e-4, 200);
-
+
-%%
+%%
-% 显示结果
+% 显示结果
-figure;
+figure;
-subplot(3, 1, 1)
+subplot(3, 1, 1)
-stem(X);
+stem(X);
-title('origin signal')
+title('origin signal')
-subplot(3, 1, 2)
+subplot(3, 1, 2)
-stem(x_hat_wl);
+stem(x_hat_wl);
-title('restored signal')
+title('restored signal')
-subplot(3, 1, 3)
+subplot(3, 1, 3)
-stem(X - x_hat_wl);
+stem(X - x_hat_wl);
-title('differ')
+title('differ')

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
+    variable x(M, N) complex
-    minimize(norm(x, 1))
+        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)