SAR多普勒中心估计Matlab代码

function fd_center=DopplerCenterEstimation_MID_jie(x,v,PRF,lambda,Rs,DeltaR) % clear; % clc % load yhf; % x=x(2001:6000,:); % v=v_part; % DeltaR=C/2/Fs; [nrn,nan]=size(x); for n=1:nrn c=x(n,2:nan).*conj(x(n,1:nan-1)); fd_center_n(n,1)=angle(mean(c))/(2*pi)*PRF; end % fdc_center=mid(fd_center_n); fdc_center=fd_center_n(1500) rs=[-nrn/2:nrn/2-1]'*DeltaR; % % % % % % ka=-2*v*v*cos(beta)/lambda./(rs+Rs); % % % % % % ka0=-2*v*v*cos(beta)/lambda./Rs; ka=-2*v*v/lambda./(rs+Rs); ka01=-2*v*v/lambda./Rs; t_nan=[-nan/2+1:nan/2]/PRF; f_slow=[-ceil(nan/2):nan-ceil(nan/2)-1]/nan*PRF; f_slow=f_slow+fdc_center; f_slow=fftshift(f_slow); gama=sqrt(1-(lambda*f_slow/2/v).^2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% estimate ka z_slt=zeros(nrn,nan); for m=1:nrn z_slt(m,:)=fft(x(m,:).*exp(-j*2*pi*fdc_center*t_nan)); z_slt(m,:)=ifft(z_slt(m,:).*exp(j*4*pi*(Rs+rs(m))/lambda*gama).*exp(-j*4*pi*(Rs+rs(m))/lambda-j*pi*1/ka(m)*f_slow.^2)).*exp(j*2*pi*fdc_center*t_nan); % z_slt(m,:)=ifft(z_slt(m,:)).*exp(j*2*pi*fdc*t_nan); end for m=1:nrn sref=exp(-j*pi*ka(m)*t_nan.^2); z_slt(m,:)=z_slt(m,:).*sref; end step=4;part_nan=nan/2*step; R_sum=zeros(part_nan,1); for m=1:nrn/2 fy1=abs(fftshift(fft(z_slt(m,[1:nan/2]),nan/2*step))); fy2=abs(fftshift(fft(z_slt(m,[nan/2+1:nan]),nan/2*step))); % fy1=(fy1(nan/2*step/2-part_nan/2+1:nan/2*step/2+part_nan/2)); % fy2=(fy2(nan/2*step/2-part_nan/2+1:nan/2*step/2+part_nan/2)); fy1=sqrt(fy1(nan/2*step/2-part_nan/2+1:nan/2*step/2+part_nan/2)); fy2=sqrt(fy2(nan/2*step/2-part_nan/2+1:nan/2*step/2+part_nan/2)); fy1=orth(fy1.').'; fy2=orth(fy2.').'; [maxR(m,1),mopt(m,1),x1,x2,R]=fastcorr(fy1,fy2); R_sum=R_sum+R; end if 1 R=R_sum; [maxR,maxm]=max(R); M=part_nan; MM=fix(M/2); point=1; xstar=0; if ((maxm~=1)&(maxm~=M*point)) f(1)=R(maxm-1); f(2)=R(maxm); f(3)=R(maxm+1); fa=(f(1)+f(3)-2*f(2))/2; fb=(f(3)-f(1))/2; fc=f(2); xstar=-fb/(2*fa); end mopt=(maxm+xstar-MM*point-1)/point; ka_deviate0=-mopt/(nan/2*step)*PRF/(nan/2/PRF) % ka_estimate=mid(ka_deviate0+ka01); ka_estimate=(ka_deviate0+ka01); ka_estimate = ka_estimate(round(length(ka_estimate)/2)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2 % for m=1:nrn % z_slt(m,:)=fft(x(m,:).*exp(-j*2*pi*fdc_center*t_nan).*exp(-j*pi*ka_deviate0*t_nan.^2)); % z_slt(m,:)=ifft(z_slt(m,:).*exp(j*4*pi*(Rs+rs(m))/lambda*gama).*exp(-j*4*pi*(Rs+rs(m))/lambda-j*pi*1/ka(m)*f_slow.^2)).*exp(j*2*pi*fdc_center*t_nan); % % z_slt(m,:)=ifft(z_slt(m,:)).*exp(j*2*pi*fdc*t_nan); % end % % for m=1:nrn % sref=exp(-j*pi*ka(m)*t_nan.^2); % z_slt(m,:)=z_slt(m,:).*sref; % end % % step=4;part_nan=nan/2*step; % % R_sum=zeros(part_nan,1); % for m=1:nrn/2 % fy1=abs(fftshift(fft(z_slt(m,[1:nan/2]),nan/2*step))); % fy2=abs(fftshift(fft(z_slt(m,[nan/2+1:nan]),nan/2*step))); % % fy1=(fy1(nan/2*step/2-part_nan/2+1:nan/2*step/2+part_nan/2)); % % fy2=(fy2(nan/2*step/2-part_nan/2+1:nan/2*step/2+part_nan/2)); % % fy1=sqrt(fy1(nan/2*step/2-part_nan/2+1:nan/2*step/2+part_nan/2)); % fy2=sqrt(fy2(nan/2*step/2-part_nan/2+1:nan/2*step/2+part_nan/2)); % fy1=orth(fy1.').'; % fy2=orth(fy2.').'; % [maxR(m,1),mopt(m,1),x1,x2,R]=fastcorr(fy1,fy2); % R_sum=R_sum+R; % end % % if 1 % R=R_sum; % [maxR,maxm]=max(R); % % M=part_nan; % MM=fix(M/2); % point=1; % xstar=0; % if ((maxm~=1)&(maxm~=M*point)) % f(1)=R(maxm-1); % f(2)=R(maxm); % f(3)=R(maxm+1); % fa=(f(1)+f(3)-2*f(2))/2; % fb=(f(3)-f(1))/2; % fc=f(2); % xstar=-fb/(2*fa); % end % mopt=(maxm+xstar-MM*point-1)/point; % ka_deviate1=-mopt/(nan/2*step)*PRF/(nan/2/PRF) % ka_estimate=mid(ka01+ka_deviate1+ka_deviate0); % end %%%%%%%%%%%%%%%%%%%%%%% the fdc of first quarter and the last quarter apeture for n=1:nrn c1=x(n,2:nan/4).*conj(x(n,1:nan/4-1)); c2=x(n,2+nan*3/4:nan).*conj(x(n,nan/4*3+1:nan-1)); fd_center_1(n,1)=angle(mean(c1))/(2*pi)*PRF; fd_center_2(n,1)=angle(mean(c2))/(2*pi)*PRF; end fdc_1=mid(fd_center_1); fdc_2=mid(fd_center_2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% estimate the DeltaRR x1=x(:,[1:nan/4]); x2=x(:,[3*nan/4+1:nan]); x1=abs(fftshift(fft(x1))); x2=abs(fftshift(fft(x2))); fy1=orth(sqrt(x1).'); fy2=orth(sqrt(x2).'); fy1=mean(fy1,2); fy2=mean(fy2,2); [maxR,moptt,x1,x2,R]=fastcorr(fy1,fy2); DeltaRR=moptt*DeltaR; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% estimate m De_fdc=fdc_1-fdc_2; m=DeltaRR*2*ka_estimate/lambda/PRF/De_fdc; m=round(m) fd_center=fdc_center-m*PRF; function [maxR,mopt,x1,x2,R,maxm]=fastcorr(x1,x2) [M,N]=size(x1); if (M==1) x1=x1.'; x2=x2.'; M=N; end MM=fix(M/2); %point=4; point=1; Vec_M=[0:M-1]'; a=fft(orth(abs(x1))); b=fft(orth(abs(x2))); R=abs(fftshift(ifft(a.*conj(b),M*point)))*point; [maxR,maxm]=max(R); xstar=0; if ((maxm~=1)&&(maxm~=M*point)) f(1)=R(maxm-1); f(2)=R(maxm); f(3)=R(maxm+1); fa=(f(1)+f(3)-2*f(2))/2; fb=(f(3)-f(1))/2; fc=f(2); xstar=-fb/(2*fa); end mopt=(maxm+xstar-MM*point-1)/point; function [x1]=mid(x) x1 = x(round(length(x)/2));