From ec0f95c6394b1d8a8823653cdbf6011d3e6e85eb Mon Sep 17 00:00:00 2001
From: Ksyer <a593510703@gmail.com>
Date: Sat, 27 Jan 2024 20:58:23 +0800
Subject: [PATCH] Add Radar simulation

---
 Radar simulation/R_d_func.m                 | 13 ++++++
 Radar simulation/R_x_func.m                 |  5 ++
 Radar simulation/T_x_func.m                 | 14 ++++++
 Radar simulation/configure_parameters.m     | 40 ++++++++++++++++
 Radar simulation/get_doppler.m              | 51 +++++++++++++++++++++
 Radar simulation/get_downconversion_pulse.m | 41 +++++++++++++++++
 Radar simulation/get_range.m                | 29 ++++++++++++
 Radar simulation/get_reflect_pulse.m        | 49 ++++++++++++++++++++
 Radar simulation/get_transmitted_pulse.m    | 40 ++++++++++++++++
 Radar simulation/main.m                     | 13 ++++++
 Radar simulation/main2.m                    | 37 +++++++++++++++
 Radar simulation/r.m                        |  4 ++
 12 files changed, 336 insertions(+)
 create mode 100644 Radar simulation/R_d_func.m
 create mode 100644 Radar simulation/R_x_func.m
 create mode 100644 Radar simulation/T_x_func.m
 create mode 100644 Radar simulation/configure_parameters.m
 create mode 100644 Radar simulation/get_doppler.m
 create mode 100644 Radar simulation/get_downconversion_pulse.m
 create mode 100644 Radar simulation/get_range.m
 create mode 100644 Radar simulation/get_reflect_pulse.m
 create mode 100644 Radar simulation/get_transmitted_pulse.m
 create mode 100644 Radar simulation/main.m
 create mode 100644 Radar simulation/main2.m
 create mode 100644 Radar simulation/r.m

diff --git a/Radar simulation/R_d_func.m b/Radar simulation/R_d_func.m
new file mode 100644
index 0000000..8e9d948
--- /dev/null
+++ b/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
diff --git a/Radar simulation/R_x_func.m b/Radar simulation/R_x_func.m
new file mode 100644
index 0000000..81f333e
--- /dev/null
+++ b/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
diff --git a/Radar simulation/T_x_func.m b/Radar simulation/T_x_func.m
new file mode 100644
index 0000000..eec5482
--- /dev/null
+++ b/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
diff --git a/Radar simulation/configure_parameters.m b/Radar simulation/configure_parameters.m
new file mode 100644
index 0000000..46daf15
--- /dev/null
+++ b/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;
diff --git a/Radar simulation/get_doppler.m b/Radar simulation/get_doppler.m
new file mode 100644
index 0000000..e756d0c
--- /dev/null
+++ b/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
diff --git a/Radar simulation/get_downconversion_pulse.m b/Radar simulation/get_downconversion_pulse.m
new file mode 100644
index 0000000..4f79479
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/Radar simulation/get_range.m b/Radar simulation/get_range.m
new file mode 100644
index 0000000..b2e36ed
--- /dev/null
+++ b/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
diff --git a/Radar simulation/get_reflect_pulse.m b/Radar simulation/get_reflect_pulse.m
new file mode 100644
index 0000000..cc91243
--- /dev/null
+++ b/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
diff --git a/Radar simulation/get_transmitted_pulse.m b/Radar simulation/get_transmitted_pulse.m
new file mode 100644
index 0000000..a7e75db
--- /dev/null
+++ b/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
diff --git a/Radar simulation/main.m b/Radar simulation/main.m
new file mode 100644
index 0000000..5a9edf4
--- /dev/null
+++ b/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);
\ No newline at end of file
diff --git a/Radar simulation/main2.m b/Radar simulation/main2.m
new file mode 100644
index 0000000..b3d3d19
--- /dev/null
+++ b/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);
diff --git a/Radar simulation/r.m b/Radar simulation/r.m
new file mode 100644
index 0000000..24659bb
--- /dev/null
+++ b/Radar simulation/r.m	
@@ -0,0 +1,4 @@
+function x = r(t)
+    global r_0 velocity;
+    x = r_0 + velocity * t;
+end