SARAutoFocus.rar_PGA_autofocus matlab_radar_成像_自聚焦成像

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ¡¡ STRIPMAP SAR¡¡AUTOFOCUS SIMULATION % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear %%%%%%%%%%%%%%% constant parameters %%%%%%%%%%%%%%%%%%% c=3e8; % propagation speed,3*10^8m/s %%%%%%%%%%%%%%% chirp signal parameters %%%%%%%%%%%%%%% fc=1.5e9; % carrier frequency,1.5GHz wc=2*pi*fc; lambda_c=c/fc; % Wavelength at carrier frequency Tr=1.5e-6; % chirp pusle duration 1.5us Br=150e6; % chirp frequebcy bandwidth 150Mhz Kr=Br/Tr; % chirp rate %%%%%%%%%%%%%% observating strip parameters %%%%%%%%%%%% Rc=3000; % Range distance to center of target area,3000m R0=150; % Target area in range is within [Rc-R0,Rc+R0] X0=200; % Target area in cross-range is within [-X0,X0] %%%%%%%%%%%% range ,r/t domain parameters %%%%%%%%%%%%%%% Nr=512; % The number of r (fast-time) domain samples r=Rc+linspace(-R0,R0,Nr); t=2*r/c; % t domain series (fast-time domain) dt=4*R0/c/Nr; % sample spacing in fast-time domain (t domain) f=linspace(-1/2/dt,1/2/dt,Nr);% f domain (FT about t) %%%%%%%%%%%Azimuth, cross range parameters x/u domian %%%%%% v=100; % speed of the SAR aircraft 100m/s theta=0; % squint angle of beam center Lsar=300; % synthesize aperture length 300m Tsar=Lsar/v; % synthesize aperture time Na=1024; % The number of x (slow-time) domain samples x=linspace(-X0,X0,Na); % x domian series (cross-range) u=x/v; % u (slow-time) domain du=2*X0/v/Na; % sample spacing in slow-time domain (u domain) fu=linspace(-1/2/du,1/2/du,Na); % fu series(FT about u domain) fdc=2*v*sin(theta)/lambda_c; % the centric Doppler frequency fr=-2*v^2/(lambda_c*Rc); % Doppler chirp rate Ba=abs(fr)*Tsar; % Doppler frequency bandwidth Ba %%%%%%%%%%%%% SAR¡¡Resolution %%%%%%%%%%%%%%%%%%%%%%%%% Dr=c/2/Br; %range resolution DX=v/Ba; %cross-range resolution %%%%%%%%%%%%% target parameters ( one point target) %%%%%%%%%%%%%%%%%%%%%%% % yn: range; xn= cross-range; fn: reflectivity rn=Rc; xn=0; fn=1 ; %%%%%%%%%%%%%% targets echo %%%%%%%%%%%%%%%%%% s_tu=zeros(Na,Nr); U=u'*ones(1,Nr); % extend to matrix T=ones(Na,1)*t; R=zeros(Na,Nr); DT=zeros(Na,Nr); R=sqrt(rn^2+(xn-v*U).^2); DT=T-2*R/c; phase=pi*Kr*DT.^2-4*pi*R/lambda_c; s_tu=fn*exp(j*phase).*(abs(DT)<Tr/2).*((abs(xn-v*U))<Lsar/2); %%%%%%%%%%%%%%% Range Compression %%%%%%%%%%%%%%% p0_t=exp(j*pi*Kr*(t-2*Rc/c).^2).*(abs(t-2*Rc/c)<Tr/2); % Range referrence signal p0_f=fftshift(fft(fftshift(p0_t))); s_fu=fftshift(fft(fftshift(s_tu.'))).'; % Range FFT src_fu=s_fu.*(ones(Na,1)*conj(p0_f)); % Range Compression src_tu=fftshift(ifft(fftshift(src_fu.'))).'; % signal after Range Compression src_tfu=fftshift(fft(fftshift(src_tu))); % Azimuth FFT and Range-Doppler domain %%%%%%%%%%%%%%% RCM Compensation %%%%%%%%%%%%%%% src_ffu=fftshift(fft(fftshift(src_fu))); % 2D spectrum after RC F=ones(Na,1)*f; FU=fu'*ones(1,Nr); p0_2f=exp(j*pi/fc^2/fr*(FU.*F).^2+j*pi*fdc^2/fc/fr*F-j*pi/fc/fr*FU.^2.*F); % RCM Compensation s2rc_ffu=(src_ffu).*p0_2f; s2rc_tfu=fftshift(ifft(fftshift(s2rc_ffu).')).'; % Range-Doppler domain s2rc_tu=fftshift(ifft(fftshift(s2rc_tfu))); % t-u domain %%%%%%%%%%%%%% Azimuth Compression %%%%%%%%%%%%%%% % p0_fu=exp(j*pi/fr*(FU-fdc).^2); % Azimuth referrence signal in Doppler domain % s2rcac_tfu=s2rc_tfu.*p0_fu; % Azimuth Compression % s2rcac_tu=fftshift(ifft(fftshift(s2rcac_tfu))); % Amizuth IFFT % %%%%%%%%%%%%%% AutoFocus by PGA algorithm %%%%%%%%%%%%%%%%%%%% % threshold=1e-3*pi; % Nw=127; % the length of window Nw % window=zeros(Na,1); % the triangle window % TriWindow=triang(Nw); % window((Na/2-(Nw-1)/2):Na/2)=TriWindow(1:(Nw+1)/2); % window(Na/2+1:Na/2+(Nw+1)/2)=TriWindow((Nw+1)/2:end); % dfr=pi; % the variety of Doppler chirp rate % count=0; % while abs(dfr)>threshold % count=count+1 % DistinctPoint=max(max(abs(s2rcac_tu))); % find the DistinctPoint % [index_a index_r]=find(abs(s2rcac_tu)==DistinctPoint); % the range unit index_r % gn=s2rcac_tu(:,index_r); % gn=circshift(gn,Na/2-index_a); % shift the DistinctPoint to the center % gn=gn.*window; % n=[-Na/2:-1 1:Na/2]; % nk=[Na/2-(Nw-1)/2:Na/2+(Nw-1)/2+1]; % nkk=-(Nw-1)/2:(Nw-1)/2+1; % gdn=-j*n'.*gn; % the differential of signal after RC,RCM % Gdk=fftshift(fft(gdn)); % Gk=fftshift(fft(gn)); % dphase=imag(Gdk(nk).*conj(Gk(nk)))./(Gk(nk).*conj(Gk(nk))); % phase gradient % [a S]=polyfit(nkk,dphase',1); % estimate the slope rate a(1) by LS % dfr=a(1)*du*(fr)^2; % the variance of Doppler chirp rate % %%%%%%%%% Azimuth Compression %%%%%%%%%% % fr=fr+dfr; % p0_fu=exp(j*pi/fr*(FU-fdc).^2); % Azimuth referrence signal in Doppler domain % s2rcac_tfu=s2rc_tfu.*p0_fu; % Azimuth Compression % s2rcac_tu=fftshift(ifft(fftshift(s2rcac_tfu))); % Amizuth IFFT % end %%%%%%%%%% Minimun Entropy Algorithm %%%%%%%%%%%%%% frc=-2*v^2/(lambda_c*Rc); % intial fr0=abs(-2*v^2/(lambda_c*Rc))*0.2; % search region [frc-fr0 frc+fr0] Ns=128; % number of sampling frs=linspace(frc-fr0,frc+fr0,Ns); H=zeros(1,Ns); for i=1:Ns %%%%%%%%%%%%%%% Azimuth Compression %%%%%%%%%%%%%%% p0_fu=exp(j*pi/frs(i)*(FU-fdc).^2); % Azimuth referrence signal in Doppler domain s2rcac_tfu=s2rc_tfu.*p0_fu; % Azimuth Compression s2rcac_tu=fftshift(ifft(fftshift(s2rcac_tfu))); % Amizuth IFFT s2rcac_X=abs(s2rcac_tu); s2rcac_X=s2rcac_X./(ones(Na,1)*sum(s2rcac_X)); H(i)=-sum(sum(s2rcac_X.*log(s2rcac_X))); end [Hmin index]=min(H); p0_fu=exp(j*pi/frs(index)*(FU-fdc).^2); % the minumum entropy Doppler chirp rate s2rcac_tfu=s2rc_tfu.*p0_fu; % Azimuth Compression s2rcac_tu=fftshift(ifft(fftshift(s2rcac_tfu))); %%%%%%%%%%%%%%%% Image show %%%%%%%%%%%%%%%%% % Show Measured Stripmap SAR Signal figure(1) G=20*log10(abs(s_tu)+1e-6); xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng); colormap(gray) image(r-Rc,x,256-cg*(G-ng)); axis([-R0 R0 -X0 X0]) grid on,axis tight xlabel('Range') ylabel('Azimuth') title('Measured Stripmap SAR Signal') % Image after range compression figure(2) G=20*log10(abs(src_tu)+1e-6); colormap(gray) xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng); image(r-Rc,x,256-cg*(G-ng)); axis([-R0 R0 -X0 X0]) grid on,axis tight xlabel('Range') ylabel('Azimuth') title('Image after RC ') % Image after RCM figure(3) G=20*log10(abs(s2rc_tu)+1e-6); colormap(gray) xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng); image(r-Rc,fu,256-cg*(G-ng)); axis([-R0 R0 -X0 X0]) grid on,axis tight xlabel('Range') ylabel('Doppler frequency') title(' Image after RCM') % show image after RC,AC figure(4) G=20*log10(abs(s2rcac_tu)+1e-6); colormap(gray) xg=max(max(G)); ng=xg-60; cg=255/(xg-ng); imagesc(r-Rc,x,cg*(G-ng).*(G>ng)); axis([-R0 R0 -X0 X0]) grid on,axis tight x