注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Mr.Right

不顾一切的去想,于是我们有了梦想。脚踏实地的去做,于是梦想成了现实。

 
 
 

日志

 
 
关于我

人生一年又一年,只要每年都有所积累,有所成长,都有那么一次自己认为满意的花开时刻就好。即使一时不顺,也要敞开胸怀。生命的荣枯并不是简单的重复,一时的得失不是成败的尺度。花开不是荣耀,而是一个美丽的结束,花谢也不是耻辱,而是一个低调的开始。

网易考拉推荐
GACHA精选

SAR成像入门(1)---如何模拟SAR数据  

2012-03-19 13:32:43|  分类: 学习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
%
% makeSARdata
%
%  make SAR project data
%
%  Written by M. A. Richards
%
clear all, close all, hold off
format compact
clc

c = 3e8;                               % velocity of light in m/sec


%######################################################
% User input section  #################################

% need a file name to store data, e.g. 'myfile' or 'temp' or 'sardata'.
% Data will then be in 'myfile.dat', etc.
file=input('Enter root file name for data and listing files: ','s');

DCR = 10;           % cross-range resolution, m
DR = 10;            % range resolution, m
Rcrp = 30000;        % range to swath center (crp: center reference point)
v = 150;             % platform velocity, m/s
fc = 10e9;           % RF frequency in Hz
lambda = c/fc;        % wavelength, m
tau = 5e-6;           % pulse length, seconds (bandwidth set by resolution)
Daz = 0.2;            % antenna azimuth size, m
thetaaz = lambda/Daz; % azimuth beamwidth, radians
BWdopp = 2*v/lambda*thetaaz;     % slow time Doppler bandwidth, Hz
Ls = 3000;            % swath depth, m
oversample_st = 3;    % slow time oversample factor; higher makes
                      % prettier pictures but larger data sets
oversample_ft = 3;    % fast time oversample factor, similar to slow time    

% Define target locations, one row of (x,R) coordinates per target
% Ranges are relative to the CRP range (Rcrp) above.

% coords = [0,0];         % a single point target at the CRP

coords = ...            % a grid of 9 point targets
    [-1000,-1000;
    -1000,0;
    -1000,+1000;
    0,-1000;
    0,0;
    0,+1000;
    +1000,-1000;
    +1000,0;
    +1000,+1000];


% End user input section #####################################
% ############################################################

% check consistency of cross-range resolution and aperture size
if (DCR < Daz/2)
    disp('Requested cross-range resolution below minimum');
    disp(['Requested cross-range resolution = ',num2str(DCR)]);
    disp(['Minimum cross-range resolution = (Daz/2) = ', num2str(Daz/2)]);
    disp('')
    return
end


% Compute required aperture time
Ta = lambda*Rcrp/2/v/DCR;
SAR数据模拟 - 完整源码 - 阿英 - Mr.Right
 Rc = (v*Ta)^2/8/Rcrp;        % nominal range curvature, m
SAR数据模拟 - 完整源码 - 阿英 - Mr.Right
 fprintf('\nRange resolution = %g meters',DR);
fprintf('\nCross-range resolution = %g meters',DCR);
fprintf('\nRange to central reference point = %g km',Rcrp/1e3);
fprintf('\nSwath length = %g m',Ls);
fprintf('\nAperture time = %g seconds',Ta);
fprintf('\nRange curvature at CRP = %g meters\n\n',Rc);

% Compute PRF, see if both range ambiguity and Doppler bandwidth
% constraints can be met

PRFgoal = oversample_st*BWdopp;     % PRF in Hz, 我以前不知道啊
PRF = min(PRFgoal,c/2/Ls);
if (PRF < PRFgoal)
    if (PRF < BWdopp)
        disp('No viable PRF');
        disp(['Minimum PRF (Doppler BW constraint) = ',num2str(BWdopp)]);
        disp(['Maximum PRF = (Swath ambiguity constraint) = ', ...
                num2str(2*Ls/c)]);
        disp('')
        return
    end
end
oversample_st = PRF/BWdopp;
PRI = (1/PRF);        % PRI in sec

% compute number of pulses, aperture positions

Npulses = round(Ta/PRI);    % # of pulses
jkl = 0:(Npulses-1);        % pulse index array
T_0 = PRI*jkl;              % relative start times of pulses, in sec
T_0 = T_0 - max(T_0)/2;     % recenters pulse times symmetrically +/- around zero
u = v*T_0;                  % aperture positions

% form radar chirp pulse

W = c/2/DR;       % chirp bandwidth, Hz
fs = oversample_ft*W;       % chirp sampling rate, Hz
s = git_chirp(tau,W,fs/W);
Ns = length(s);

fprintf('\nPulse length = %g microseconds\n',tau/1e-6)
fprintf('Chirp bandwidth = %g Mhz\n',W/1e6)
fprintf('Time-bandwidth product = %g\n',W*tau)
fprintf('Fast time sampling rate = %g Msamples/sec\n',fs/1e6)

figure(1)
plot((1e6/fs)*(0:length(s)-1),[real(s) imag(s)])
title('Real and Imaginary Parts of Chirp Pulse')
xlabel('time (usec)')
ylabel('amplitude')
grid

