【限免】线性调频信号的脉冲压缩及二维分离SAR成像算法【附MATLAB代码】

文章来源：微信公众号：EW Frontier

QQ交流群：949444104

程序一

对线性调频信号进行仿真，输出其时频域的相关信息，并模拟回波信号，

对其进行脉冲压缩和加窗处理。

实验记录：

1.线性调频信号时域包络、相位；实部、虚部

2.线性调频信号频谱幅频、相频特性；实部、虚部

3.两个目标回波的时域和频域波形

4.信号通过匹配滤波器的输出结果（脉冲压缩）。

5.用Hamming窗抑制脉冲压缩结果副瓣

%%  基本参数 clc;clear all;close all;
T = 10e-6; % LFM周期/脉宽 10usB = 60e6; % LFM带宽 60Mhzfs = 100e6; % 采样率 100MHzK = B/T;
%%  模拟发射信号n = round(15*T*fs);t = linspace(-10*T, 10*T,n);
lfmT = rectpuls(t,T).*exp(1j*pi*K*t.^2);lfmF = fftshift(fft(fftshift(lfmT)));f = linspace(-fs,fs,n);
%% 时域绘图figure();plot(diff(phase(lfmT)));title('LFM信号的时间-频率变化趋势图');xlabel('时间');ylabel('频率');xlim([7200,7800])% 包络figure();subplot(2,2,1);plot(t,abs(lfmT));title('LFM信号时域包络');xlabel('t/s');ylabel('幅度');xlim([-1e-5,1e-5])ylim([-0.5,1.5])% 相位subplot(2,2,2);plot(t,phase(lfmT));title('LFM信号时域相位');xlabel('t/s');ylabel('相位');xlim([-5e-6,5e-6])% 实部subplot(2,2,3);plot(t,real(lfmT));title('LFM信号时域实部');xlabel('t/s');ylabel('幅度');xlim([-1.5e-6,1.5e-6]);ylim([-1,1]);% 虚部subplot(2,2,4);plot(t,imag(lfmT));title('LFM信号时域虚部');xlabel('t/s');ylabel('幅度');xlim([-1.5e-6,1.5e-6]);ylim([-1,1]);%% 频域绘图figure();subplot(2,2,1);plot(f,abs(lfmF));title('LFM信号幅频特性');xlabel('Hz');ylabel('幅度');
subplot(2,2,2);plot(unwrap(angle(lfmF)));title('LFM信号相频特性');xlabel('Hz');ylabel('相位');
subplot(2,2,3);plot(f,real(lfmF));title('LFM信号频谱实部');xlabel('Hz');ylabel('幅度');    xlim([-3e7,3e7]);
subplot(2,2,4);plot(f,imag(lfmF));title('LFM信号频谱虚部');xlabel('Hz');ylabel('幅度');    xlim([-3e7,3e7]);%%  模拟回波信号
% 延迟为50us和51us的两个信号( 5000点和5100点,t从4000开始存储 )%延迟t1=50e-6;t2=51e-6;
%由于实际信号时间从0开始,时间轴重新定义echo1=rectpuls((t-t1),T).*exp(1j*pi*K*(t-t1).^2);echo2=rectpuls((t-t2),T).*exp(1j*pi*K*(t-t2).^2);echo=echo1+echo2;
figure();subplot(2,2,1)plot(t,abs(echo1),'r');hold on;plot(t,abs(echo2),'b');title('两个回波信号');xlabel('t/s');ylabel('幅度');xlim([3.5e-5,7e-5])ylim([0,1.5])
subplot(2,2,2)plot(t,real(echo));title('回波信号时域实部');xlabel('t/s');ylabel('幅度');xlim([4.9e-5 5.2e-5])
subplot(2,2,3)plot(f,fftshift(abs(fft(echo))));title('回波信号幅频特性');xlabel('Hz');ylabel('幅度');
subplot(2,2,4)plot(t,imag(echo));title('回波信号时域虚部');xlabel('t/s');ylabel('幅度');xlim([4.9e-5 5.2e-5])
lfmT = rectpuls(t,T).*exp(1j*pi*K*t.^2);lfmF = fftshift(fft(fftshift(lfmT)));f = linspace(-fs,fs,n);

脉冲压缩(匹配滤波）

方法1：

把信号变到频域（2048点的FFT)[-fs/2,fs/2]

频域相乘 H(w) = *S(w)

逆傅里叶变换

方法2：

驻留相位原理获得传递函数

% 这里采用方法1%回波信号的频域形式echo_F=fftshift(fft(echo));%回波信号的匹配滤波器Hf=fftshift(fft(conj(lfmT)));%脉压结果Pc_F=echo_F.*Hf;%脉压时域结果Pc_T=ifftshift(ifft(Pc_F));

