#脉冲压缩仿真及源码

之前做过脉压的仿真，最近温习一下
假设雷达发射线性调频信号，带宽17.2MHz,脉冲宽度42us
该程序主要解决了三个问题：

1. 距离雷达30km有一个静止点目标，实现标准脉压，画出滤波器幅度、相位响应，脉压后信号的幅度、相位

2. 画出加Hamming窗后的滤波器幅度、相位响应，以及加窗脉压后信号的幅度、相位。

3. 当目标发生运动、产生的多普勒频率为带宽的5%，10%。20%时，分别画出标准脉压后信号的幅度和相位

注意事项

1.频域匹配滤波器实现有三种方法（这里采用第2种方法）：

1)将时间反褶后的发射复制脉冲取复共轭，然后做补零FFT
2)发射复制脉冲做补零FFT，然后取复共轭（不用反褶）
3)根据设定的线型调频特性，直接在频域生成匹配滤波器

2.频域加窗时需要考虑弃置区的问题
因为st的频谱如下图所示

而频域滤波器的频谱如下图

均未经fftshift处理
所以匹配滤波器点乘后，做ifft需要ifftshift将峰值拉回中间，此时图像图下图所示

去除两边的多余数据，即可得到频域加窗脉压后的信号幅度
3.加入多普勒的信息后，总的结果如下图所示

4.其余的注意事项均写在程序中，给大家参考

具体代码实现

%%
W = 17.2e6; %带宽
T = 42e-6; %时宽
Fs = 10 * W;
N = T * Fs;
t = linspace(-T/2 ,T/2,N);
st = exp(1j * pi * W / T * t.^2);% 发射信号形式

R = 30e3; % 目标距离
c = 3e8; %光速
tao = 2 * R /c; %双程延时 200us     

sr = exp(1j * pi * W / T * (t-tao).^2);% 接收信号形式

h = exp(-1j * pi * W / T * fliplr(t.^2));% 匹配滤波器形式

% h1 = conj(fliplr(st));
figure
subplot(211)
plot(t,real(h))
xlabel('I(t)')
subplot(212)
plot(t,imag(h))
xlabel('Q(t)')
title('滤波器的幅度响应')
figure
plot(t,angle(h))
title('滤波器的相位响应')

% y_time = conv(sr,h);%
% 如果采用这种形式，由于回波时延远大于脉宽，
% 回波信号的相位与匹配滤波器相位没有可以匹配得上的，如果是无限长的信号当然可以
% 这个时候只能利用st与h的卷积然后人为挪动到应该在的位置
y_time = conv(st,h);% 时域做卷积 
figure
plot(linspace(-T,T,2*N-1)*1e6+200,20*log10(abs(y_time/max(y_time))))
hold on;plot([199 201], [-3 -3],'r')
axis([199 201 -30 1]);grid on;
xlabel('时间,us'),ylabel('功率，dB'),title('时域脉压后的幅度')
% 3dB带宽 7220-7228index区间 0.0465us

figure
plot(linspace(-T,T,2*N-1)*1e6+200,abs(angle(y_time/max(y_time))))
axis([199 201 -5 5]);grid on;
xlabel('时间,us'),ylabel('角度，rad'),title('时域脉压后的相位')

% 频域实现脉压
h_fre = fft(h,16384);%之所以取16384是要大于2N-1即可
st_fre = fft(st,16384);% 一定要注意做fft补零，时域补零，频域插值，如果不补零，那么
% 很有可能看不到副瓣-13dB，因为插值没有那么密
y_fre = st_fre .* h_fre;
y_time1 = ifft(y_fre);
y_time1_norm = y_time1/max(y_time1);% 进行归一化

figure
plot(linspace(-T,T,2*N-1)*1e6+200,20*log10(abs(y_time1_norm(1:2*N-1))))
% 如果作图不取前1:2N-1个点，第2N点后边的值，ifft后非常小，从时域也可以理解，两个信号
% 卷积后总共信号长度就是2N-1，从频域做也应该如此
axis([199 201 -30 1]);grid on;
xlabel('时间,us'),ylabel('功率，dB'),title('频域脉压后的幅度')
figure
plot(linspace(-T,T,2*N-1)*1e6+200,abs(angle(y_time1_norm(1:2*N-1))))
axis([199 201 -5 5]);grid on;
xlabel('时间,us'),ylabel('角度，rad'),title('频域脉压后的相位')

%% 时域加窗
hamm = hamming(N);
h_hamm = h .* hamm';% hamming窗加权的匹配滤波器

figure
subplot(211)
plot(t,real(h_hamm))
xlabel('I(t)')
subplot(212)
plot(t,imag(h_hamm))
xlabel('Q(t)')
title('加窗匹配滤波器的幅度响应')
figure
plot(t,angle(h_hamm))
title('加窗匹配滤波器的相位响应')
% 加窗与否对滤波器的相位无影响，因为滤波器是实数，点乘对相位没有影响

y_time_hamm = conv(st,h_hamm);% 时域做卷积 
figure
plot(linspace(-T,T,2*N-1)*1e6+200,20*log10(abs(y_time_hamm/max(y_time_hamm))))
hold on;plot([199 201], [-3 -3],'r')
axis([199 201 -70 1]);grid on;
xlabel('时间,us'),ylabel('功率，dB'),title('时域加窗脉压后的幅度')
% 3dB带宽对应 7218-7230   0.0698us
figure
plot(linspace(-T,T,2*N-1)*1e6+200,abs(angle(y_time_hamm/max(y_time_hamm))))
axis([199 201 -5 5]);grid on;
xlabel('时间,us'),ylabel('角度，rad'),title('时域加窗脉压后的相位')

