Add Radar simulation

.gitignore

@@ -32,3 +32,9 @@ docs/site/
 Manifest.toml
 *.asv
 *.pdf
+
+0 - example
+*.mat
+*.fig
+*.tif
+*.bmp

1 - experiment/Chirp/Test.m

@@ -0,0 +1,20 @@
+N = 1e3 + 1;
+f_c = 1e6;
+t = linspace(0, pi / f_c, N);
+f_s = 1 / N;
+freqs = (0:N-1) / f_s;
+x = exp(-1j * 2 * pi * f_c * t);
+
+subplot(211);
+plot(t, imag(x));
+title("x(t) = sin(t)");
+xlabel("t");
+ylabel("x(t)");
+
+x_fft = abs(fft(imag(x)));
+subplot(212);
+plot(t, x_fft);
+title("X(f) = fft(x(t))");
+xlabel("f");
+ylabel("X(f)");
+

1 - experiment/Chirp/compress.m

@@ -0,0 +1,17 @@
+clc; clear;
+
+configure_parameters;
+global s_T s_R
+
+compressed_signal = conv(s_R, fliplr(s_T), 'same');
+
+figure(1);
+subplot(311);
+plot(t, s_T);
+title("Chirp signal");
+subplot(312);
+plot(t, s_R);
+title("Received Chirp signal");
+subplot(313);
+plot(t, compressed_signal);
+title("Matcher filter");

1 - experiment/Chirp/configure_parameters.m

@@ -0,0 +1,21 @@
+global f_s T_r simu_time t start_freq end_freq c len r_0 freqs s_T s_R
+
+f_s = 1e4;
+T_r = 1 / f_s;
+simu_time = 1;
+t = linspace(0, simu_time, f_s);
+start_freq = 0;
+end_freq = 1e7;
+c = 3e8;
+r_0 = 8e6;
+
+s_T = chirp(t, start_freq, simu_time, end_freq, 'linear') * 10;
+len = length(s_T);
+freqs = linspace(-f_s / 2, f_s / 2, len);
+
+s_R = zeros(size(s_T));
+delay = 2 * r_0 / c;
+delay_samples = round(delay * len) + 1;
+s_R(delay_samples:len) = 0.8 * s_T(1:len-delay_samples + 1);
+s_R = awgn(s_R, 1);
+

1 - experiment/Chirp/frequency_domain.m

@@ -0,0 +1,39 @@
+clc; clear;
+
+configure_parameters;
+global s_T s_R c
+
+% 频域上的匹配滤波
+s_T_fft = fft(s_T);
+s_R_fft = fft(s_R);
+s_R_fft_conj = conj(s_R_fft);
+
+matched_filter = ifft(s_T_fft .* s_R_fft_conj);
+[~, peak] = max(matched_filter);
+t_delay = 1 - peak / length(matched_filter);
+d = (c * t_delay) / 2;
+
+fprintf("r = %f, predict d = %f\n", r_0, d);
+fprintf("eps = %f\n", abs(r_0 - d));
+
+figure(1);
+subplot(311);
+plot(freqs, abs(s_T_fft));
+title("Chirp signal");
+subplot(312);
+plot(freqs, abs(s_R_fft));
+title("Received Chirp signal");
+subplot(313);
+plot(freqs, abs(s_T_fft .* s_R_fft_conj));
+title("Matcher filter");
+
+figure(2);
+subplot(311);
+plot(t, s_T);
+title("Chirp signal");
+subplot(312);
+plot(t, s_R);
+title("Received Chirp signal");
+subplot(313);
+plot(t, matched_filter);
+title("After matched filter");

1 - experiment/Chirp/time_domain.m

@@ -0,0 +1,24 @@
+clc; clear;
+
+configure_parameters;
+global s_T s_R c
+
+% 时域上的匹配滤波
+matched_filter = conv(s_R, fliplr(s_T), 'same');
+[~, peak] = max(matched_filter);
+t_delay = peak / length(matched_filter) - 0.5;
+d = (c * t_delay) / 2;
+
+fprintf("r = %f, predict d = %f\n", r_0, d);
+fprintf("eps = %f\n", abs(r_0 - d));
+
+figure(1);
+subplot(311);
+plot(t, s_T);
+title("Chirp signal");
+subplot(312);
+plot(t, s_R);
+title("Received Chirp signal");
+subplot(313);
+plot(t, matched_filter);
+title("Matcher filter");