%用Hamming窗抑制副瓣%制作汉明窗Hm = [zeros(1,2300) hamming(9900)' zeros(1,2800)];%频域加窗Pcw_F = Hm.*Pc_F;Pcw_T=ifftshift(ifft(Pcw_F));
figure();subplot(2,2,1)plot(t,abs(Pc_T));title('脉冲压缩后的时域波形');xlim([4.5e-5 5.5e-5])
subplot(2,2,2)plot(f,abs(Pc_F));title('脉冲压缩后的频谱');
subplot(2,2,3)plot(f,abs(Pcw_T));title('脉压加窗后的时域波形');xlim([4.5e7 5.5e7])
subplot(2,2,4)plot(f,abs(Pcw_F));title('脉压加窗后的频谱');

程序一仿真结果

程序二

先根据点目标分布，计算出对应的延时，再根据表达式计算回波数据进行仿真。再对回波数据

进行距离向和方位向上的脉冲压缩，得出二维图像。

1.二维回波信号幅度、相位

2.距离向脉冲压缩结果的二维等高线图

3.点目标成像结果；二位等高线图，距离和方位剖面

4.用Hamming窗抑制成像结果副瓣

%%  基本参数和配置 
clc;clear all;close all;v_c =3e+8;%光速T =10e-6;%发射脉冲时间Br=60e6;%距离向带宽lamda=0.03;%波长f0=v_c/lamda;%载频vx=150;%雷达平台运动速度R0=15e3;%场景中心最短斜距Kr=Br/T;%调频斜率Nr=2048;%距离向采样点数,必须要大于T*BrFr=100e6;%距离向采样频率deta_t=1/Fr;%距离向采样时间间隔tr=2*R0/v_c+((0:Nr-1)-Nr/2)*deta_t;%距离向采样时间轴fr=((-Nr/2):(Nr/2-1))/Nr*Fr;PRF=100;%PRFNa=60;%方位向采样点数ta=((0:(Na-1))-Na/2)/PRF;T_sar=Na/PRF;%方位向采样时间fa=((-Na/2):(Na/2-1))/Na*PRF;Ka=2*vx^2/(lamda*R0);

%% 根据目标点计算二维回波信号%%计算放置点的参数dot_num_a=1; % 方位向点个数deta_a=100; % 方位向点间距dot_num_r=3; % 距离向点个数deta_r=600; % 距离向点间距dot_xy_cell=cell(1,dot_num_a);
middle_point_r=ceil(dot_num_r/2);middle_point_a=ceil(dot_num_a/2);line_x=vx*ta;line_y=zeros(1,Na);
for i_dot_num_a=1:dot_num_adot_xy=zeros(dot_num_r,2);for i_dot_num_r=1:dot_num_rdot_xy(i_dot_num_r,2)=(i_dot_num_r-middle_point_r)*deta_r;dot_xy(i_dot_num_r,1)=(i_dot_num_a-middle_point_a)*deta_a;enddot_xy_cell{1,i_dot_num_a}=dot_xy;end
slant_range_cell=cell(1,dot_num_a);%计算每个点在所有方位的斜距for i_dot_num_a=1:dot_num_aslant_range=zeros(dot_num_r,Na);%singledot_xy=dot_xy_cell{1,i_dot_num_a};for i_dot_num_r=1:dot_num_rslant_range(i_dot_num_r,:)=sqrt((line_y-(R0+dot_xy(i_dot_num_r,2))).^2+(line_x-dot_xy(i_dot_num_r,1)).^2);%???endslant_range_cell{1,i_dot_num_a}=slant_range;end%计算每个点在所有方位的时延t_delay_cell=cell(1,dot_num_a);for i_dot_num_a=1:dot_num_aslant_range=slant_range_cell{1,i_dot_num_a};t_delay=slant_range*2/v_c;%singlet_delay_cell{1,i_dot_num_a}=t_delay;endecho_data=zeros(Nr,Na);tr_matrix=repmat(tr.',1,Na);%矩阵化处理for i_dot_num_a=1:dot_num_at_delay= t_delay_cell{1,i_dot_num_a};for i_dot_num_r=1:dot_num_rt_delay_matrix=repmat(t_delay(i_dot_num_r,:),Nr,1);%矩阵化处理echo_data=echo_data+rect((tr_matrix-t_delay_matrix)/T).*exp(1j*pi*Kr*(tr_matrix-t_delay_matrix).^2).*exp(-1j*2*pi*f0*t_delay_matrix);endendecho_data_phase=unwrap(angle(echo_data));
% 输出图像figure;subplot(2,1,1)imagesc(abs(echo_data));title('二维信号的幅度');subplot(2,1,2)imagesc(echo_data_phase);title('二维信号的相位');%% 脉冲压缩%距离向脉冲压缩Echo=circshift(fft2(circshift(echo_data,[-Nr/2,-Na/2])),[Nr/2,Na/2]);ref=exp(1j*pi*fr.^2/Kr);ref_matrix=repmat(ref.',1,Na);% 距离向二维匹配滤波器COMP=Echo.*ref_matrix; % 距离向匹配滤波后频域
%方位向脉冲压缩fr_matrix=repmat(fr.',1,Na);fa_matrix=repmat(fa,Nr,1);Haz=exp(-1j*pi.*(fa.^2)./Ka);Haz_matrix=repmat(Haz,Nr,1); % 方位向二维匹配滤波器SAC=COMP.*Haz_matrix; % 再经方位向脉冲压缩后

%     [a b]=find(SAC==max(max(SAC))); % 找出最大值
%     figure;plot(abs(SAC(:,31)));
%     figure;plot(abs(SAC(1599,:)));sac=circshift(ifft2(circshift(SAC,[0,0])),[Nr/2,Na/2]); % 时域图像figure;subplot(2,1,1)imagesc(abs(sac));title('点目标成像时域结果')subplot(2,1,2)imagesc(abs(SAC));title('点目标成像频域结果')
%% 生成Hamming窗window=hamming(1151)*hamming(37).';window_hamming=zeros(2048,60);xx=10;window_hamming(450+xx:1600+xx,13:49)=window;WINDOW_data=window_hamming.*SAC;window_data=circshift(ifft2(circshift(WINDOW_data,[0,0])),[Nr/2,Na/2]);figure;subplot(2,1,1)imagesc(abs(window_data));title('加Hamming窗后时域')subplot(2,1,2)imagesc(abs(WINDOW_data));title('加Hamming窗后频域')
%% 剖面图i_x=1; % 第i个点i_y=1;
x=1025+300*(i_x-1);y=31+400*(i_y-1);pou=sac(x-10:x+9,y-10:y+9); % 截取pou_range=pou(1:20,11);pou_azimuth=pou(11,1:20);range_db=ifft([fft(pou_range);zeros(980,1)]);tr_db=linspace(tr(x-10),tr(x+10),1000);
azimuth_db=ifft([fft((pou_azimuth).');zeros(980,1)]);ta_db=linspace(ta(y-10),ta(y+10),1000);
POU=fft2(pou);POU1=zeros(1000,1000);POU1(1:20,1:20)=POU;pou1=ifft2(POU1);a=max(abs(pou1));b=max(a);pou1_db=20*log10(abs(pou1)/b);c=max(pou1_db);dgx1=[-3 -13 -20 -30];figure;subplot(2,1,1)contour(pou1_db,dgx1);title('-3 -13 -20 -30(db)不加窗时二维图像');xlabel('方位向(m)');ylabel('距离向(m)');
% 加窗pou2=window_data(x-10:x+9,y-10:y+9);
pou_range2=pou2(1:20,11);pou_azimuth2=pou2(11,1:20);range_db2=ifft([fft(pou_range2);zeros(980,1)]);tr_db=linspace(tr(x-10),tr(x+10),1000);
azimuth_db2=ifft([fft((pou_azimuth2).');zeros(980,1)]);ta_db=linspace(ta(y-10),ta(y+10),1000);
POU=fft2(pou2);POU1=zeros(1000,1000);POU1(1:20,1:20)=POU;pou1=ifft2(POU1);a=max(abs(pou1));b=max(a);pou1_db=20*log10(abs(pou1)/b);c=max(pou1_db);dgx2=[-3 -13 -20 -30];subplot(2,1,2)contour(pou1_db,dgx2);title('-3 -13 -20 -30(db)加窗时二维图像');xlabel('方位向(m)');ylabel('距离向(m)');

% 加窗对比输出figure;subplot(2,2,1)plot(tr_db,20*log10(abs(range_db)/max(abs(range_db))));title('不加窗的距离像')subplot(2,2,2)plot(ta_db,20*log10(abs(azimuth_db)/max(abs(azimuth_db))));title('不加窗的方位像')subplot(2,2,3)plot(tr_db,20*log10(abs(range_db2)/max(abs(range_db2))));title('加窗的距离像')subplot(2,2,4)plot(ta_db,20*log10(abs(azimuth_db2)/max(abs(azimuth_db2))));title('加窗的方位像')

程序二仿真结果