雷达成像 Matlab 仿真 2 —— 脉冲压缩与加窗

1. 脉冲压缩

1.1 原理

脉冲压缩的实质是匹配滤波，匹配滤波器的系统函数 $H(j\omega )$ 与输入信号的频谱特性 $S_{i}(j\omega) $ 的关系为：
$$H(j\omega )=k\times S_i^{\ast }(j\omega )e^{-j\omega t_0}$$

幅频特性与输入信号一致
$$|H(j\omega )|=|S_i^{\ast }(j\omega )|$$

相位特性与输入信号相反，同时叠加了时延 $t_0$
$$\varphi (j\omega )=-[\theta (\omega )+\omega t_0]=-[-\frac{\pi }{K}{f}^2+\frac{\pi }{4}]-\omega t_0$$

这样，匹配滤波能够使 输出信号在采样时刻 $t_0$ 处达到最大，相位特性首先使高频成分的延迟更小，低频成分延迟更大，这样使得高频和低频成分集中在一块叠加，宏观上表现为峰值增大、脉冲变窄；相位特性第二部分 $-\omega t_0$ 引入了延迟，让叠加后的峰值在 $t_0$ 时刻达到最大，便于采样。

匹配滤波器的时域冲激响应（ 时域共轭翻转，频域共轭）：
$$h(t)=k\times s_i^{\ast }(t_0-t)$$

仿真时可以设置 $t_0=0$，matlab中 $h(t)$ 可任意表示为 conj(flip(st))

脉冲压缩输入信号宽度为 $T$，调频斜率为 $K$ 的 LFM 信号，输出信号近似 $sinc$ 函数，主瓣宽度 $\tau$ 输出信号宽度为 $\frac{1}{B}$（主瓣的两个零点分别为 $-\frac{1}{B}$ 和 $\frac{1}{B}$，取 3dB 位置处的宽度为 $\frac{0.886}{B}$ ，近似为 $\tau =\frac{1}{B}$），脉冲压缩比$D$ 是输入脉冲宽度和输出脉冲宽度的比值：
$$D=T\times B$$

---

1.2 程序实现脉冲压缩

1. 时域卷积

t = linspace(-T/2, T/2, N);
Si = exp(1i*pi*K*t.^2);     % chirp signal
Ht = exp(-1i*pi*K*t.^2);    % matched filter   s*(-t)
Sot = conv(Si, Ht);         % chirp signal after matched filter

2. 频域相关

匹配滤波相当于回波信号与发射信号做相关运算，参考

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

---

完整程序

clear; clc;							% 清空变量和输出
set(0,'defaultfigurecolor', 'w')	% 设置图片白底背景

生成输入信号 $S_i$

%% parameters
T = 10e-6;      % Pulse duration 10us
B = 30e6;       % Bandwidth 30MHz
K = B / T;      % chirp slope
Fs = 10 * B;    % sampling frequency
Ts = 1 / Fs;    % sampling spacing
N = T / Ts;     % Number of samples

%% Chirp after matched filter 
t = linspace(-T/2, T/2, N);
Si = exp(1i*pi*K*t.^2);     % chirp signal

通过卷积函数 conv(X, Y) 进行脉冲压缩

Ht = exp(-1i*pi*K*t.^2);    % matched filter   s*(-t)
Sot = conv(Si, Ht);         % chirp signal after matched filter

将匹配滤波输出转换为用 $dB$ 做单位，绘制波形

L = 2 * N - 1;
t1 = linspace(-T, T, L);

Z = abs(Sot);               % Amplitude of signal
Z = Z / max(Z);             % normalize
Z = 20*log10(Z + 1e-6);     % dB

subplot(2, 1, 1)
plot(t1 * B, Z, 'black', 'LineWidth', 1);
axis([-15, 15, -50, inf]);
grid on;
xlabel('Time in sec \times \itB');
ylabel('Amplitude, dB');
title('Chirp signal after matched filter');