1 - experiment/Circulant Matrix/circulant_matrix.m

@@ -0,0 +1,47 @@
+%%
+clc;
+clear;
+
+N = 100;
+eps = 1e-10;
+F_n = dftmtx(N);
+
+for i = 1: N
+    for j = 1: N
+        fprintf(string(F_n(i, j)) + ", ");
+    end
+    fprintf("\b\b\n");
+end
+
+fprintf("\n");
+
+% 检查傅里叶矩阵的正确性
+assert(abs(F_n(2, 2) ^ 2 - F_n(2, 3)) < eps);
+assert(abs(F_n(2, 2) ^ N - 1) < eps);
+
+%% 检查循环卷积
+x = randn(N, 1);
+z = randn(N, 1);
+p = cconv(x, z, N);
+
+% 检查循环卷积公式的正确性
+for i = 1: N
+    tmp = 0;
+    for j = 1: N
+        tmp = tmp + x(j) * z(mod(i-j+N, N)+1);
+    end
+    assert(abs(tmp - p(i)) < eps);
+end
+
+% 检查 x * z = A(x)z 的正确性
+A_x = toeplitz([x(1) fliplr(x(2:end)')], x);
+assert(is_vector_equal(p, A_x' * z, eps));
+
+% 检查卷积定理的正确性
+x_hat = fft(x);
+z_hat = fft(z);
+xz_hat = fft(p);
+assert(is_vector_equal(x_hat .* z_hat, xz_hat, eps));
+
+
+

1 - experiment/Circulant Matrix/is_vector_equal.m

@@ -0,0 +1,5 @@
+function flg = is_vector_equal(A, B, eps)   
+    diff = abs(A - B);
+    max_diff = max(diff);
+    flg = max_diff <= eps;
+end

Block RIP/Phase transitions/block_phase_transition.m

@@ -0,0 +1,39 @@
+function pt = block_phase_transition(N, M)
+    tau_min = 0;
+    tau_max = 10;
+    tau_interval = 0.05;
+    tau_range = tau_min:tau_interval:tau_max;
+    
+    K_min = 1;
+    K_max = 25;
+    K_range = K_min:K_max;
+    
+    pt = zeros(25, 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_block_integral(tau, 4);
+        end
+        save I.mat
+    end
+    
+    
+    for K_idx = 1:length(K_range)
+        K = K_range(K_idx);
+        f_set = zeros(length(tau_range), 1);
+    
+        for tau_idx = 1:length(tau_range)
+            tau = tau_range(tau_idx);
+            f_set(tau_idx) = 1/2 * (K * (2 * M + tau ^ 2) + (N - K) * I(tau_idx));
+        end
+    
+        pt(K_idx) = min(f_set);
+    end
+end

Block RIP/Phase transitions/block_phase_transition_0.m

@@ -0,0 +1,40 @@
+function pt = block_phase_transition_0()
+    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 = 25;
+    s_b_range = s_b_min:s_b_max;
+    
+    pt = zeros(25, 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_block_integral(tau, 4);
+        end
+    
+    end
+    
+    save I.mat
+    
+    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
+    
+        pt(s_b_idx) = min(f_set);
+    end
+end

Block RIP/Phase transitions/block_phase_transition_2.m

@@ -0,0 +1,20 @@
+function pt = block_phase_transition_2(N, M)
+    K_min = 1;
+    K_max = 25;
+    K_range = K_min:K_max;
+    pt = zeros(25, 1);
+    d = N;
+    m = M;
+    
+    for K_idx = 1:length(K_range)
+        K = K_range(K_idx);
+        s = K;
+        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);
+        pt(K_idx) = n;
+    end
+end

Block RIP/Phase transitions/calc_block_integral.m

@@ -0,0 +1,15 @@
+function I = calc_block_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) * x ^ (m - 1) / (2 ^ (m/2 - 1) * gamma(m/2));
+        I = double(int(f, [tau, +inf]));
+    end
+
+end

Block RIP/Phase transitions/draw.m

@@ -0,0 +1,52 @@
+clc;
+clear;
+
+% epi = 0;
+% epi = 0.1;
+% epi = 0.5;
+epi = 1;
+CONTOURF = false;
+
+filename = strcat("FAR_block_sparsity_phase_epi_", string(epi), ".mat");
+load(filename);
+
+figure;
+x_range = 1:size(result, 2);
+y_range = 1:size(result, 1);
+
+if CONTOURF
+    handler = contourf(x_range, y_range, result);
+else
+    handler = pcolor(x_range, y_range, result);
+    shading interp;
+    colorbar;
+    xlim([1, size(result, 2)]);
+    ylim([1, size(result, 1)]);
+    box on;
+end
+
+ax = gca;
+grid(ax, 'on');
+set(ax, 'Visible', 'on');
+
+pt = block_phase_transition_2(128, 4);
+l = line(1:size(pt), pt);
+l.Color = "w";
+l.LineWidth = 5;
+
+lgd = legend("", "Theoretical curve");
+set(lgd, 'Location', 'southeast');
+set(lgd, 'TextColor', 'white');
+set(lgd, 'Box', 'off');
+set(lgd, 'Color', 'none');
+
+xlabel("Number of targets ($K$)", 'Interpreter', 'latex');
+ylabel("Number of measurement ($n$)", 'Interpreter', 'latex');
+
+% ax.XDisplayLabels = nan(size(ax.XDisplayData));
+% ax.YDisplayLabels = nan(size(ax.YDisplayData));
+% ax.ColorbarVisible = 'on';
+% ax.GridVisible = 'off';
+% ax.Colormap = parula(100);
+
+grid on;

Block RIP/Phase transitions/theoretic.m

@@ -0,0 +1,8 @@
+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);
+end

Radar simulation/FAR multi target/config_parameters.m

@@ -0,0 +1,38 @@
+global M N K d epsilon f_c Delta_f c scatter_coef B f_s T_p T_r lambda delta_t max_t range_t len freqs slow_len num_pulse
+
+N = 128; % 脉冲个数
+M = 5; % 频点个数
+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;
+
+T_p = 1e-6 / 3; % 单载频脉冲下的采样周期 / 脉冲宽度
+f_s = 100 / T_p; % 快时间采样率
+T_r = T_p * 6;
+
+lambda = c / f_c; % 雷达工作波长
+
+% 仿真时间
+num_pulse = 180;
+delta_t = 1e-2 * T_p;
+max_t = num_pulse * 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 = 1e4;

Radar simulation/FAR multi target/get_range.m

@@ -0,0 +1,22 @@
+function [range, range_idx] = get_range(s_T, s_R, r_0)
+    global T_r delta_t len c range_t;
+
+    % 只取第一个 T_r 的数据计算
+    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)];
+
+    s_T_fft = fft(s_T_first, len);
+    s_R_fft = fft(s_R_first, len);
+
+    % 发射信号和回波信号做相关
+    p = ifft(s_R_fft .* conj(s_T_fft));
+    norm_p = real(p).^2 + imag(p).^2;
+    [~, range_idx] = max(norm_p);
+
+    range = range_t(range_idx) * c / 2;
+    % fprintf("predict range = %f\n", range);
+    % fprintf("origin range = %f\n", r_0);
+    % fprintf("err = %f\n", range - r_0);
+    % fprintf("\n");
+end

Radar simulation/FAR multi target/get_velocity.m

@@ -0,0 +1,47 @@
+function velocity = get_velocity(s_T, s_R, origin_velocity)
+    global T_r delta_t c f_c max_t slow_len;
+    freq_s_T = 1 / T_r;
+
+    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);
+
+    idx = 1;
+
+    while abs(s_R(idx)) == 0
+        idx = idx + 1;
+    end
+
+    for i = 0:num_slow - 1
+
+        while abs(s_R(idx + i * k)) == 0
+            idx = idx + 1;
+        end
+        s_R_slow(i + 1) = s_R(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 = freq_s_T - freq_s_R;
+    velocity = (c * f_d) / (2 * f_c);
+
+    fprintf("predict velocity: %f\n", velocity);
+    fprintf("origin velocity: %f\n", origin_velocity);
+    fprintf("err = %f\n", velocity - origin_velocity);
+    fprintf("\n");
+
+end

Radar simulation/FAR multi target/main.m

@@ -0,0 +1,86 @@
+%%
+clc;
+clear;
+
+config_parameters;
+
+range = [250, 140, 167, 128, 12];
+velocity = [16, 76, 155, 125, 463];
+
+global C_n f_n;
+C_n = zeros(M, len);
+f_n = zeros(M, len);
+
+for n_idx = 1:M
+
+    for t_idx = 1:len
+        C_n(n_idx, t_idx) = floor(rand * (M - 1));
+        f_n(n_idx, t_idx) = f_c + C_n(n_idx, t_idx) * Delta_f;
+    end
+
+end
+
+T_x = zeros(M, len);
+R_x = zeros(M, len);
+R_d = zeros(M, len);
+
+for n_idx = 1:M
+
+    for t_idx = 1:len
+        t = range_t(t_idx);
+        T_x(n_idx, t_idx) = T_x_func(t);
+        R_x(n_idx, t_idx) = R_x_func(t, range(n_idx), velocity(n_idx));
+        R_d(n_idx, t_idx) = R_d_func(t, range(n_idx), velocity(n_idx));
+    end
+
+end
+
+save temp.mat
+
+%%
+
+clc
+clear
+load temp.mat
+
+pred_range = zeros(1, M);
+range_idx = zeros(1, M);
+pred_velocity = zeros(1, M);
+
+for n_idx = 1:M
+
+    [pred_range(n_idx), range_idx(n_idx)] = get_range(T_x(n_idx, :), R_x(n_idx, :), range(n_idx));
+    pred_velocity(n_idx) = get_velocity(T_x(n_idx, :), R_x(n_idx, :), velocity(n_idx));
+
+end
+
+%%
+r_range = 1:500;
+v_range = 1:500;
+
+[X, Y] = meshgrid(r_range, v_range);
+pred_z = zeros(length(r_range), length(v_range));
+origin_z = zeros(length(r_range), length(v_range));
+
+for n_idx = 1:M
+    pred_z(round(pred_range(n_idx)), round(pred_velocity(n_idx))) = 1;
+    origin_z(round(range(n_idx)), round(velocity(n_idx))) = 1;
+end
+
+figure(1)
+subplot(211);
+mesh(X, Y, pred_z);
+title("Predict range-velocity reconstruction under M = 5");
+xlabel("Range (m)");
+ylabel("Velocity (m/s)");
+
+subplot(212);
+mesh(X, Y, origin_z);
+title("Origin range-velocity data under M = 5");
+xlabel("Range (m)");
+ylabel("Velocity (m/s)");
+
+
+
+
+

Radar simulation/FAR single target/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/FAR single target/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/FAR single target/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/FAR single target/config_parameters.m

@@ -19,6 +19,8 @@ T_p = 1e-6 / 3; % 单载频脉冲下的采样周期 / 脉冲宽度
 f_s = 100 / T_p; % 快时间采样率
 T_r = T_p * 10;
 
+r_0 = 30; % 初始距离 r(0)
+velocity = 3e2; % 目标速度（假设目标做匀速直线运动）
 lambda = c / f_c; % 雷达工作波长
 
 % 仿真时间

Radar simulation/FAR single target/get_doppler.m

@@ -9,19 +9,19 @@ function v = get_doppler(s_T, s_R)
     s_R_slow = zeros(1, slow_len);
     k = floor(T_r / delta_t);
 
-    idx = 1;
+    init_idx = 1;
 
-    while abs(s_R(idx)) == 0
-        idx = idx + 1;
+    while abs(s_R(init_idx)) == 0
+        init_idx = init_idx + 1;
     end
 
     for i = 0:num_slow - 1
 
-        while abs(s_R(idx + i * k)) == 0
-            idx = idx + 1;
+        while abs(s_R(init_idx + i * k)) == 0
+            init_idx = init_idx + 1;
         end
 
-        s_R_slow(i + 1) = s_R(idx + i * k);
+        s_R_slow(i + 1) = s_R(init_idx + i * k);
     end
 
     s_R_fft = fft(s_R_slow);

Radar simulation/FAR single target/main.m

@@ -0,0 +1,37 @@
+%%
+clc;
+clear;
+
+config_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 false
+    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/FAR single target/r.m

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

Radar simulation/csdn/Main.m

@@ -0,0 +1,14 @@
+global c f_c B T_r N F_s M k lambda
+
+c = 3e8; %光速
+f_c = 70e9;  %发射信号载频 中心频率
+B = 500e6;  %发射信号带宽
+T_r = 1e-5;  %扫频时间 也就是周期
+N = 256;  %采样点
+F_s = round(N / T_r); %采样率
+M = 256;  %chirp的数目
+k = B/T_r; %chirp斜率
+lambda = c / (f_c - B/2);
+
+ksy_main;
+draw;

Radar simulation/csdn/calc.m

@@ -0,0 +1,38 @@
+% https://blog.csdn.net/Xiao_Jie123/article/details/115296169
+
+%% 超参数
+c = 3e8; %光速
+f_c = 76.5e9;  %发射信号载频 中心频率
+B = 500e6;  %发射信号带宽
+T_r = 10e-6;  %扫频时间 也就是周期
+N = 256;  %采样点
+F_s = 25.6e6; %采样率
+M = 256;  %chirp的数目
+k = B/T_r; %chirp斜率
+index = 1:1:N;  %产生点向量
+IF_mat = zeros(M,N);  %存储带有噪声的中频信号
+%% 发射信号参数
+AT = 10;  %发射信号增益
+t = 0:1/F_s:T_r-1/F_s;  %时间向量 确定256个点在一个Tr中的每个时刻
+t = t - T_r/2; %将fc作为中心频率
+
+%% 回波信号参数
+distance = 50;  %目标距离雷达50m的距离
+t_d = 2 * distance / c;  %目标距离雷达的延迟
+velocity = -20;  %目标距雷达的相对速度为30m/s
+f_d = 2 * (f_c - B/2) * velocity / c;  %多普勒频移
+AR = 0.8;  %回波信号衰减的比例值
+
+%% 生成数据
+for i = 1:1:M  %chirp的循环
+    s_T = AT*exp((1i*2*pi)*(f_c*(t+i*T_r)+k/2*t.^2)); %发射信号
+    s_R = AR*AT*exp((1i*2*pi)*((f_c-f_d)*(t-t_d+i*T_r)+k/2*(t-t_d).^2));  %回波信号
+    %% 求回波信号的共轭
+    s_R_conj = conj(s_R);  %求回波信号的共轭
+    %% 求中频信号
+    IF = s_T .* s_R_conj;  %求中频信号
+    SNR = 10;  %信噪比
+    IF_with_Noise = awgn(IF,SNR,'measured');  %给中频信号加高斯白噪声，在添加噪声的时候，要进行能量的测量
+    IF_mat(i,:) = IF_with_Noise;  %将带有噪声的中频信号保存
+end
+save('Ego_vehicle.mat', 'IF_mat');  %进行数据的保存

Radar simulation/csdn/draw.m

@@ -0,0 +1,44 @@
+% https://blog.csdn.net/Xiao_Jie123/article/details/115296169
+
+clc; clear;
+
+global c T_r F_s k lambda
+
+%% 加载数据
+% IF_mat = cell2mat(struct2cell(load('Ego_vehicle.mat','IF_mat')));
+IF_mat = cell2mat(struct2cell(load('Ego_vehicle_ksy.mat','IF_mat')));
+[N, M] = size(IF_mat);
+%% 生成窗
+range_win = hamming(N);  %生成range窗
+doppler_win = hamming(M);  %生成doppler窗
+%% range fft
+for i = 1:1:N
+    temp = IF_mat(i,:) .* range_win';
+    % temp_fft = fftshift(fft(temp,N));
+    temp_fft = fft(temp,N);
+    IF_mat(i,:) = temp_fft;
+end
+%% doppler fft
+for j = 1:1:M
+    temp = IF_mat(:,j) .* doppler_win;
+    temp_fft = fftshift(fft(temp,M));
+    IF_mat(:,j) = temp_fft;
+end
+%% 画图
+figure;
+distance_temp = (-N/2:N/2 - 1) * F_s * c / N / 2 / k;
+% distance_temp = (0:N - 1) * F_s * c / N / 2 / k;
+speed_temp = (-M / 2:M / 2 - 1) * lambda / T_r / M / 2;
+[X,Y] = meshgrid(distance_temp,speed_temp);
+mesh(X,Y,(abs(IF_mat)));
+xlabel('距离(m)');
+ylabel('速度(m/s)');
+zlabel('信号幅值');
+title('2维FFT处理三维视图');
+
+figure;
+speed_temp = -speed_temp;
+imagesc(distance_temp,speed_temp,abs(IF_mat));
+title('距离-多普勒视图');
+xlabel('距离(m)');
+ylabel('速度(m/s)');

Radar simulation/csdn/ksy_main.m

@@ -0,0 +1,39 @@
+%% 超参数
+global c f_c B T_r N F_s M k lambda
+index = 1:1:N;  %产生点向量
+IF_mat = zeros(M,N);  %存储带有噪声的中频信号
+t = 0:1/F_s:T_r-1/F_s;  %时间向量 确定256个点在一个Tr中的每个时刻
+t = t - T_r/2; %将fc作为中心频率
+
+dist = 10;  %目标距离雷达50m的距离
+t_d = 2 * dist / c;  %目标距离雷达的延迟
+velocity = -20;  %目标距雷达的相对速度为30m/s
+f_d = 2 * f_c * velocity / c;  %多普勒频移
+% f_d = 2 * (f_c - B/2) * velocity / c;  %多普勒频移
+scat_coef = 0.8;  %回波信号衰减的比例值
+
+s_T = chirp(t, f_c, t(end), f_c + B);
+pad = round(t_d * F_s);
+s_R = [zeros(1, pad), scat_coef * s_T(1: end - pad)];
+
+for i = 1: M
+    s_T = 1*exp((1i*2*pi)*(f_c*(t+i*T_r)+k/2*t.^2)); %发射信号
+    s_R = scat_coef*exp((1i*2*pi)*((f_c-f_d)*(t-t_d+i*T_r)+k/2*(t-t_d).^2));  %回波信号
+
+    IF = s_T .* conj(s_R);
+    % SNR = 10;
+    % IF_Noise = awgn(IF, SNR, 'measured');
+    IF_Noise = IF;
+    IF_mat(i, :) = IF_Noise;
+end
+
+save('Ego_vehicle_ksy.mat','IF_mat');
+
+% figure;
+% subplot(211);
+% plot(t, s_T);
+% ylim([-1.3, 1.3]);
+% 
+% subplot(212);
+% plot(t, s_R);
+% ylim([-1.3, 1.3]);

Radar simulation/multi target/main.m

@@ -1,50 +0,0 @@
-%%
-clc;
-clear;
-
-config_parameters;
-
-target_num = 2;
-range = [250, 140];
-velocity = [3e2, 15];
-
-global C_n f_n;
-C_n = zeros(target_num, len);
-f_n = zeros(target_num, len);
-
-for n_idx = 1:target_num
-
-    for t_idx = 1:len
-        C_n(n_idx, t_idx) = floor(rand * (M - 1));
-        f_n(n_idx, t_idx) = f_c + C_n(n_idx, t_idx) * Delta_f;
-    end
-
-end
-
-T_x = zeros(target_num, len);
-R_x = zeros(target_num, len);
-R_d = zeros(target_num, len);
-
-for n_idx = 1:target_num
-
-    for t_idx = 1:len
-        t = range_t(t_idx);
-        T_x(n_idx, t_idx) = T_x_func(t);
-        R_x(n_idx, t_idx) = R_x_func(t, range(n_idx), velocity(n_idx));
-        R_d(n_idx, t_idx) = R_d_func(t, range(n_idx), velocity(n_idx));
-    end
-
-end
-
-%%
-
-r = zeros(1, target_num);
-range_idx = zeros(1, target_num);
-doppler = zeros(1, target_num);
-
-for n_idx = 1:target_num
-
-    [r(n_idx), range_idx(n_idx)] = get_range(T_x(n_idx, :), R_x(n_idx, :));
-    doppler(n_idx) = get_doppler(T_x(n_idx, :), R_d(n_idx, :));
-
-end