%--------------------------------------------------------%
% 正側(cè)視SAR點(diǎn)目標(biāo)仿真 
% 點(diǎn)目標(biāo)參數(shù):方位向距離為0,斜距為50000,后向反射系數(shù)為1
%--------------------------------------------------------%
clc;clear;close all;
%--------------------------------------------------------%
%Parameter--constant   恒參
C=3e8;                 %指的是光速3e8
%Parameter--radar characteristics    具有雷達(dá)特征的參數(shù)
lambda=0.032;
v = 150;
Kv = 2;                %每走1米發(fā)射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    多普勒頻率調(diào)制率
Bd = abs(Kd*Tsar);     %doppler frequency modulation bandwidth   多普勒帶寬
PRF = Kv*v;            %pulse repetition frequency  脈沖重復(fù)頻率
PRT = 1/PRF;           %pulse repetition time       脈沖重復(fù)周期
ds = PRT;              %sample spacing in slow-time domain   在慢時間范圍內(nèi)的取樣間隔
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代表的是點(diǎn)的數(shù)目，即分成n-1等分，步長應(yīng)當(dāng)是(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   脈沖持續(xù)時間
Br = 100e6;            %chirp frequency modulation bandwidth 100MHz   chirp信號頻率調(diào)制帶寬100MHz ,回波信號帶寬為B=Kr*Tr
Kr = Br/Tr;            %chirp slope   chirp信號斜率
Fsr = 150e6;           %sample frequency in fast-time domain    在快時間范圍內(nèi)的取樣頻率
dt = 1/Fsr;            %sample spacing in fast-time domain      在快時間范圍內(nèi)的取樣間隔
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 回波信號持續(xù)時間為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     分辨率參數(shù)
DY = C/2/Br;            %range resolution       距離向分辨率
DX = D/2;               %cross-range resolution 方位向分辨率
%Parameter--point targets  點(diǎn)目標(biāo)參數(shù)
Ptarget = [0,R0,1];     %target position and RCS  點(diǎn)目標(biāo)參數(shù):方位向距離為0,斜距為50000,后向反射系數(shù)為1
disp('Parameters:')
disp('Sampling Rate in fast-time domain');disp(Fsr/Br)  %奈奎斯特采樣率-----采樣率必須高于最高信號頻率的兩倍，試驗(yàn)結(jié)果為 2.0469
disp('Sampling Number in fast-time domain');disp(M)     %距離向取樣點(diǎn)數(shù) 
disp('Sampling Rate in slow-time domain');disp(PRF/Bd)  %奈奎斯特采樣率-----采樣率必須高于最高信號頻率的兩倍????????????????????????
disp('Sampling Number in slow-time domain');disp(N)     %方位向取樣點(diǎn)數(shù) 
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);      %全零數(shù)組   N為方位向采樣點(diǎn)，M為距離向采樣點(diǎn)
Pa = v*sn;              %sample position to target in slow-time domain  定位在慢時間段的樣點(diǎn)位置
R = sqrt(R0^2+(Pa).^2); %sample slant range   抽樣點(diǎn)（即約為發(fā)射脈沖的位置）距離點(diǎn)目標(biāo)的距離
delay = 2*R/C;          %time delay in transmitted pulse   發(fā)射脈沖的時間延遲
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));      %經(jīng)過檢波之后的連續(xù)回波信號抽樣化后的矩陣序列的指數(shù)部分 
Srnm = Ptarget(1,3)*exp(phase).*(abs(Pr)<=Tr/2).*(abs(Pa'*ones(1,M))<=Lsar/2); 
%sig : target echo model 點(diǎn)目標(biāo)回波原始數(shù)據(jù)（將原始連續(xù)回波信號抽樣化后的數(shù)據(jù)，這樣利于作基2-快速傅立葉變換）
%--------------------------------------------------------%
%--------------------------------------------------------%
% Range compression          %距離壓縮
fr = -1/2:1/M:(1/2-1/M);     %M個抽樣點(diǎn)將其化為M-1等分  
fr = fr*Fsr;                 %此處抽樣頻率約等于信號帶寬（前面有結(jié)果=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'); %距離向先進(jìn)行FFT，后進(jìn)行IFFT    
row = C*tm/2;                   %行？？？？？？？？？？？？？？？？？？？/？？？？？？？？？？？？？？？？？？？
col = v*sn;                     %列
colormap(gray);                 %colormap將當(dāng)前顏色表置為Matlab內(nèi)置的灰度圖
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個抽樣點(diǎn)將其化為N-1等分
fa = fa*PRF;                    %為什么要用 PRF而不用Bd ？？？？？？？？？？？？？？？？？？？？？？
deteR = 2/C*(lambda^2*R0*fa.^2/8/v^2)*Fsr;     % 距離徙動所占的采樣點(diǎn)數(shù),不一定為整數(shù)值,所以接下來要進(jìn)行插值
N_interp = 8;                   %N_interp=8表示以8倍進(jìn)行插值
N_ex = 1000;                    %N_ex必須大于距離徙動所占的采樣點(diǎn)數(shù)
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');%經(jīng)過插值之后
%--------------------------------------------------------%
%--------------------------------------------------------%
% Azimuth compression  方位壓縮
fa = -1/2:1/N:(1/2-1/N);
fa = fa'*PRF;      %注意要轉(zhuǎn)置，為了下面第（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');