%% 频域加窗
hamm = hamming(N);
% hamm = kaiser(N,2.5)
st_hamm = st .* hamm';
h_hamm_fre = conj(fft(st_hamm,16384));%figure;plot(abs(h_hamm_fre))
h_hamm1 = ifft(h_hamm_fre);
figure
subplot(211)
plot(t,real(h_hamm1(end-N+1:end)))%取后边长度为N的序列
xlabel('I(t)')
subplot(212)
plot(t,imag(h_hamm1(end-N+1:end)))
xlabel('Q(t)')
title('频域加窗匹配滤波器的幅度响应')
figure
plot(t,angle(h_hamm1(end-N+1:end)))
title('频域加窗匹配滤波器的相位响应')

y_fre_hamming = st_fre .* h_hamm_fre;

y_time2 = ifftshift(ifft(y_fre_hamming));
y_time2_norm = y_time2/max(y_time2);% 进行归一化
y_time2_norm1 = y_time2_norm((16384-2*N)/2+2:16384-(16384-2*N)/2);% 去除弃置区
% 作图真实数据是2N-1个点，通过观察总的图形，确认弃置区

figure
plot(linspace(-T,T,2*N-1)*1e6+200,20*log10(abs(y_time2_norm1)))
axis([199 201 -70 1]);grid on;
xlabel('时间,us'),ylabel('功率，dB'),title('频域加窗脉压后的幅度')
figure
plot(linspace(-T,T,2*N-1)*1e6+200,abs(angle(y_time2_norm1)))
axis([199 201 -5 5]);grid on;
xlabel('时间,us'),ylabel('角度，rad'),title('频域加窗脉压后的相位')

%% 加入多普勒分析
fd1 = 0.05 * W;  %多普勒频移 带宽的5%
% st_fd1 = exp(1j * pi * W / T * (fd1 * t + t.^2));%
% 这种建模方法是错的，W/T是线性调频率，是2次相位的系数
st_fd1 = exp(1j * pi *(fd1 * t + W / T * t.^2));% 加入多普勒效应的信号

y_time_fd1 = conv(st_fd1,h);% 时域做卷积 
figure
plot(linspace(-T,T,2*N-1)*1e6+200,20*log10(abs(y_time_fd1/max(y_time_fd1))))
% hold on;plot([197 199], [-3 -3],'r')
axis([195 200 -40 1]);grid on;
xlabel('时间,us'),ylabel('功率，dB'),title('时域脉压后的幅度，fd=0.05W')
% 3dB带宽 7220-7228index区间 0.0465us

figure
plot(linspace(-T,T,2*N-1)*1e6+200,abs(angle(y_time_fd1/max(y_time_fd1))))
axis([195 200 -5 5]);grid on;
xlabel('时间,us'),ylabel('角度，rad'),title('时域脉压后的相位，fd=0.05W')


fd2 = 0.1 * W;  %多普勒频移 带宽的5%
% st_fd1 = exp(1j * pi * W / T * (fd1 * t + t.^2));%
% 这种建模方法是错的，W/T是线性调频率，是2次相位的系数
st_fd2 = exp(1j * pi *(fd2 * t + W / T * t.^2));% 加入多普勒效应的信号

y_time_fd2 = conv(st_fd2,h);% 时域做卷积 
figure
plot(linspace(-T,T,2*N-1)*1e6+200,20*log10(abs(y_time_fd2/max(y_time_fd2))))
% hold on;plot([197 199], [-3 -3],'r')
axis([195 200 -40 1]);grid on;
xlabel('时间,us'),ylabel('功率，dB'),title('时域脉压后的幅度，fd=0.1W')
% 3dB带宽 7220-7228index区间 0.0465us

figure
plot(linspace(-T,T,2*N-1)*1e6+200,abs(angle(y_time_fd2/max(y_time_fd2))))
axis([195 200 -5 5]);grid on;
xlabel('时间,us'),ylabel('角度，rad'),title('时域脉压后的相位，fd=0.1W')


fd3 = 0.2 * W;  %多普勒频移 带宽的5%
% st_fd1 = exp(1j * pi * W / T * (fd1 * t + t.^2));%
% 这种建模方法是错的，W/T是线性调频率，是2次相位的系数
st_fd3 = exp(1j * pi *(fd3 * t + W / T * t.^2));% 加入多普勒效应的信号

y_time_fd3 = conv(st_fd3,h);% 时域做卷积 
figure
plot(linspace(-T,T,2*N-1)*1e6+200,20*log10(abs(y_time_fd3/max(y_time_fd3))))
% hold on;plot([197 199], [-3 -3],'r')
axis([195 200 -40 1]);grid on;
xlabel('时间,us'),ylabel('功率，dB'),title('时域脉压后的幅度，fd=0.2W')

figure
plot(linspace(-T,T,2*N-1)*1e6+200,abs(angle(y_time_fd3/max(y_time_fd3))))
axis([195 200 -5 5]);grid on;
xlabel('时间,us'),ylabel('角度，rad'),title('时域脉压后的相位，fd=0.2W')