文章目录
脉冲压缩能够将信号压缩变窄，使得在时域上有重叠的两个信号能够分离，使用脉冲压缩和加窗抑制旁瓣，能够实现多目标的检测。
- 多目标信号的回波
使用 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
第二行代码生成了一个 2 × R m i n c \frac{2\times R_{min}}{c}c 2 ×R m i n 到 2 × R m a x c \frac{2\times R_{max}}{c}c 2 ×R m a x 的时间序列（观测窗口），间隔为 T S T_S T S ，大小为 1 × N w i d 1 \times N_{wid}1 ×N w i d
第八行代码中， ones(M, 1)*t 得到一个 M × N w i d M \times N_{wid}M ×N w i d 的矩阵，每一行代表一个时间序列， delay'*ones(1, Nwid) 为时间延迟，故第八行的每一列相当于 t − d e l a y i t – delay_i t −d e l a y i 的时间序列
每一个回波信号的数学表达式为：
t d i = t − 2 × R i c S r i ( t ) = S ( t d i ) ∗ R e c t ( t d i T ) td_i = t-2 \times \frac{R_i}{c}\ S_{r_i}(t) = S(td_i) * Rect(\frac{td_i}{T})t d i =t −2 ×c R i S r i (t )=S (t d i )∗R e c t (T t d i )
其中，2 × R i c 2 \times \frac{R_i}{c}2 ×c R i 为第 i 个目标的距离引起的时间延迟，这个延迟将反映在回波信号的延迟上，S r i ( t ) S_{r_i}(t)S r i (t ) 乘上 R e c t ( t d i T ) Rect(\frac{td_i}{T})R e c t (T t d i ) 的原因是发射信号是脉冲信号，每次发射一个脉冲宽度
最后一行代码为计算回波信号，根据上面的分析就能理解最后一行代码了（最后一行中，第 i 行乘上第 i 个目标的 RCS后求和，矩阵运算），如下图所示：
- 脉冲压缩
使用频域相关计算回波信号 S r t S_{rt}S r t 与发射 S t S_t 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 S_t S t 乘上 hann 窗，再与回波信号做相关运算，实现加窗的脉冲压缩（也可以先计算发射信号与窗函数相乘后，求 S t × w i n S_t \times win S t ×w i n 的共轭翻转，得到匹配滤波的脉冲响应函数，将脉冲响应函数与回波做卷积）
%% 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);
- 绘制图像
%% 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)');
Original: https://blog.csdn.net/qq_41140138/article/details/124913440
Author: 我有两颗糖
Title: 雷达成像 Matlab 仿真 3 —— 多目标检测
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/629172/
转载文章受原作者版权保护。转载请注明原作者出处！