雷达成像 Matlab 仿真 3 —— 多目标检测

脉冲压缩能够将信号压缩变窄，使得在时域上有重叠的两个信号能够分离，使用脉冲压缩和加窗抑制旁瓣，能够实现多目标的检测。

1. 多目标信号的回波

使用 matlab 时，能用矩阵运算就尽量不用 for 循环，因为矩阵运算能够加速运算。
当雷达探测范围内在某方向上有多个目标时，看可以用矩阵的每一行存储每一个目标的回波信号，所有行的回波信号叠加得到实际接收到的回波。
接下来写一个检测多目标脉冲压缩函数

1.1 参数设置

首先设置脉冲宽度、带宽、探测的最大距离和最小距离、点目标的距离（向量）和点目标的 RCS。接着计算一些需要用到的变量值：Twid, Fs,Ts 和 Nwid

function LFM_radar(T, B, Rmin, Rmax, R, RCS)
% inputs:
%   T: pulse duration
%   B: chirp frequency modulation bandwidth
%   Rmin, Rmax: range bin
%   R: position of ideal point targets
%   RCS: radar cross section of ideal targets

clear; clc;
set(0,'defaultfigurecolor', 'w');

%% default parameters
if nargin == 0
    T = 10e-6;
    B = 30e6;
    Rmin = 10000;
    Rmax = 15000;
    R = [11000, 12000, 12005, 13000, 13020];
    RCS = [1 1 1 1 1];
end

%% parameters
C = 3e8;                    % speed of light
K = B/T;                    % chirp slope
Rwid = Rmax - Rmin;         % recieved window in meter
Twid = 2 * Rwid / C;        % recieced window in second
Fs = 5 * B;                 % sampling frequency
Ts = 1 / Fs;                % sampling spacing
Nwid = ceil(Twid / Ts);     % recieved window in sampling number

1.2 生成回波信号

%% generate echo wave
t = linspace(2*Rmin/C, 2*Rmax/C, Nwid);         % received window
                                                % open window at 2Rmin/C
                                                % close window at 2Rmax/C

M = length(R);                                  % number of targets
delay = 2*R/C;                                  % delay of point targets
td = ones(M, 1)*t - delay'*ones(1, Nwid);       % time -tau of targets
Srt = RCS*(exp(1i*pi*K*td.^2).*(abs(td)<T/2));  % echo from all targets

第二行代码生成了一个 $\frac{2\times R_{min}}{c}$ 到 $\frac{2\times R_{max}}{c}$ 的时间序列（观测窗口），间隔为 $T_S$，大小为 $1\times N_{wid}$

第八行代码中，ones(M, 1)*t 得到一个 $M\times N_{wid}$ 的矩阵，每一行代表一个时间序列，delay'*ones(1, Nwid) 为时间延迟，故第八行的每一列相当于 $t-dela{y}_i$ 的时间序列

每一个回波信号的数学表达式为：
$$\begin{gathered} t{d}_i=t-2\times \frac{R_i}{c} \\ {S}_{r_i}(t)=S(t{d}_i)\ast Rect(\frac{t{d}_i}{T}) \end{gathered}$$
其中，$2\times \frac{R_i}{c}$ 为第 i 个目标的距离引起的时间延迟，这个延迟将反映在回波信号的延迟上，$S_{r_i}(t)$ 乘上 $Rect(\frac{t{d}_i}{T})$ 的原因是发射信号是脉冲信号，每次发射一个脉冲宽度

最后一行代码为计算回波信号，根据上面的分析就能理解最后一行代码了（最后一行中，第 i 行乘上第 i 个目标的 RCS后求和，矩阵运算），如下图所示：

---

2. 脉冲压缩

使用频域相关计算回波信号 $S_{rt}$ 与发射 $S_t$ 做相关运算计算匹配滤波（也可以用发射信号的翻转共轭与回波卷积）：

Sot = fftshift(ifft(fft(Srt, Nfft).*conj(fft(St, Nfft))));		% corr
Sot = conv(Srt, conj(fliplr(St)));		% conv

---

2.1 未加窗

%% Pulse compression radar using FFT and IFFT
Nfft = 2^nextpow2(Nwid + Nwid-1);       % number needed to compute linear
                                        % convolve length = Nwid + Nwid-1

Nchirp = ceil(T / Ts);                  % pulse duration in number
t0 = linspace(-T/2, T/2, Nchirp);       % discrete time of T

St = exp(1i*pi*K*t0.^2);                % transmit signal
Sot = fftshift(ifft(fft(Srt, Nfft).*conj(fft(St, Nfft))));

N0 = Nfft/2 - Nchirp/2;                 % valid series of fft in positive axis
Z = abs(Sot(N0 : N0+Nwid-1));           % FFT series
max_Z = max(Z);
Z = Z / max(Z);
Z = 20*log10(Z + 1e-6);

注：倒数第四行代码还不理解，如果有大佬理解可以留言哈！

---

2.2 加 hann 窗

将发射信号 $S_t$ 乘上 hann 窗，再与回波信号做相关运算，实现加窗的脉冲压缩（也可以先计算发射信号与窗函数相乘后，求 $S_t\times win$ 的共轭翻转，得到匹配滤波的脉冲响应函数，将脉冲响应函数与回波做卷积）

%% hann window
St = exp(1i*pi*K*t0.^2);                % transmit signal
win = hann(Nchirp)';
St_w = St.*win;
Sot = fftshift(ifft(fft(Srt, Nfft).*conj(fft(St_w, Nfft))));