%% Zoom
N0 = 3 * Fs / B;
t2 = -N0*Ts : Ts : N0*Ts;

subplot(2, 1, 2)
plot(t2 * B, Z(N-N0 : N+N0), 'black', 'LineWidth', 1);
grid on;
set(gca, 'Xtick', [-3, -2, -1, -0.5, 0, 0.5, 1, 2, 3]);
set(gca, 'Ytick', [-50, -13.4, -4, 0]);
xlabel('Time in sec \times \itB');
ylabel('Amplitude, dB');
title('Chirp signam after matched filter (zoom)');

结果：得到图像如下：

分析：输入信号的时域宽度为 $10\mu s$，输出信号宽度为 $1/30MHz$，故脉压比 $D=10\mu s\times 30MHz=300$，也就是说信号经过匹配滤波器后，宽度变为原来的 $1/300$！

注：频域实现脉冲压缩的例子参考 脉冲压缩检测多目标

---

2. 加窗

加窗是为了抑制旁瓣，防止旁瓣过高出现伪距目标甚至遮挡真实目标。

加窗的本质为对匹配滤波器的脉冲响应函数加权

不加窗前，匹配滤波器的时域脉冲响应为：发射信号的共轭翻转；
加窗后，匹配滤波器的时域脉冲响应为：发射信号加窗后，再取共轭翻转。

---

