omegaka.m

波数域成像算法

% OmegaKA.m   仿真点目标SAR成像原理
%-------------------------------------------------
% 产生雷达脉冲，经过单点反射之后接收数据
% 数据存放在一个二维数组中，即距离向和方位向
% 采用精确或近似的Omega-K成像算法
% 先经过二维FFT，将信号变换到二维频域，做RFM，对参考距离上的目标完成精确成像
% 再沿距离向频率轴进行Stolt插值，完成differential压缩，最后做二维IFFT，转换到图像域 
%%================================================================
clear;clc;close all;
%%================================================================
%%Parameter--constant
C=3e8;                            %propagation speed
%%Parameter--radar characteristics
Fc=5.3*1e9;                          %carrier frequency 5.3GHz
lambda=C/Fc;                 %wavelength 
%%Parameter--target area
Xmin=-50;                          %target area in azimuth is within[Xmin,Xmax]
Xmax=50;                                  
Yc=20000;                      %center of imaged area 20km
Y0=2000;                          %target area in slant range is within[Yc-Y0,Yc+Y0]
%%Parameter--orbital information
V=150;                            %SAR velosity 150 m/s
H=6000;                          %height 6 km
R0=Yc;
%%Parameter--antenna
thita=0*pi/180;         %squint angle
D=2.5;                                %antenna length in azimuth direction
R0c=R0/cos(thita);        %the slant range to the target at the time it is illuminated by the beam center
Lsar=lambda*R0c/D/cos(thita);         %SAR integration length
Tsar=Lsar/V;                   %SAR integration time
%%Parameter--slow-time domain
Ka=-2*V^2*cos(thita)^3/lambda/R0;    %doppler frequency modulation rate
Ba=abs(Ka*Tsar);           %doppler frequency modulation bandwidth
PRF=Ba*1.25;                         %pulse repitition frequency
PRT=1/PRF;                   %pulse repitition time
ds=PRT;                         %sample spacing in slow-time domain
Nslow=ceil((Xmax-Xmin+Lsar)/V/ds); %sample number in slow-time domain
Nslow=2^nextpow2(Nslow);              %for fft
% Nslow=20000;
sn=linspace((Xmin-Yc*tan(thita)-Lsar/2)/V,(Xmax-Yc*tan(thita)+Lsar/2)/V,Nslow);%discrete time array in slow-time domain
PRT=(Xmax-Xmin+Lsar)/V/Nslow;    %refresh
PRF=1/PRT;
ds=PRT;
%%Parameter--fast-time domain
Tr=5e-6;                         %pulse duration 5us
Br=100e6;                        %chirp frequency modulation bandwidth 100MHz


Fsr=1.2*Br;                        %sampling frequency in fast-time domain
dt=1/Fsr;                         %sample spacing in fast-time domain
if ((Yc-Y0)*tan(thita)-Lsar/2)>0
    Rmin=sqrt((Yc-Y0)^2+((Yc-Y0)*tan(thita)-Lsar/2)^2);    
else
    Rmin=abs(Yc-Y0);
end
Rmax=sqrt((Yc+Y0)^2+((Yc+Y0)*tan(thita)+Lsar/2)^2);
Nfast=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt);%sample number in fast-time domain
Nfast=2^nextpow2(Nfast);                   %for fft
% Nfast=1300;
tm=linspace(2*Rmin/C-Tr/2,2*Rmax/C+Tr/2,Nfast)-(2*Rmin/C-Tr/2); %discrete time array in fast-time domain
Rmin1=abs(Rmin);
Rmax1=sqrt((Yc+Y0)^2+((Yc+Y0)*tan(thita)+Lsar/2)^2)/2+(Yc-Y0+sqrt((2*Y0).^2+((Yc+Y0)*tan(thita)+Lsar/2)^2))/2;
tm1=linspace(2*Rmin1/C-Tr/2,2*Rmax1/C+Tr/2,Nfast)-(2*Rmin1/C-Tr/2); %discrete time array in fast-time domain

dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast;    %refresh
Fsr=1/dt;
Ntr = floor(Tr*Fsr);	    % total points of pulse
half_Ntr = floor(Tr*Fsr/2);	    % total points of half pulse
%%Parameter--point targets
Ntarget=1;                            %number of targets
% format [x, y, reflectivity]
% Ptarget=[
% %     -3*600*tan(thita),Yc-1800/2,1               %position of targets
%          -2*600*tan(thita),Yc-1200,1
%          -600*tan(thita),Yc-600,1
%          0,Yc,1
%          600*tan(thita),Yc+600,1
%          2*600*tan(thita),Yc+1200,1
% %          3*600*tan(thita),Yc+1800/2,1
%          ] 
Ptarget=[
         0,Yc,1
         ]; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% run start time
disp(strcat('program begins at : ',datestr(now)));
run_start = rem(now,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%产生原始数据
%
Kr=Br/Tr;                        %chirp slope
% for kkk=-2:2
%     if kkk==0
%         kkk=1;
%     else
%     Kr=1/2*kkk*Kr;
jamzip=[0,20000,200];
recvsig = zeros(Nslow,Nfast,'single');
% recvsig1 = zeros(Nslow,Nfast);
% h = waitbar(0,'SAR回波生成');
for k=1:1:Ntarget
    sigma=Ptarget(k,3);
    Dslow=sn*V-Ptarget(k,1);
    R=sqrt(Dslow.^2+Ptarget(k,2)^2);
    R00=sqrt((Dslow-1000).^2+(Ptarget(k,2)-1000)^2);
    tau=(2*R-R00)/C-(2*Rmin/C-Tr/2);
    Dfast=single(ones(Nslow,1)*tm-tau'*ones(1,Nfast));
    phase=single(pi*Kr*Dfast.^2-(2*pi/lambda)*((2*R-R00)'*ones(1,Nfast)));
    recvsig=recvsig+sigma*exp(j*phase);
%    recvsig=recvsig+sigma*exp(j*phase).*(-Tr/2<Dfast&Dfast<Tr/2).*((abs(Dslow+Ptarget(k,2)*tan(thita))<Lsar/2)'*ones(1,Nfast));
    %     waitbar(k/Ntarget)
%     Dslow1=sn*V-Ptarget(k,1);
%     R1=sqrt(Dslow.^2+Ptarget(k,2)^2)/2+sqrt((sn*V).^2+jamzip(2)^2)/2+sqrt((Ptarget(k,2)-jamzip(2)).^2+Ptarget(k,1).^2)/2;
%     tau1=2*R1/C-(2*Rmin/C-Tr/2);    Dfast1=ones(Nslow,1)*tm1-tau1'*ones(1,Nfast);
%     phase1=pi*Kr*Dfast1.^2-(4*pi/lambda)*(R1'*ones(1,Nfast));
%     recvsig1=recvsig1+sigma*exp(j*phase1).*(-Tr/2<Dfast1&Dfast1<Tr/2).*((abs(Dslow+Ptarget(k,2)*tan(thita))<Lsar/2)'*ones(1,Nfast));
end
% for k=1:1:Ntarget
%     sigma=Ptarget(k,3);
%     Dslow=sn*V-Ptarget(k,1);
%     R=sqrt(Dslow.^2+(Ptarget(k,2)+300)^2);
%     tau=2*R/C-(2*Rmin/C-Tr/2);
%     Dfast=ones(Nslow,1)*tm-tau'*ones(1,Nfast);
%     phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,Nfast));
%     recvsig1=recvsig1+sigma*exp(j*phase).*(-Tr/2<Dfast&Dfast<Tr/2).*((abs(Dslow+Ptarget(k,2)*tan(thita))<Lsar/2)'*ones(1,Nfast));
%     %     waitbar(k/Ntarget)
% %     Dslow1=sn*V-Ptarget(k,1);
% %     R1=sqrt(Dslow.^2+Ptarget(k,2)^2)/2+sqrt((sn*V).^2+jamzip(2)^2)/2+sqrt((Ptarget(k,2)-jamzip(2)).^2+Ptarget(k,1).^2)/2;
% %     tau1=2*R1/C-(2*Rmin/C-Tr/2);    Dfast1=ones(Nslow,1)*tm1-tau1'*ones(1,Nfast);
% %     phase1=pi*Kr*Dfast1.^2-(4*pi/lambda)*(R1'*ones(1,Nfast));
% %     recvsig1=recvsig1+sigma*exp(j*phase1).*(-Tr/2<Dfast1&Dfast1<Tr/2).*((abs(Dslow+Ptarget(k,2)*tan(thita))<Lsar/2)'*ones(1,Nfast));
% end
%  recvsig1(1:Nslow-1,:)=recvsig1(2:Nslow,:);

% recvsig=2*recvsig1+recvsig;
clear recvsig1 phase1 ;
% close(h);
% 输出原始数据二维图形
clear phase Dfast Dfast1;
figure(1);
% imagesc(real(recvsig(100:400,1730:2360)));
imagesc(real(recvsig(100:400,1300:2800)));
colormap(gray);
axis xy;
% axis auto;
% axis normal;
%  axis xy;
title('原始回波'),xlabel('range direction'),ylabel('azimuth direction');
recvsig1=recvsig;
clear recvsig;

% title('paired-echo signal dispersed by dot target'),xlabel('range direction'),ylabel('azimuth direction');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%对方位向和距离向进行FFT 将信号转换到二维频域
%  for kkk=-1:2
%      if kkk==0
%        kkk=1
%      else
% Kr=(kkk/2)*Kr;
Kr=Br/Tr;  
recvsig = zeros(Nslow,Nfast);
recvsig=recvsig1;
for nr = 1:Nfast
	recvsig(:,nr) =fftshift(fft(recvsig(:,nr)));  
end
for na = 1:Nslow
	recvsig(na,:) = fftshift(fft(recvsig(na,:)));   
end
% figure(2);
% subplot(2,1,1),imshow(abs(recvsig(:,:)),[ ]); %,axis xy;
% title('amplitude signal in two-dimensional frequency domain '),xlabel('range freqency'),ylabel('Doppler freqency');
% subplot(2,1,2),imshow(angle(recvsig(:,:)),[ ]); %,axis xy;
% title('phase signal in two-dimensional frequency domain '),xlabel('range freqency'),ylabel('Doppler freqency');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%做RFM
fyitac=2*V*sin(thita)/lambda; %reference azimuth frequency
ftao=linspace(-Fsr/2,Fsr/2,Nfast);%kx域序列
if (((fyitac+Ba/2)-floor(fyitac/PRF)*PRF)<PRF/2)&(((fyitac-Ba/2)-floor(fyitac/PRF)*PRF)>-PRF/2)
    fyita=linspace(-PRF/2,PRF/2,Nslow)+floor(fyitac/PRF)*PRF;
    Nm=Nslow;
else
    Nm=floor((Nslow-Ba/PRF*Nslow)/2+((fyitac+Ba/2)-floor(fyitac/PRF)*PRF-PRF/2)*Nslow/PRF);    
    fyita=linspace(-PRF/2,PRF/2,Nslow)+floor(fyitac/PRF)*PRF;
    fyita(1:Nm)=fyita(1:Nm)+PRF;
end
for na = 1:Nslow
    recvsig(na,:)=recvsig(na,:).*exp(j*4*pi*R0*sqrt(((Fc+ftao)/C).^2-(fyita(1,na)/V)^2/4)+j*pi/Kr*ftao.^2-j*2*pi*ftao*(2*Rmin/C-Tr/2));  
 end
% figure(3);
% subplot(2,1,1),imshow(abs(recvsig(:,:)),[ ]); %,axis xy;
% title('amplitude signal in two-dimensional frequency domain after RFM'),xlabel('range freqency'),ylabel('Doppler freqency');
% subplot(2,1,2),imshow(angle(recvsig(:,:)),[ ]); %,axis xy;
% title('phase signal in two-dimensional frequency domain after RFM'),xlabel('range freqency'),ylabel('Doppler freqency');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%做stolt插值
kx=ftao;%kx域序列
ky=(fyita)*C/2/V;%ky域序列
KX=ones(Nslow,1)*kx;%扩充为矩阵
KY=ky'*ones(1,Nfast);%扩充为矩阵
DKX=abs(sqrt((KX+Fc).^2+KY.^2)-(KX+Fc));%待插值点的弯曲量
P=8/2;%8点sinc插值
NDKX=fix(DKX*Nfast/Fsr);
Nd=NDKX(Nm,1);
B=[recvsig recvsig(:,1:Nd+P)];
clear KX KY;
recvsig=zeros(Nslow,Nfast);%初始化
% h=waitbar(0,'Stolt插值');
for na=1:Nslow
    for nr=P:Nfast
       NN=NDKX(na,nr)+(-P+1:P);
       recvsig(na,nr)=sinc(DKX(na,nr)*Nfast/Fsr-NN)*B(na,nr+NN)';%利用采样定理恢复待插值点的数值
    end
%     waitbar(na/Nslow) 
end
% close(h);
% figure(4);
% subplot(2,1,1),imshow(abs(recvsig(:,:)),[ ]); %,axis xy;
% title('amplitude signal in two-dimensional frequency domain after Stolt interpolation'),xlabel('range freqency'),ylabel('Doppler freqency');
% subplot(2,1,2),imshow(angle(recvsig(:,:)),[ ]); %,axis xy;
% title('phase signal in two-dimensional frequency domain after Stolt interpolation'),xlabel('range freqency'),ylabel('Doppler freqency');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%对距离向进行IFFT 将信号转换到RD域
for na = 1:Nslow
    recvsig(na,:) = fftshift(ifft(fftshift(recvsig(na,:))));
end
% figure(5);
% subplot(2,1,1),imshow(abs(recvsig(:,:)),[ ]); %,axis xy;
% title('signal after range compression'),xlabel('range direction'),ylabel('Doppler direction');
% subplot(2,1,2),imshow(angle(recvsig(:,:)),[ ]); %,axis xy;
% title('signal after range compression'),xlabel('range direction'),ylabel('Doppler direction');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%对方位向进行IFFT 将信号转换到二维图像域
% Dfyita=sqrt(1-(C*fyita/2/V/Fc).^2);
for nr = 1:Nfast
%     recvsig(:,nr)=recvsig(:,nr).*exp(-j*2*pi*(nr-Nfast/2)*dt*Fc*Dfyita)';
	recvsig(:,nr) = ifft(fftshift(recvsig(:,nr)));  
end
%输出二维图形
G=20*log10(abs(recvsig)+1);
gm=max(max(G));
recvsig=255/gm*G;
figure;imagesc(255-abs(recvsig(100:400,1300:2800)));
% imshow(abs(recvsig(:,:)),[]); %,axis xy;

colormap(gray);
axis xy;
title('点目标成像'),xlabel('range direction'),ylabel('Doppler direction');
figure;mesh(recvsig(:,:)); figure(gcf)
%      end
%  end
% figure(9);
% mesh(Ro,((1:Nfast)-386)*PRT*V,abs(recvsig(:,:)));
% title('signal after range & azimuth compression'),xlabel('range direction (m)'),ylabel('azimuth direction (m)');

% mymax = max(max(abs(recvsig(:,:))));  
% 
% figure(8);
% contour(Ro,((1:Nfast)-386)*PRT*V,abs(recvsig(:,:)),12);
% hold on;
% contour(sqrt((range_t*C/2).^2-H.^2),V*azim_t,abs(recvsig(:,:)),[mymax,mymax]);
% hold off;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% run end time
% disp(strcat('program finishs at : ',datestr(now)));
% run_end = rem(now,1);
% % display total run time
% fprintf('\nTotal run time is %f seconds.\n',(run_end-run_start)*3600*24);