sar3.m

经典算法--距离多普勒算法

%--------------------------------------------------------%
% 正侧视SAR点目标仿真 
% 点目标参数:方位向距离为0,斜距为50000,后向反射系数为1
%--------------------------------------------------------%
clc;clear;close all;
%--------------------------------------------------------%
%Parameter--constant   恒参
C=3e8;                 %指的是光速3e8
%Parameter--radar characteristics    具有雷达特征的参数
lambda=0.032;
v = 150;
Kv = 2;                %每走1米发射2个脉冲即2个/m
R0 = 50000;
D = 4;
Lsar = lambda*R0/D;    %合成孔径长度
Tsar = Lsar/v;         %合成孔径时间
%Parameter--slow-time domain:Azimuth   慢时间范围：方位向
Kd = -2*v^2/lambda/R0; %doppler frequency modulation rate    多普勒频率调制率
Bd = abs(Kd*Tsar);     %doppler frequency modulation bandwidth   多普勒带宽
PRF = Kv*v;            %pulse repetition frequency  脉冲重复频率
PRT = 1/PRF;           %pulse repetition time       脉冲重复周期
ds = PRT;              %sample spacing in slow-time domain   在慢时间范围内的取样间隔
N = ceil(Lsar/v/ds);   %朝正无穷方向舍入 
N = 2^nextpow2(N);     %for fft      nextpow2(x)----取大于并最接近x的2次幂.
sn = linspace(-Lsar/2/v,Lsar/2/v,N);
%discrete time array in slow-time domain   linspace(x0,x1,n)线性等分向量：n代表的是点的数目，即分成n-1等分，步长应当是(x1-x0)/(n-1)
PRT = Lsar/v/N;        %refresh    更新
PRF = 1/PRT;
ds = PRT;
%Parameter--fast-time domain:Range   慢时间范围：距离向
Tr = 5e-6;             %pulse duration 10us   脉冲持续时间
Br = 100e6;            %chirp frequency modulation bandwidth 100MHz   chirp信号频率调制带宽100MHz ,回波信号带宽为B=Kr*Tr
Kr = Br/Tr;            %chirp slope   chirp信号斜率
Fsr = 150e6;           %sample frequency in fast-time domain    在快时间范围内的取样频率
dt = 1/Fsr;            %sample spacing in fast-time domain      在快时间范围内的取样间隔
Rmin = R0;             %R的最小值
Rmax = sqrt(R0^2+(Lsar/2)^2);       %R的最大值
M = ceil(2*(Rmax-Rmin)/C/dt+Tr/dt); %sample number in fast-time domain 回波信号持续时间为2*(Rmax-Rmin)/C+Tr???????????????????????
M = 2^nextpow2(M);                  %for fft    基2-FFT 
tm = linspace(2*Rmin/C-Tr/2,2*Rmax/C+Tr/2,M);   %discrete time array in fast-time domain????????????????????????????????????
dt = (2*Rmax/C+Tr-2*Rmin/C)/M;                  %refresh    更新
Fsr = 1/dt;
%Parameter--resolution     分辨率参数
DY = C/2/Br;            %range resolution       距离向分辨率
DX = D/2;               %cross-range resolution 方位向分辨率
%Parameter--point targets  点目标参数
Ptarget = [0,R0,1];     %target position and RCS  点目标参数:方位向距离为0,斜距为50000,后向反射系数为1
disp('Parameters:')
disp('Sampling Rate in fast-time domain');disp(Fsr/Br)  %奈奎斯特采样率-----采样率必须高于最高信号频率的两倍，试验结果为 2.0469
disp('Sampling Number in fast-time domain');disp(M)     %距离向取样点数 
disp('Sampling Rate in slow-time domain');disp(PRF/Bd)  %奈奎斯特采样率-----采样率必须高于最高信号频率的两倍????????????????????????
disp('Sampling Number in slow-time domain');disp(N)     %方位向取样点数 
disp('Range Resolution');disp(DY)
disp('Azimuth Resolution');disp(DX) 
disp('Synthetic Aperture Length');disp(Lsar)
disp('Synthetic Aperture Time');disp(Tsar)
disp('Position of targets');disp(Ptarget)
%--------------------------------------------------------%
%--------------------------------------------------------%
% Generate the raw signal data
Srnm = zeros(N,M);      %全零数组   N为方位向采样点，M为距离向采样点
Pa = v*sn;              %sample position to target in slow-time domain  定位在慢时间段的样点位置
R = sqrt(R0^2+(Pa).^2); %sample slant range   抽样点（即约为发射脉冲的位置）距离点目标的距离
delay = 2*R/C;          %time delay in transmitted pulse   发射脉冲的时间延迟
Pr = ones(N,1)*tm-delay'*ones(1,M); %sample time for transmitted pulse
phase = j*pi*Kr*(Pr).^2-j*4*pi/lambda*(R'*ones(1,M));      %经过检波之后的连续回波信号抽样化后的矩阵序列的指数部分 
Srnm = Ptarget(1,3)*exp(phase).*(abs(Pr)<=Tr/2).*(abs(Pa'*ones(1,M))<=Lsar/2); 
%sig : target echo model 点目标回波原始数据（将原始连续回波信号抽样化后的数据，这样利于作基2-快速傅立叶变换）
%--------------------------------------------------------%
%--------------------------------------------------------%
% Range compression          %距离压缩
fr = -1/2:1/M:(1/2-1/M);     %M个抽样点将其化为M-1等分  
fr = fr*Fsr;                 %此处抽样频率约等于信号带宽（前面有结果=Fsr/Br= 2.0469）
filter_r = fftshift(exp(j*pi*fr.^2/Kr));     %matched filter in range derection在距离向的匹配滤波
for i=1:N
Srnm(i,:)=ifft(fft(Srnm(i,:)).*filter_r);
end
disp('End of range compression'); %距离向先进行FFT，后进行IFFT    
row = C*tm/2;                   %行？？？？？？？？？？？？？？？？？？？/？？？？？？？？？？？？？？？？？？？
col = v*sn;                     %列
colormap(gray);                 %colormap将当前颜色表置为Matlab内置的灰度图
figure(1);
imagesc(row,col,255-abs(Srnm)); %显示亮度图像，此处为灰度图象         255-abs(Srnm)？？？？？？？？？？？？？？？？？？？？？？？？
xlabel('\rightarrow\itRange in meters'),ylabel('\itAzimuth in meters\leftarrow'),
title('2-d time domain after range compression'),
%--------------------------------------------------------%
%--------------------------------------------------------%
% Azimuth fft 方位向的FFT
for i=1:M
Srnm(:,i)=fftshift(fft(Srnm(:,i)));
end
disp('End of azimuth fft');     %方位向的FFT
%--------------------------------------------------------%
%--------------------------------------------------------%
% RCMC
fa = -1/2:1/N:(1/2-1/N);        %N个抽样点将其化为N-1等分
fa = fa*PRF;                    %为什么要用 PRF而不用Bd ？？？？？？？？？？？？？？？？？？？？？？
deteR = 2/C*(lambda^2*R0*fa.^2/8/v^2)*Fsr;     % 距离徙动所占的采样点数,不一定为整数值,所以接下来要进行插值
N_interp = 8;                   %N_interp=8表示以8倍进行插值
N_ex = 1000;                    %N_ex必须大于距离徙动所占的采样点数
for i=1:N
newp = round(((1:M)+deteR(i)-1)*N_interp+1); %round表示四舍五入取整
Srnm_interp=zeros(size(1:M*N_interp+N_ex));
Srnm_interp(1:M*N_interp)=interp(Srnm(i,:),N_interp); %interp表示插值,详见help interp
Srnm(i,:)=Srnm_interp(newp);
end
disp('End of RCMC');
colormap(gray);
figure(2);
imagesc(row,col,255-abs(Srnm));
xlabel('\rightarrow\itRange in meters'),ylabel('\itAzimuth in meters\leftarrow'),
title('r-d domain after interpolation');%经过插值之后
%--------------------------------------------------------%
%--------------------------------------------------------%
% Azimuth compression  方位压缩
fa = -1/2:1/N:(1/2-1/N);
fa = fa'*PRF;      %注意要转置，为了下面第（4）步
filter_a = exp(-j*pi*fa.^2/abs(Kd));    %matched filter in azimuth derection方位向的匹配滤波
for i=1:M
Srnm(:,i)=ifft(fftshift(Srnm(:,i).*filter_a));%(4)
end
disp('End of azimuth compression');% 完成Azimuth  compression 方位压缩

colormap(gray);
figure(3);
imagesc(row,col,255-abs(Srnm));
xlabel('\rightarrow\itRange in meters'),ylabel('\itAzimuth in meters\leftarrow'),
title('2-d time domain after azimuth compression');
figure(4);
mesh(row(450:550),col(450:550),abs(Srnm(450:550,450:550)));axis tight;
xlabel('\rightarrow\itRange in meters'),ylabel('\itAzimuth in meters\leftarrow'),
title('2-d time domain after azimuth compression');