clear;
%仿真合成孔徑雷達信號
colormap(gray(256))
cj=sqrt(-1);
pi2=2*pi;

c=3e8;              % 傳播速度
fc=200e6;           % 頻率
lambda=c/fc;        % 波長
k=pi2/lambda;       % 波數
Xc=1.e3;            % 距離目標中心的距離

L=400;              % 合成孔徑大小為2*L
Y0=100;             % 目標區域位于[Yc-Y0,Yc+Y0]
Yc=0;               % 距離目標中心的跨距離

theta_c=atan(Yc/Xc);  %到目標中心的斜視角
Rc=sqrt(Xc^2+Yc^2);   %到目標中心的斜視距離
kus=2*k*sin(theta_c); %Doppler頻率偏移

Xcc=Xc/(cos(theta_c)^2);       
DY=(Xcc*lambda)/(4*L);         % 交叉距離分辨率
L_min=max(Y0,L);               % 補零的孔徑

% 壓縮信號的u域參數和陣列
du=(Xcc*lambda)/(4*(Y0+L));    % 孔徑域的采樣間隔
duc=(Xcc*lambda)/(4*Y0);       %壓縮信號孔徑域的采樣間隔
mc=2*ceil(L_min/duc);            
uc=duc*(-mc/2:mc/2-1);        
dkuc=pi2/(mc*duc);                % ku域的采樣間隔
kuc=dkuc*(-mc/2:mc/2-1);          % kuc陣列
%
dku=dkuc;                         % ku域的采樣間隔

%合成孔徑信號的u域參數和陣列
m=2*ceil(pi/(du*dku));            %樣本數
du=pi2/(m*dku);
u=du*(-m/2:m/2-1);                % 合成孔徑陣列
ku=dku*(-m/2:m/2-1);              % ku陣列
%
ntarget=2;             % 目標數
%目標的坐標和反射率
yn(1)=.7*Y0;         fn(1)=0.5;
yn(2)=yn(1)-4*DY;    fn(2)=1;

%測量的回波信號
s=zeros(1,mc);         
for i=1:ntarget;
 dis=sqrt(Xc^2+(Yc+yn(i)-uc).^2);
 s=s+fn(i)*exp(-cj*2*k*dis).*(abs(uc) <= L);
end;
figure(1);
plot(uc,real(s));
xlabel('合成孔徑 u');
ylabel('幅值 A');
title('孔徑信號的實部');
axis([uc(1) uc(mc) 1.1*min(real(s)) 1.1*max(real(s))]);

sig=s(1:50);
%計算Wigner-Ville分布
[tfr,t,f]=tfrwv(sig',1:length(sig),256);
%顯示Wigner-Ville分布等高線
figure(2);
mesh(t,f,tfr);
xlabel('時間 t');
ylabel('頻率 f');
zlabel('幅值 A');

%Hough變換
[WH,rho,theta]=htl(tfr);
figure(3);
mesh(rho,theta,WH');
xlabel('極半徑 rho');
ylabel('角度 theta');
zlabel('幅值 A');