wk.m

基于matlab的几个kalman滤波以及多传感器融合和wk算法的程序

%w-k SAR 成像算法，正侧视，掠射角任意。成像区域的方位向宽度可调。无拼贴过程。
clear
clc
time0=cputime;
%线形调频信号参数
fc=500e+6;
B=250e+6;
tp=3e-7;
K=B/tp;
fs=3.5*B;
dt=1/fs;
N=ceil(tp/dt);

%观察区域定义
c=3e+8;%m/s
dx=5;%m
dz=c/2*dt;

depression=0*pi/180;
Hr0=200;                               %the center horizontal range
Hrn=Hr0-10;                            %the nearest horizontal range
Hrf=Hr0+10;                            %the farest horizontal range
Sr0=Hr0/cos(depression);               %the center slant range
Srn=Hrn/cos(depression);               %the nearest slant range
Srf=Hrf/cos(depression);               %the farest slant range

ka=1;                                  %weight factor
pa=5;                                %azimuth resolution
thetamin=asin(ka*c/(4*pa*(fc+B/2)));   %the minimum semi beamangle.about 8.6 degree
thetamed=asin(ka*c/(4*pa*fc));         %the medium semi beamangle.about 12.4 degree
thetamax=asin(ka*c/(4*pa*(fc-B/2)));   %the maximum semi beamangle.about 22.0 degree
theta=thetamed;                        %处理角 ---min or med or max
halfaperture=Srf*tan(theta);           %天线波束方位向覆盖区域半宽度
halfaperturewidth=fix(halfaperture/dx);%
imghalfaperture=10;%（米）成像区域半宽度
imghalfaperturewidth=fix(imghalfaperture/dx)+1;%采样点个数
delaymax=2*sqrt(Srf^2+(dx*imghalfaperturewidth)^2)/c+tp/2;%距离波门结束时刻； 
delaymin=2*Srn/c-tp/2;%距离波门起始时刻；
z0=Sr0;
shift=round(2*(z0-Srn)/(c*dt));
slicelen=ceil((delaymax-delaymin)/dt);%the number of samples within a range slice
rawlen=slicelen+N-1;%the length of output after match-filtering
rpclen=round(2*(sqrt(Srf^2+(dx*imghalfaperturewidth)^2)-Srn)/c/dt);%the length of processing
%rpclen=rawlen;

%Range reference signal
t=[0:N-1]*dt-tp/2;
s0=exp(j*2*pi*(fc+0.5*K*t).*t);
xf=zeros(1,rawlen);
xf(1:N)=s0;
xf=conj(fft(xf));

