www.pudn.com > huibo.rar > huibo.m, change:2012-06-21,size:3333b


tic 
clear all; 
c=3e8;                                 
fc=75e9;    
lamda=c/fc;  
 
% X0=200;                                 
Rc=865000; 
theta=23/180*pi; 
Rr=Rc*sin(theta); 
Z=Rc*cos(theta); 
  
 
% Lsar=lamda*Rc/6; 
beta=66*pi/180; 
% vr=0; 
% vr=0; 
% va=0; 
% vh=0; 
% R0=150;                                
  
 
B=500; 
alpha=-6*pi/180; 
Z1=Z; 
Y1=0; 
Z2=Z+B*sin(alpha); 
Y2=B*cos(alpha); 
 
% R0=Rc*cos(theta); 
% B1=B*cos(theta-alpha); 
% d=B1/R0/sin(theta); 
%  
% fen=lamda/2/d; 
  
 
Tr=4e-6;                           
Br=50e6;                             
Kr=Br/Tr;                             
Nr=512; 
fs=100e6; 
dt=1/fs;                             
  
 
v=7098.2;                                    
Lsar=lamda*Rc/6;                             
Na=512; 
fp=2*v/6*2;                                            %4.5e3; 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     载入地形数据       %%%%%%%%%%%%%%%%%%%%%%%% 
load zhui.mat; 
% load duo.mat; 
% load ping.mat; 
 
% mesh(plane); 
% figure 
% mesh(plane); 
 
tmp=size(find(plane)); 
Axis=zeros(3,tmp(1));  
flagrow=1; 
for flagp=1:size(plane(:,1))                                                   %%% flagp 是plane的行数 
    tmp=max(abs(plane(flagp,:)));                                              %%% tmp是取plane中每一行的最大值 为了寻找不等于0的 
    if tmp 
        tmpp=find(plane(flagp,:));                                             %%% tmpp是记录 plane中 在flagp行的不等于0的数的位置 
        tmpa=size(tmpp);                                                       %%% tmpa(2)是记录 plane中 在flagp行的不等于0的数的个数 
        tmpa=tmpa(2); 
        flagrow=flagrow+tmpa; 
        Axis(1,flagrow-tmpa:flagrow-1)=(flagp-30)*3; 
        Axis(2,flagrow-tmpa:flagrow-1)=(tmpp-20)*3; 
        Axis(3,flagrow-tmpa:flagrow-1)=plane(flagp,tmpp)*1.015;%*1.00487443452357*4/5; 
    end 
end 
  
y=Axis(2,:)*cos(beta); 
z=-Axis(2,:)*sin(beta); 
y=y+Axis(3,:)*sin(beta); 
z=z+Axis(3,:)*cos(beta); 
Axis(2,:)=y; 
Axis(3,:)=z;                        
  
 
echo=zeros(512,512); 
ttmp=(-256:255)*dt; 
tr=ttmp+2*Rc/c;               
% p=2*Rc/c; 
x=(-256:255)/fp*v; 
for flaga=1:Na 
    for flagdot=1:flagrow-1 
        rd=sqrt((Rr-Y1+Axis(2,flagdot))^2+(x(flaga)-Axis(1,flagdot))^2+(Z1-Axis(3,flagdot))^2); 
        td=2*rd/c; 
        wr=rectpuls(tr-td,Tr);         
        wa=rectpuls(x(flaga)-Axis(1,flagdot),Lsar); 
        echo(flaga,:)=echo(flaga,:)+exp(i*2*pi*(-fc*td+1/2*Kr*(tr-td).^2)).*wr*wa; 
    end 
end 
figure  
imagesc(abs(echo)); 
save planechoh1.mat echo; 
disp 1 
%%%%%%%%%%%%%%%%     第二组回波   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
echo=zeros(512,512); 
ttmp=(-256:255)*dt; 
tr=ttmp+(sqrt((Rr-Y2)^2+Z2^2)+Rc)/c;                  
p=(sqrt((Rr-Y2)^2+Z2^2)+Rc)/c; 
x=(-256:255)/fp*v; 
for flaga=1:Na 
    for flagdot=1:flagrow-1 
        rd2=sqrt((Rr-Y2+Axis(2,flagdot))^2+(x(flaga)-Axis(1,flagdot))^2+(Z2-Axis(3,flagdot))^2); 
%         rd1=sqrt((Rr-Y1+Axis(2,flagdot))^2+(x(flaga)-Axis(1,flagdot))^2+(Z1-Axis(3,flagdot))^2); 
        td=2*rd2/c; 
        wr=rectpuls(tr-td,Tr);          
        wa=rectpuls(x(flaga)-Axis(1,flagdot),Lsar); 
        echo(flaga,:)=echo(flaga,:)+exp(i*2*pi*(-fc*td+1/2*Kr*(tr-td).^2)).*wr*wa; 
    end 
end 
figure  
imagesc(abs(echo)); 
save planechoh2.mat echo; 
disp 2 
toc