fprintf('\nWe are simulating %g pulses at an RF of %g GHz',Npulses,fc/1e9)
fprintf('\nand a PRF of %g kHz, giving a PRI of %g usec.\n\n',PRF/1e3,PRI/1e-6)

Ntargets = size(coords,1);

fprintf('\nThere are %g targets with the following locations relative to the CRP:',Ntargets)
for nt = 1:Ntargets
    fprintf('\n  x=%+8.1f m, R=%+8.1f m',coords(nt,1),coords(nt,2));
end
fprintf('\n\n')

% Precompute range vs. pulse number.  Radar coordinates given by x=u,r=0;
% target coordinates relative to the CRP are in array 'coords'.

range = zeros(Ntargets,Npulses);
for nt = 1:Ntargets
    for m = 1:Npulses
SAR成像入门(1)---如何模拟SAR数据 - 阿英 - Mr.Right
 
% 我今天在这个问题上纠结了整整一天,原因是没看出来散射点都在阴影部分的X-O-Yg平面上,Rcrp是椭圆长条的中心啊
range(nt,m) = sqrt((coords(nt,1)-u(m))^2 + (Rcrp+coords(nt,2))^2);
% 以上生成SAR数据的操作其实就是模型变换的过程,将模型坐标系中的坐标点转换到全局坐标系下的坐标
 
 % 我以前不知道Rcrp就是Rp,类似于到坐标原点, coords(nt,1) - u(m) 是x轴方向上P点到 x = u这条线的距离
% Rcrp 加上 coords(nt,2) 的效果就是相对于yg轴偏移coords(nt,2)
SAR成像入门(1)---如何模拟SAR数据 - 阿英 - Mr.Right
     end
end
figure(2)
plot(u,range)
xlabel('aperture position (m)')
ylabel('range (m)')
title('Range vs. Aperture Position for each Scatterer')
grid

Rmin = min(min(range));
Rmax = max(max(range));

% check to see if min and max range fall in the defined swath
if ( (Rmin < Rcrp-Ls/2) | (Rmax > Rcrp+Ls/2) )   % Rcrp就是swath椭圆长条的中心啊
    disp('Target ranges outside of swath!')
    disp(['Swath limits are ',num2str(Rcrp-Ls/2),' to ',num2str(Rcrp+Ls/2)])
    disp(['Min and max range are ',num2str(Rmin),' and ',num2str(Rmax)])
    return
end

% Define the range window
T_out = [2*(Rcrp-Ls/2)/c, 2*(Rcrp+Ls/2)/c+tau];  % 时间窗对应的距离窗

fprintf('\nThe range window limits are %g to %g usec.\n', ...
    T_out(1)/1e-6,T_out(2)/1e-6)
fprintf('\nThe range window starts at %g km.',(Rcrp-Ls/2)/1e3)
fprintf('\nThe range window ends at %g km.',(Rcrp+Ls/2)/1e3)
fprintf('\nThe unambiguous range interval is %g km.\n\n',c/2/PRF/1e3)


% Now build the fast time/slow time data matrix
disp('Now forming signal matrix')

delta_t = 1/fs;                        % sampling interval (sec)
t_y = [ T_out(1):delta_t:T_out(2) ]';  % output sampling times (sec)
T_p = Ns*delta_t;                      % length of input pulse (sec)

% ensure that all vectors are column vectors

s=s(:);  T_0=T_0(:);

% determine the quadratic phase modulation parameters for
% later interpolation of pulse samples

t_s = delta_t*[0:(Ns-1)]';
s_ph = unwrap(angle(s));   % s是线性调频chirp信号
warning off MATLAB:polyfit:RepeatedPointsOrRescale  % kills an annoying MATLAB warning about
% the polyfit function
q = polyfit(t_s,s_ph,2);    % 2次多项式拟合chirp信号的相位

