clc;
clear;
close all;

%% 参数设置
sigma_n = 0.1;
gamma = 0.5;

% 信号参数
B = 5e5;        %信号带宽
Tp = 100e-6;    %脉宽100us
fs = 2 * B;     %采样频率
Ts = 1 / fs;    %采样周期
K = B / Tp;     %线性调频率
fc = 1e8;       %载波频率

Tr = 1e-3;

t = 0: 1/fs: Tr - 1/fs;
t2 = 0: 1/fs/2: Tr - 1/fs/2;

c = 3e8; % 光速

distance_max = (Tr-Tp) * c / 2;

target_scattering = [0.8, 1, 0.9]; %扩展目标各点散射强度


%% 生成矩阵 A
% 生成发射信号 signal_t 及 
N = Tr * fs;
N_high = Tp * fs;

signal_t  = zeros(1, N);
signal_td  = zeros(1, N);
for i = 1:N_high
    tp = (i - 1) * (1 / fs);
    signal_t(1, i) = exp(1j*2*pi*(fc*tp+0.5*K*tp.^2));
    tp2 = (i - 0.5) * (1 / fs);
    signal_td(1, i + 1) = exp(1j*2*pi*(fc*tp2+0.5*K*tp2.^2));
end

A = generate_matrix_new(signal_t.', signal_td.');

J = A'*A;
lambda_J=eig(J);

figure(1)
spy(A)


%% 生成回波 y
% 设置目标 - 扩展目标，由三个点组成
distance1 = 52000;
tau1 = distance1 * 2 / c;
n_tau1 = round(tau1 * fs);

alpha1 = 0.6; % 扩展目标整体散射强度
signal_r1_t1  = zeros(1, N);
signal_r2_t1  = zeros(1, N);
signal_r3_t1  = zeros(1, N);

for i = 1:N
    temp = i - n_tau1;
    if temp >= 1 && temp <= N_high
        signal_r1_t1(1, i) = alpha1 * target_scattering(1) * signal_t(1, temp);
        if i + 1 <= N
            signal_r2_t1(1, i + 1) = alpha1 * target_scattering(2) * signal_t(1, temp);
        end
        if i + 1 <= N
            signal_r3_t1(1, i + 2) = alpha1 * target_scattering(3) * signal_t(1, temp);
        end
    end
end
signal_r1 = signal_r1_t1 + signal_r2_t1 + signal_r3_t1;

% 回波
signal_r = signal_r1;

% 加入噪声
noise = random('Normal', 0, sigma_n/sqrt(2), 1, length(signal_r)) + 1j * random('Normal', 0, sigma_n/sqrt(2), 1, length(signal_r));
signal_r_n = signal_r + noise;

y = signal_r_n.';
figure(2)
subplot(211);
plot(t, real(signal_t));
title('发射信号')
xlabel('时间')
ylabel('幅度')

subplot(212);
plot(t, real(y))
title('接收信号(y)')
xlabel('时间')
ylabel('幅度')

%% 理论 x
distance_node = round((distance1 * 2 / c) * fs);
x_t = zeros(1, 2 * N);
for i = 1: length(target_scattering)
    x_t((distance_node + i) * 2 - 1) = alpha1 * target_scattering(i);
end
x = x_t.';

figure(3)
plot(t2, x);
title('目标散射点(x)')
xlabel('时间')
ylabel('幅度')



%% 验证
y_t = A * x;
y_r = signal_r.';

figure(4)
subplot(311);
plot(real(y_r))
title('实际回波')

subplot(312);
plot(real(y_t));
title('计算结果')

subplot(313);
plot(abs(y_r - y_t));
title('差异')