%Place Num targets:tar(target number,(azimuth,range,re(amplitude),im(amplitud),gamma,alpha) 
%gammay与alpha表示目标散射的方位向特性，即目标回波被限制在%(gamma-alpha/2,gamma+alpha/2)范围内。
%Two way:the random field as rnf==1 and the special position as rnf==0
Num=7;
tar=zeros(Num,4);
tar(1,1)=0.0;tar(1,2)=Sr0+2.0;tar(1,3)=1.0;tar(1,4)=.0;
tar(2,1)=0.0;tar(2,2)=Sr0+0.0;tar(2,3)=1.0;tar(2,4)=.0;
tar(3,1)=0.0;tar(3,2)=Sr0-2.0;tar(3,3)=1.0;tar(3,4)=.0;
tar(4,1)=4.0;tar(4,2)=Sr0+1.0;tar(4,3)=1.0;tar(4,4)=.0;
tar(5,1)=4.0;tar(5,2)=Sr0-1.0;tar(5,3)=1.0;tar(5,4)=.0;
tar(6,1)=-4.0;tar(6,2)=Sr0+1.0;tar(6,3)=1.0;tar(6,4)=.0;
tar(7,1)=-4.0;tar(7,2)=Sr0-1.0;tar(7,3)=1.0;tar(7,4)=.0;
if (min(tar(:,1))<-imghalfaperturewidth*dx)|(max(tar(:,1))>imghalfaperturewidth*dx)
   error('ERROR AZIMUTH!!!!!!!!!!!!');
end
if (min(tar(:,2))<Srn)|(max(tar(:,2))>Srf)
   error('ERROR RANGE!!!!!!!!!!!!');
end

%Forming echo,i.e rawdata
rpchalfaperturewidth=halfaperturewidth+imghalfaperturewidth;%雷达方位向积累半宽度
rawdatawidth=2*rpchalfaperturewidth+1;
rpcdata=zeros(rawdatawidth,rpclen);%the data after range compression
for i=1:rawdatawidth
  temp=zeros(1,rawlen);
  for m=1:Num
    dazimuth=(i-1-fix(rawdatawidth/2))*dx-tar(m,1);
    if(abs(dazimuth)>tar(m,2)*tan(theta))%雷达在位置i时，目标m在波束照射范围外；
      temp=temp;
    else
      delaym=2*sqrt(dazimuth^2+tar(m,2)^2)/c-tp/2-delaymin;%雷达在位置i时，目标m的延时；
      t0=fix(delaym/dt)+1;t1=fix((delaym+tp)/dt);tm=[t0*dt-delaym:dt:t1*dt-delaym]-tp/2;
      temp(t0:t1)=temp(t0:t1)+exp(j*2*pi*(fc+0.5*K*tm).*tm);
    end
  end
  temp=ifft(xf.*fft(temp));
  rpcdata(i,:)=temp(1:rpclen);       %Range compression
end%for i=1:rawdatawidth
rpcdata=fft(rpcdata,[],2).*exp(j*2*pi*(ones(rawdatawidth,1)*[0:rpclen-1])*shift/rpclen);
rpcdata=ifft(rpcdata,[],2);
time1=cputime;
scenetime=time1-time0 

%Azimuth compression(w-k algorithm)
time0=cputime;
w=zeros(1,rpclen);
kz=zeros(1,rpclen);
kx=zeros(rawdatawidth,1);
dw=2*pi*fs/rpclen;
dkx=2*pi/(dx*rawdatawidth)
w=[0:rpclen-1]*dw;
kz=2/c*w;
kx(1:rpchalfaperturewidth+1)=[0:rpchalfaperturewidth]*dkx;
kx(rpchalfaperturewidth+2:rawdatawidth)=[-rpchalfaperturewidth:-1]*dkx;
rpcdata=fft2(rpcdata);
img=zeros(rawdatawidth,rpclen);
img(1,:)=rpcdata(1,:).*exp(-j*(z0-Srn)*kz);
for i=2:rawdatawidth,
   wcz=c/2*sqrt(kz.^2+kx(i)^2);%1*rpclen vector
   weight=sinc(ones(rpclen,1)*wcz/dw-[0:rpclen-1]'*ones(1,rpclen));%rpclen*rpclen matrix
   temp=rpcdata(i,:)*weight;%1*rpclen vector
   img(i,:)=temp.*exp(-j*(z0*2/c*wcz-Srn*kz))*c.*abs(kz)./(2*wcz);
end
img=ifft2(img);%img=img/max(max(abs(img)));
imgrangesamples=fix((Srf-Srn)/(c/2*dt))+1;
imgazimuthsamples=2*imghalfaperturewidth+1;
img=conj((img(halfaperturewidth+1:halfaperturewidth+imgazimuthsamples,1:imgrangesamples))');
time1=cputime;
imagingtime=time1-time0
kx0=round(imgazimuthsamples/2);
for i=1:imgrangesamples
   img(i,:)=img(i,:).*exp(j*2*pi*kx0*[0:imgazimuthsamples-1]/imgazimuthsamples);%这时,kx0对应kx=0
end
%figure;plotbw(1,img);grid;xlabel('azimuth');ylabel('range');
figure;
image([-imghalfaperturewidth,imghalfaperturewidth]*dx,Srn+[0:imgrangesamples-1]*dz,...
                                        63*abs(img)/max(max(abs(img))))
axis xy;grid;xlabel('azimuth(m)');ylabel('range(m)');