% check result using correlation coefficient
sfit = polyval(q,t_s);   % sfit与s_ph应该接近
if (s_ph'*sfit)/norm(s_ph)/norm(sfit) < 0.99    % 相关系数是这样算的,我以前不会啊!
    disp('pulse phase is not linear or quadratic')
    disp('')
    return
end

%
%---  Form (initially empty) output matrix ---
%
Nr = length(t_y);       % output samples in a matrix
y = zeros(Nr,Npulses);    %Nr是快时间维度也叫距离维,Npulses是慢时间维度也叫多普勒维,横向距离维度

for m = 1:Npulses      % loop over aperture positions
    if (mod(m,100)==1)
        disp(['  ... processing pulse #',int2str(m),' of ',int2str(Npulses)])
    end
    for n = 1:Ntargets   % loop over targets
        r = range(n,m);
        
        % Compute start and end time of reflected pulse at receiver,
        % ensure that it falls entirely within the range (time) window
        tmin = 2*r/c; tmax = tmin + T_p;    % T_p是脉冲宽度, tmin和tmax原来是 时间观察窗口,即“时间窗”
        if ( (tmin < T_out(1)) | (tmax > T_out(2)) )
            fprintf('\nEcho from target #%g at range %g km',n,r/1e3)    % r 是距离窗
            fprintf('\nis COMPLETELY OUT OF the range window')
            fprintf('\non pulse #%g.\n',m)
            return
        end
        
        % Figure out which sample locations in the output grid contain
        % reflected pulse
        t_vals = t_y - tmin;   % t_y是快时间
        n_out = find( t_vals >= 0  &  t_vals < T_p );
        
        % Place pulse into output matrix.
        % Stop-and-hop assumed, so only phase modulation is due to geometry.
% All pulses have unit amplitude.
% There is no adjustment for R^4, either.
      
        y(n_out,m) = y(n_out,m) + ...
            exp( -j*2*pi*fc*tmin ) ...          % fc是载频,tmin在这里用,我以前不知道啊,tmin相当于一个时间起点
            .* exp( j*polyval(q,t_vals(n_out)) );
        
    end  % end of loop over targets
end    % end of loop over pulses

figure(3)
imagesc(1:Npulses,t_y,real(y))
xlabel('pulse number')
ylabel('fast time (sec)')
title('Real Part of Data Matrix')


% Save the data matrix in specified file.
% Save the student version in the mystery file.
% Also save all parameter value displays in corresponding file

data_file=[file,'.mat'];
mystery_file=[file,'_mys.mat'];
listing_file=[file,'.lis'];

eval(['save ',data_file,' c tau W fs s Npulses PRF PRI T_out fc', ...
        ' Ta v lambda Rc u Rcrp Ls DR DCR Ntargets coords y t_y Nr']);

eval(['save ',mystery_file,' c tau W fs s Npulses PRF PRI T_out fc', ...
        ' Ta v lambda Rc u Rcrp Ls DR DCR y t_y Nr']);

fid=fopen(listing_file,'w');   % listing file

fprintf(fid,['\rDESCRIPTION OF DATA IN FILE ',file,'.mat AND ',file,'_mys.mat\r\r']);
fprintf(fid,'\nRange resolution = %g meters',DR);
fprintf(fid,'\nCross-range resolution = %g meters',DCR);
fprintf(fid,'\nRange to central reference point = %g km',Rcrp/1e3);
fprintf(fid,'\nAperture time = %g seconds',Ta);
fprintf(fid,'\nRange curvature at CRP = %g meters\n\n',Rc);
fprintf(fid,'\rPulse length = %g microseconds\r',tau/1e-6);
fprintf(fid,'\rChirp bandwidth = %g Mhz\r',W/1e6);
fprintf(fid,'Time-bandwidth product = %g\n',W*tau);
fprintf(fid,'\rWe are simulating %g pulses at an RF of %g GHz',Npulses,fc/1e9);
fprintf(fid,'\r     and a PRF of %g kHz, giving a PRI of %g usec.',PRF/1e3,PRI/1e-6);
fprintf(fid,'\nRange resolution = %g meters',DR);
fprintf(fid,'\nCross-range resolution = %g meters',DCR);
fprintf(fid,'\nRange to central reference point = %g km',Rcrp/1e3);
fprintf(fid,'\nSwath length = %g m',Ls);
fprintf(fid,'\nRange curvature at CRP = %g meters\n\n',Rc);
fprintf(fid,'\rThe range window limits are %g to %g usec.\r', ...
    T_out(1)/1e-6,T_out(2)/1e-6);
fprintf(fid,'\rThe range window starts at %g km.',c*T_out(1)/2/1e3);
fprintf(fid,'\rThe range window ends at %g km.',c*T_out(2)/2/1e3);
fprintf(fid,'\nAperture time = %g seconds',Ta);
fprintf(fid,'\rThere are %g scatterers with the following coordinates:', ...
    Ntargets);
for nt = 1:Ntargets
    fprintf(fid,'\r  x=%+8.1f m, R=%+8.1f m',coords(nt,1),coords(nt,2) );
end

fclose(fid);

fprintf(['\n\nData is in file ',data_file])
fprintf(['\nStudent data is in file ',mystery_file])
fprintf(['\nListing is in file ',listing_file,'\n\n'])


% -------------- 用到的函数————————————————
function   x = git_chirp( T, W, p )
%CHIRP      generate a sampled chirp signal
%    X = git_chirp( T, W, <P> )
%      X:  N=pTW samples of a "chirp" signal
%            exp(j(W/T)pi*t^2)   -T/2 <= t < +T/2
%      T:  time duration from -T/2 to +T/2
%      W:  swept bandwidth from -W/2 to +W/2
%    optional:
%      P:  samples at P times the Nyquist rate (W)
%            i.e., sampling interval is 1/(PW)
%            default is P = 1
%
if nargin < 3
  p = 1;
end    
J = sqrt(-1);
%--------------
delta_t = 1/(p*W);
N = round( p*T*W );    %--same as T/delta_t, but rounded
nn = [0:N-1]';
% x = exp( J*pi*W/T * (delta_t*nn - T/2).^2 );  % old version
x = exp( J*pi*W/T * (delta_t*nn - (N-1)/2/p/W).^2 );  % symmetric version

  评论这张
 
阅读(1360)| 评论(2)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016