时域卷积：
使用 conv 函数，对回波信号 $S_{rt}$ 和加窗后的匹配滤波器脉冲响应函数 conj(fliplr(St.*win')) 做卷积运算

win = hann(Ns);
Ht_w = conj(fliplr(St.*win'));

L = 2 * Ns - 1;             % length of signal after convolve
t1 = linspace(-T, T, L);
Sot_w = conv(St, Ht_w);     % using convolve

上面的 hann(Ns) 为长度为 Ns 的汉宁窗，默认为列向量！

---

频域卷积

win = blackman(Ns);
Ht_w = conj(fliplr(St.*win'));      		% matched filter impulse response
L = 2 * Ns - 1;                     		% length of signal after convolve
Sot_w = ifft(fft(St, L).*fft(Ht_w, L));		% convolve in frequency domain

上面的 blackman(Ns) 为长度为 Ns 的汉宁窗，默认为列向量！

---

相关实现
让回波信号 $S_r(t)$ 与发射信号 $S_t(t)\times w(t)$ 做相关运算，直接得到输出（频域实现）

win = hann(length_St)';
St_w = St.*win;

% corelation in frequency domain
Sot = fftshift(ifft(fft(Srt, Nfft).*conj(fft(St_w, Nfft))));

注意互相关运算的输出与序列的先后顺序有关，下面两种写法不等价：

cor_xy = fftshift(ifft(fft(X).*conj(fft(Y))))
cor_yx = fftshift(ifft(fft(Y).*conj(fft(X))))

---

完整程序的输出结果如下图：

可见，加窗后能有效抑制旁瓣，但缺点是带来了主瓣的展宽（因为加窗后的匹配滤波器脉冲响应函数与发射信号的共轭翻转不完全匹配，因此主瓣变宽，不利于采样时刻对准），具体程序见文章底部。

---

3.完整程序

3.1 脉冲压缩完整程序

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

%% parameters
T = 10e-6;      % Pulse duration 10us
B = 30e6;       % Bandwidth 30MHz
K = B / T;      % chirp slope
Fs = 10 * B;    % sampling frequency
Ts = 1 / Fs;    % sampling spacing
N = T / Ts;     % Number of samples

%% Chirp after matched filter 
t = linspace(-T/2, T/2, N);
Si = exp(1i*pi*K*t.^2);     % chirp signal
Ht = exp(-1i*pi*K*t.^2);    % matched filter   s*(-t)
Sot = conv(Si, Ht);         % chirp signal after matched filter

L = 2 * N - 1;
t1 = linspace(-T, T, L);

Z = abs(Sot);               % Amplitude of signal
Z = Z / max(Z);             % normalize
Z = 20*log10(Z + 1e-6);     % dB

Z1 = abs(sinc(B.*t1));      % sinc function
Z1 = 20*log10(Z1 + 1e-6);   % dB

subplot(2, 1, 1)
plot(t1 * B, Z, 'black', 'LineWidth', 1); hold on;
plot(t1 * B, Z1, 'black.', 'MarkerSize', 12);
axis([-15, 15, -50, inf]);
grid on;
legend('emulational', 'sinc');
xlabel('Time in sec \times \itB');
ylabel('Amplitude, dB');
title('Chirp signal after matched filter');


%% Zoom
N0 = 3 * Fs / B;
t2 = -N0*Ts : Ts : N0*Ts;

subplot(2, 1, 2)
plot(t2 * B, Z(N-N0 : N+N0), 'black', 'LineWidth', 1); hold on;
plot(t2 * B, Z1(N-N0 : N+N0), 'black.', 'MarkerSize', 12);
axis([-inf, inf, -50, inf]); grid on;
set(gca, 'Xtick', [-3, -2, -1, -0.5, 0, 0.5, 1, 2, 3]);
set(gca, 'Ytick', [-50, -13.4, -4, 0]);
xlabel('Time in sec \times \itB');
ylabel('Amplitude, dB');
legend('emulational', 'sinc');
title('Chirp signam after matched filter (zoom)');

3.2 加窗的完整程序

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

%% parameters
T = 10e-6;                  % Pulse duration 10us
B = 30e6;                   % Bandwidth 30MHz
K = B / T;                  % chirp slope
Fs = 10 * B;                % sampling frequency
Ts = 1 / Fs;                % sampling spacing
Ns = T / Ts;                % Number of samples

%% Chirp after matched filter 
t = linspace(-T/2, T/2, Ns);
St = exp(1i*pi*K*t.^2);     % chirp signal
Ht = conj(fliplr(St));      % matched filter impulse response
Sot = conv(St, Ht);         % chirp signal after matched filter

Z = abs(Sot);               % Amplitude of signal
Z_max = max(Z);
Z = Z / max(Z);             % normalize by max(Z_w)
Z = 20*log10(Z + 1e-6);     % dB

%% plot chirp after matched filter
L = 2 * Ns - 1;
t1 = linspace(-T, T, L);

subplot(3, 2, 1)
plot(t1 * B, Z, 'black', 'LineWidth', 1); hold on;
axis([-15, 15, -50, inf]);
grid on;
xlabel('(1) Time in sec \times \itB');
ylabel('Amplitude, dB');
title('Chirp signal after matched filter');

N0 = 4 * Fs / B;
t2 = -N0*Ts : Ts : N0*Ts;

subplot(3, 2, 2)
plot(t2 * B, Z(Ns-N0 : Ns+N0), 'black', 'LineWidth', 1);
axis([-inf, inf, -50, inf]); grid on;
set(gca, 'Xtick', [-4 -3 -2 -1 -0.5 0 0.5 1 2 3 4]);
set(gca, 'Ytick', [-50 -13.4 -4 0]);
xlabel('(2) Time in sec \times \itB');
ylabel('Amplitude, dB');
title('Chirp signam after matched filter (zoom)');


%% Hann window (method 1)
win = hann(Ns);
Ht_w = conj(fliplr(St.*win'));

L = 2 * Ns - 1;             % length of signal after convolve
t1 = linspace(-T, T, L);
Sot_w = conv(St, Ht_w);     % using convolve


%% Hann window (method 2)
win = hann(Ns)';
St_w = St.*win;
% xcorr(X, Y) = fftshift(ifft(fft(X).*conj(fft(Y))))    corelation in frequency domain
Sot_w = fftshift(ifft(fft(St, L).*conj(fft(St_w, L))));


Z_hann_w = abs(Sot_w);                          % Amplitude of signal
Z_hann_w = Z_hann_w / Z_max;                    % normalize by using max(Z)
Z_hann_w = 20*log10(Z_hann_w + 1e-6);           % dB

subplot(3, 5, 6)
plot(t*1e6, hann(Ns), 'black', 'LineWidth', 1); 
grid on;
xlabel('(3) Time in usec');
ylabel('Amplitude');
title('Hanning window');

subplot(3, 5, [7, 8])
plot(t1 * B, Z, 'color', [0.7 0.6 0.3], 'LineWidth', 1); hold on;
plot(t1 * B, Z_hann_w, 'color', [0.5 0.7 0.9], 'LineWidth', 1);
axis([-300, 300, -120, inf]);
grid on;
legend('not windowed', 'windowed');
xlabel('(4) Time in sec \times \itB');
ylabel('Amplitude, dB');
title('Side lobe suppression with Hanning window');

% ZOOM
N0 = 4 * Fs / B;
t2 = -N0*Ts : Ts : N0*Ts;

subplot(3, 5, [9 10])
plot(t2 * B, Z(Ns-N0 : Ns+N0), 'color', [0.7 0.6 0.3], 'LineWidth', 1.5); hold on;
plot(t2 * B, Z_hann_w(Ns-N0 : Ns+N0), 'color', [0.5 0.7 0.9], 'LineWidth', 1.5);
grid on;
legend('not windowed', 'windowed');
xlabel('(5) Time in sec \times \itB');
ylabel('Amplitude, dB');
title('zoom');
axis([-inf, inf, -60, inf]); grid on;
set(gca, 'Xtick', [-4 -3 -2 -0.80 0 0.80 2 3 4]);
set(gca, 'Ytick', [-60 -37.6 -10 -6 0]);

%% Blackman window
win = blackman(Ns);
Ht_w = conj(fliplr(St.*win'));      % matched filter impulse response
L = 2 * Ns - 1;                     % length of signal after convolve
t1 = linspace(-T, T, L);
Sot_w = ifft(fft(St, L).*fft(Ht_w, L));         % convolve in frequency domain

Z_blackman_w = abs(Sot_w);                      % Amplitude of signal
Z_blackman_w = Z_blackman_w / Z_max;            % normalize by using max(Z)
Z_blackman_w = 20*log10(Z_blackman_w + 1e-6);	% dB


subplot(3, 5, 11)
plot(t*1e6, blackman(Ns), 'black', 'LineWidth', 1); 
grid on;
xlabel('(6) Time in usec');
ylabel('Amplitude');
title('Blackman window');

subplot(3, 5, [12 13])
plot(t1 * B, Z, 'color', [0.7 0.6 0.3], 'LineWidth', 1); hold on;
plot(t1 * B, Z_blackman_w, 'color', [0.5 0.7 0.9], 'LineWidth', 1);
axis([-300, 300, -120, inf]);
grid on;
legend('not windowed', 'windowed');
xlabel('(7) Time in sec \times \itB');
ylabel('Amplitude, dB');
title('Side lobe suppression with Blackman window');

% ZOOM
N0 = 4 * Fs / B;
t2 = -N0*Ts : Ts : N0*Ts;

subplot(3, 5, [14 15])
plot(t2 * B, Z(Ns-N0 : Ns+N0), 'color', [0.7 0.6 0.3], 'LineWidth', 1.5); hold on;
plot(t2 * B, Z_blackman_w(Ns-N0 : Ns+N0), 'color', [0.5 0.7 0.9], 'LineWidth', 1.5);
grid on;
legend('not windowed', 'windowed');
xlabel('(8) Time in sec \times \itB');
ylabel('Amplitude, dB');
title('zoom');
axis([-inf, inf, -100, inf]); grid on;
set(gca, 'Xtick', [-4 -3 -2 -0.94 0 0.94 2 3 4]);
set(gca, 'Ytick', [-100 -65.7 -11.53 -7.53 0]);