cubo.m

rang doppler imaging and motion compensation中的源代码

function cubo(w,vt,NI,n)

%   The following file generates images using
%   subapertures and the Gabor Wavelet Transform method.
%
%   Additional files needed are:
%    Genttf02.m
%	  Wtv02img.m
%
%   cubo(w,vt,NI,n)
%
%   w = angular velocity
% Try w = 0.1
%   vt = translational velocity
% Try values close to 0 (i.e. vt = 0.02)   
%   NI = Number of images to be obtained using subapertures
%        (i.e., number of subapertures)
% Try NI = 4
%   n = width of Gaussian window in gabor transform
%       standard deviation
% Try n = 12

tiempo0 = clock;
c=3.0e8;
f0=3.0e9;
bw=225.0e6;		%.5*225.0e6;
pulses=64;
burst=256;
T1=(pulses-1)/bw; 	%pulse duration
PRF= 5.4e3;		% 2*5.4e3;
T2=1/PRF; 		%PRI pulse repetition interval
%w=0.1;			% original value 0.10; 
%vt=10;
R0=0.0e3;
A=1.0;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sampling times Si
% Sik array of 64x64 elements
%

i = 1:pulses*burst;
Sik = T1/2+2*R0/c+(i-1)*T2;

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Multiply sampling times Si by angular velocity w 
% new twist: make angular velocity variable from .001 to .04 in steps of
% 9.5215e-6
%
  wt = Sik*w;
  t_vel = Sik*vt;

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

 i = 1:pulses;
 deltaf = (i-1)*1/T1;

% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

factor(burst,pulses) = 0;
k = ( -4*pi*(f0+deltaf) )/c;
for i = 1:burst;
   factor(i,:) = k;
end
   
   
%		T25602.M
% Synthetic ISAR Processing or Stepped-Frequency imaging
% Summary of method:
% 1.  Use stepped-freq waveforms to obtain time history of targets reflectivity
%     sampled in frequency domain.
% 2.  Correct spectral samples for radar system phase and amplitude ripple.
% 3.  Correction of spectral samples for velocity produced distortion.
% 4.  Conversion of the time history of target s corrected samples of freq domain
%     reflectivity into a time history of the target s synthetic range profile
% 5.  Conversion of target s range profile history into a slant range versus
%     cross range image.
% 6.  Moving window image processing to obtain an animated real-time image display
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 	deltaf - freq step
% 	c - 	speed of light c=3.0e8 m/s
% 	fi - 	frequency of pulse i(Hz), RF freq. inside pulse is varied. 
% 	w - 	angular velocity (rad/sec)
% 	R - 	initial range (m)
% 	N - 	number of pulses in kth burst
% 	T1 - 	pulse width
% 	T2 - 	Pulse Repetition Interval
% 	Sik - 	sampling time  Sik=T1/2+2*R0/c+i*T2+k*n*T2, i=pulse index
%
%	Generating positions or values of R (range)
% 	R=v*Sik, instantaneous position Ri=R0-v*Sik
% 	Sampling times for first burst
%

col=64;
row=256;
g(256,64)=0.0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Include angle theta1 for scatterers from wt term, t=Sik
% 
%

% disp('1.- Generating arg1,arg2,... arrays.')
% scatterers order 1,2,3...
theta=[0.0,1.64,3.57,5.18,6.84,7.96,10.61,12.87,15.37,19.79,26.56...
      ,36.86,51.34,70.01,90.00,101.65,110.69,115.03,115.98,110.66,99.76...
      ,106.48,113.31,118.07,124.58,124.81,131.30,132.53,136.70,141.31...
      ,147.27,157.53,162.07,163.98,165.97,160.80,157.03,154.22,156.42...
      ,159.33,164.62,170.50,169.93,175.22,176.62...
	,-1.64,-3.57,-5.18,-6.84,-7.96,-10.61,-12.87,-15.37,-19.79,-26.56...
      ,-36.86,-51.34,-70.01,-90.00,-101.65,-110.69,-115.03,-115.98,-110.66,-99.76...
      -106.48,-113.31,-118.07,-124.58,-124.81,-131.30,-132.53,-136.70,-141.31...
      ,-147.27,-157.53,-162.07,-163.98,-165.97,-160.80,-157.03,-154.22,-156.42...
      ,-159.33,-164.62,-170.50,-169.93,-175.22,-176.62];
theta=theta+(90.0);