N0 = Nfft/2 - Nchirp/2;                 % valid series of fft in positive axis
Z_w = abs(Sot(N0 : N0+Nwid-1));           % FFT series
Z_w = Z_w / max_Z;
Z_w = 20*log10(Z_w + 1e-6);

---

3. 绘制图像

%% Figures
figure(2)
subplot(3, 1, 1)
plot(t*1e6, real(Srt), 'black', 'LineWidth', 0.8); axis tight;
xlabel('(1) Time in \itu sec');
ylabel('Amplitude');
title('Radar echo without compression');

subplot(3, 1, 2)
plot(t*C/2, Z, 'black', 'LineWidth', 0.8)
axis([10000, 15000, -100, 0])
xlabel('(2) Range in meters');
ylabel('Amplitude in dB');
title('Radar echo after compression');

subplot(3, 1, 3)
plot(t*C/2, Z_w, 'black', 'LineWidth', 0.8)
axis([10000, 15000, -100, 0])
xlabel('(3) Range in meters');
ylabel('Amplitude in dB');
title('Radar echo after compression (windowed)');

输出结果如下：

---

完整程序：

function LFM_radar(T, B, Rmin, Rmax, R, RCS)
% inputs:
%   T: pulse duration
%   B: chirp frequency modulation bandwidth
%   Rmin, Rmax: range bin
%   R: position of ideal point targets
%   RCS: radar cross section of ideal targets

clear; clc;
set(0,'defaultfigurecolor', 'w');

%% default parameters
if nargin == 0
    T = 10e-6;
    B = 30e6;
    Rmin = 10000;
    Rmax = 15000;
    R = [11000, 12000, 12005, 13000, 13020];
    RCS = [1 1 1 1 1];
end

%% parameters
C = 3e8;                    % speed of light
K = B/T;                    % chirp slope
Rwid = Rmax - Rmin;         % recieved window in meter
Twid = 2 * Rwid / C;        % recieced window in second
Fs = 5 * B;                 % sampling frequency
Ts = 1 / Fs;                % sampling spacing
Nwid = ceil(Twid / Ts);     % recieved window in sampling number


%% generate echo wave
t = linspace(2*Rmin/C, 2*Rmax/C, Nwid);         % received window
                                                % open window at 2Rmin/C
                                                % close window at 2Rmax/C
M = length(R);                                  % number of targets
delay = 2*R/C;                                  % delay of point targets
td = ones(M, 1)*t - delay'*ones(1, Nwid);       % time -tau of targets
Srt = RCS*(exp(1i*pi*K*td.^2).*(abs(td)<T/2));  % echo from all targets

%% Pulse compression radar using FFT and IFFT
Nfft = 2^nextpow2(Nwid + Nwid-1);       % number needed to compute linear
                                        % convolve length = Nwid + Nwid-1
Srw = fft(Srt, Nfft);                   % FFT of echo wave

Nchirp = ceil(T / Ts);                  % pulse duration in number
t0 = linspace(-T/2, T/2, Nchirp);       % discrete time of T
St = exp(1i*pi*K*t0.^2);                % transmit signal



%% not windowed
Sw = fft(St, Nfft);
Sot = fftshift(ifft(Srw.*conj(Sw)));    % convolve using FFT

% Sot = fftshift(ifft(fft(Srt, Nfft).*conj(fft(St, Nfft))));    % convolve using FFT

N0 = Nfft/2 - Nchirp/2;                 % valid series of fft in positive axis
Z = abs(Sot(N0 : N0+Nwid-1));           % FFT series
max_Z = max(Z);
Z = Z / max(Z);
Z = 20*log10(Z + 1e-6);

%% hann window
St = exp(1i*pi*K*t0.^2);                % transmit signal
win = hann(Nchirp)';
St_w = St.*win;
% corelation in frequency domain
% fftshift(ifft(fft(X).*conj(fft(Y))))
Sot = fftshift(ifft(fft(Srt, Nfft).*conj(fft(St_w, Nfft))));


% different
sum(abs(conj(fft(St_w, Nfft))-fft(conj(fliplr(St_w)), Nfft)))
% Sot = fftshift(ifft(fft(Srt, Nfft).*fft(conj(fliplr(St_w)), Nfft)));  bad


N0 = Nfft/2 - Nchirp/2;                 % valid series of fft in positive axis
Z_w = abs(Sot(N0 : N0+Nwid-1));           % FFT series
Z_w = Z_w / max_Z;
Z_w = 20*log10(Z_w + 1e-6);

%% Figures
figure(2)
subplot(3, 1, 1)
plot(t*1e6, real(Srt), 'black', 'LineWidth', 0.8); axis tight;
xlabel('(1) Time in \itu sec');
ylabel('Amplitude');
title('Radar echo without compression');

subplot(3, 1, 2)
plot(t*C/2, Z, 'black', 'LineWidth', 0.8)
axis([10000, 15000, -100, 0])
xlabel('(2) Range in meters');
ylabel('Amplitude in dB');
title('Radar echo after compression');

subplot(3, 1, 3)
plot(t*C/2, Z_w, 'black', 'LineWidth', 0.8)
axis([10000, 15000, -100, 0])
xlabel('(3) Range in meters');
ylabel('Amplitude in dB');
title('Radar echo after compression (windowed)');