Matlab simulation of radar imaging 3 - multi-target detection

Pulse compression can narrow the signal , So that two signals with overlap in the time domain can be separated , Use pulse compression and windowing to suppress side lobes , It can realize multi-target detection .

1. Echo of multi-target signal

Use matlab when , If you can use matrix operation, try not to for loop , Because matrix operation can speed up the operation .
When there are multiple targets in a certain direction within the radar detection range , Look, you can use each row of the matrix to store the echo signal of each target , The echo signals of all lines are superimposed to get the actually received echo .
Next, write a pulse compression function to detect multiple targets

1.1 Parameter setting

First set the pulse width 、 bandwidth 、 Maximum and minimum distance of detection 、 Point the distance of the target （ vector ） And point target RCS. Then calculate some variable values that need to be used ：Twid, Fs,Ts and 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 signal

%% 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

The second line of code generates a $\frac{2\times R_{min}}{c}$ To $\frac{2\times R_{max}}{c}$ Time series of （ Observation window ）, The interval is $T_S$, The size is $1\times N_{wid}$

In the eighth line of code ,ones(M, 1)*t Get one $M\times N_{wid}$ Matrix , Each line represents a time series ,delay'*ones(1, Nwid) For time delay , Therefore, each column in the eighth row is equivalent to $t-dela{y}_i$ Time series of

The mathematical expression of each echo signal is ：
$$\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}$$
among ,$2\times \frac{R_i}{c}$ For the first time i The time delay caused by the distance between two targets , This delay will be reflected in the delay of the echo signal ,$S_{r_i}(t)$ Multiply $Rect(\frac{t{d}_i}{T})$ The reason is that the transmitted signal is a pulse signal , One pulse width at a time

The last line of code is to calculate the echo signal , According to the above analysis, we can understand the last line of code （ In the last line , The first i Line multiplied by No i A goal of RCS After the sum , Matrix operations ）, As shown in the figure below ：

---

2. Pulse compression

Use frequency domain correlation to calculate the echo signal $S_{rt}$ And launch $S_t$ Do correlation operation to calculate matched filter （ It can also be used to convolute the inverted conjugate of the transmitted signal with the echo ）：

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

---

2.1 Without windows

%% 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);

notes ： The penultimate line of code is not understood , If there are big guys understand, you can leave a message ！

---

2.2 Add hann window

Will send a signal $S_t$ Multiply hann window , Then do correlation calculation with the echo signal , Realize windowed pulse compression （ You can also calculate the transmission signal and multiply it by the window function first , seek $S_t\times win$ Conjugate flip of , Get the impulse response function of matched filter , Convolute the impulse response function with the echo ）

%% 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. The plot

%% 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)');

The output is as follows ：

---

The whole program ：

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)');