%
r=[13.5,12.8,11.85,11.0,10.0,9.0,8.0,7.125,6.25,5.25,4.45,3.75,3.125,2.85...
   ,3.70,4.85,6.20,7.45,8.90,8.55,8.625,9.125,9.80,10.20,10.75,11.45,10.65...
   ,9.5,8.75,8.0,7.40,7.85,6.50,7.25,8.25,9.125,10.25,11.50,12.50,12.75...
   ,12.45,12.125,10.875,12.0,12.75,...
   12.8,11.85,11.0,10.0,9.0,8.0,7.125,6.25,5.25,4.45,3.75,3.125,2.85...
   ,3.70,4.85,6.20,7.45,8.90,8.55,8.625,9.125,9.80,10.20,10.75,11.45,10.65...
   ,9.5,8.75,8.0,7.40,7.85,6.50,7.25,8.25,9.125,10.25,11.50,12.50,12.75...
   ,12.45,12.125,10.875,12.0,12.75];
r=r.*1.0; %2.5;

%
theta_rad=theta.*0.017455329;  	   	% 0.017455329 rad/deg
%


for scat=1:3:89
     
       arg1 = (wt - theta_rad(scat) );
       rngterm1 = R0 - t_vel - r(scat)*sin(arg1);
       
       range1 = reshape(rngterm1,pulses,burst);
       range1 = range1.';
   
	    clear rngterm1;
%
	    phase1 = factor.* range1;
        clear range1;
%
	g1=exp(j*phase1);
%
	clear phase1;
	g=g+g1;
        clear g1;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Matrix G is the time history of the target's reflectivity 
%     measured at discrete frequencies in the RF domain. Complex echo 
%     samples are represented by amplitude and phase. every complex echo 
%     sampled is collected from each target echo pulse (one sample per 
%     transmitted pulse).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     The inverse Fourier transform from each of the N bursts of n freq 
%     samples is taken, row X row to generate N new rows of synthetic range
%     profiles each with n synyhetic range cells, this matrix is called the time-
%     history of the target's range profilesignature.
%     Generate inverse FFT of each vector above
%

th_256=ifft(g');
th_256=th_256';
org_img256 = fft(th_256);
org_img256 = fftshift(org_img256);
org_img256 = abs(org_img256);
org_img256 = org_img256/max(max(org_img256));
org_img256 = (1-org_img256);
colormap(gray(128))
X=0:42/63:42;
Y=0:48/63:48;
image(X,Y,(128*org_img256));
title(['ISAR image generated via the Fourier transform', ' w = ',num2str(w),' vt = ',num2str(vt)] )
ylabel('cross range (m)')

% Compute Minimum Entropy

mec = fft(th_256);
mec = fftshift(mec);
mec = abs(mec);
mec = mec/(sum(sum(mec)));
mec = mec.*exp(1-mec);
mec = sum(sum(mec));

xlabel(['slant range (m).    MEV = ',num2str(mec)])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate same image as above but using a Hamming window
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure
[gs1 gs2] = size(g);

ven = hamming(gs1);
for i=1:gs2
   gh(:,i) = g(:,i).*ven;
end

th_256h=ifft(gh');
th_256h=th_256h';

[gs1 gs2] = size(th_256h);
ven = hamming(gs1);
for i=1:gs2
   gh(:,i) = g(:,i).*ven;
end

org_img256h = fft(th_256h);
org_img256h = fftshift(org_img256h);
org_img256h = abs(org_img256h);
org_img256h = org_img256h/max(max(org_img256h));
org_img256h = (1-org_img256h);
colormap(gray(128))
X=0:42/63:42;
Y=0:48/63:48;
image(X,Y,(128*org_img256h));
title(['ISAR image generated via the Fourier transform and Hamming window', ' w = ',num2str(w),' vt = ',num2str(vt)] )
ylabel('cross range (m)')

% Compute Minimum Entropy

mec = fft(th_256);
mec = fftshift(mec);
mec = abs(mec);
mec = mec/(sum(sum(mec)));
mec = mec.*exp(1-mec);
mec = sum(sum(mec));

xlabel(['slant range (m).    MEV = ',num2str(mec)])

clear org_img256h gh ven gs1 gs2 th_256h

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate a sequence of ISAR images using the Fourier transform
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

th_256=ifft(g');
th_256=th_256';

len = floor(256/NI);
for i = len:len:len*NI
     org_img256 = fft(th_256(((i-len+1):i),:));
     org_img256 = fftshift(org_img256);
     org_img256 = abs(org_img256);
     org_img256 = org_img256/max(max(org_img256));
     org_img256 = (1-org_img256);
     colormap(gray(128))
     figure
     image((128*org_img256));
     title(['ISAR image generated via the Fourier transform and using subapertures', ' w = ',num2str(w),' vt = ',num2str(vt)] )
     xlabel('slant range')
     ylabel('cross range')
     % Compute Minimum Entropy
     mec = fft(th_256(((i-len+1):i),:));
     mec = fftshift(mec);
     mec = abs(mec);
     mec = mec/(sum(sum(mec)));
     mec = mec.*exp(1-mec);
     mec = sum(sum(mec));
     xlabel(['slant range (m).    MEV = ',num2str(mec)])
end

genttf02
%
%		END  OF T25602.M 
%