sar_matlabcode.m

SAR (synthetic aperture radar) image processing matlab code

The code in MATLAB referred to:
%****************************************************
function [cref,fcref]=rng_ref(nfft,fs,pulsedur,slope)
%
% routine to compute ERS chirp and its fourier transform
%
% input
% fs - sampling frequency, ts=1./fs
% pulsedur - pulse duration
% slope - chirp slope
%
% set the constants and make npts be odd
%
npts=floor(fs*pulsedur);
ts=1./fs;
if(mod(npts,2.0) == 0.0)
npts=npts+1;
end
%
% compute the reference function
%
npt2=floor(npts/2.);
t=ts*(-npt2:npt2);
phase=pi*slope*t.*t;
cref1=exp(i*phase);
%
% pad the reference function to nfft
%
cref=[cref1,zeros(1,nfft-npts)]';
%
% compute the fourier transform
%
fcref=fft(cref)/nfft;
%****************************************************
function [cazi,fcazi]=azi_ref(nazi,PRF,fdc,fr)
%
% routine to compute ERS azimuthal chirp and its fourier transform
%
% input
% nazi - number of points in azimuth
% PRF - pulse repitition frequency, ts=1./fs
% fdc - doppler centriod frequency
% fr - doppler frequency rate
%
% set the constants and make npts be odd
%
npts=min(nazi-1,1296);
ts=1./PRF;
if(mod(npts,2.0) == 0.0)
npts=npts+1;
end
%
% compute the azimuth chirp function
%
npt2=floor(npts/2.);
t=ts*(-npt2:npt2);
phase=-2.*pi*fdc*t+pi*fr*t.*t;
cazi1=exp(i*phase);
%
% pad the function to nfft
%
cazi=[cazi1(npt2:npts),zeros(1,nazinpts),
cazi1(1:npt2-1)]';
%
% compute the fourier transform
%
fcazi=fft(cazi)/nazi;
%****************************************************
function [csar,nrow,ncol]=read_rawsar(sar_file)
%
% routine to read and unpack ERS SAR data in DPAF format
%
fid=fopen(sar_file,'r');
%
% read the bytes
%
[sar,nsar]=fread(fid,'uchar');
sar=reshape(sar,11644,nsar/11644);
nrow=nsar/11644;
ncol=5616;
st=fclose(fid);
%
% extract the real and imaginary parts
% and remove the mean value
%
csar=complex(sar(413:2:11643,:)-15.5,sar(414:2:11644,:)-15.5);
%
%****************************************************
% matlab script to focus ERS-2 signal data
%
% set some constants for e2_10001_2925
%
% range parameters
%
rng_samp_rate = 1.896e+07;
pulse_dur = 3.71e-05;
chirp_slope = 4.1779e+11;
%
% azimuth parameters
%
PRF=1679.902394;
radar_wavelength=0.0566666;
SC_vel=7125.;
%
% compute the range to the radar reflectors
%
near_range=829924.366;
dr=3.e08/(2.*rng_samp_rate);
range=near_range+2700*dr;
%
% use the doppler centroid estimated from the data and the
% doppler rate from the spacecraft velocity and range
%
fdc=284;
fr=2*SC_vel*SC_vel/(range*radar_wavelength);
%
% get some sar data
%
[cdata,nrow,ncol] = read_rawsar('out2.raw');
%
% make a colormap
%
map=ones(21,3);
for k=1:21;
level=0.05*(k+8);
level=min(level,1);
map(k,:)=map(k,:).*level;
end
colormap(map);
%
% image the raw data
%
figure(1)
subplot(2,2,1),imagesc(abs(cdata'));
xlabel('range')
ylabel('azimuth')
title('unfocussed raw data')
axis([2600,2900,1000,1200])
%
% generate the range reference function
%
[cref,fcref]=rng_ref(ncol,rng_samp_rate,pulse_dur,chirp_slope);
%
% take the fft of the SAR data
%
fcdata=fft(cdata);
%
% multiply by the range reference function
%
cout=0.*fcdata;
for k=1:nrow;
ctmp=fcdata(:,k);
ctmp=fcref.*ctmp;
cout(:,k)=ctmp;
end
clear cdata
%
% now take the inverse fft
%
odata=ifft(cout);
clear cout
%
% plot the image and the reflector locations
%
x0=[2653.5,2621];
x0=x0+90;
y0=[20122,20226];
y0=y0-19500+427;
figure(1)
hold
subplot(2,2,2),imagesc(abs(odata'));
plot(x0,y0,'o')
xlabel('range')
ylabel('azimuth')
title('range compressed')
axis([2600,2900,1000,1200])
%
% use this for figure 2 as well
%
figure(2)
colormap(map);
subplot(2,2,1),imagesc(abs(odata'));
hold on
plot(x0,y0,'o')
xlabel('range')
ylabel('azimuth')
title('range compressed')
axis([2600,2900,1000,1200])
%
% generate the azimuth reference function
%
[cazi,fcazi]=azi_ref(nrow,PRF,fdc,fr);
%
% take the column-wise fft of the range-compressed data
%
fcdata=fft(odata');
%
% multiply by the azimuth reference function
%
cout=0.*fcdata;
for k=1:ncol;
ctmp=fcdata(:,k);
ctmp=fcazi.*ctmp;
cout(:,k)=ctmp;
end
%
% now take the inverse fft and plot the data
%
odata=ifft(cout);
figure(2)
subplot(2,2,2),imagesc(abs(odata));
hold on
plot(x0,y0,'o')
xlabel('range')
ylabel('azimuth')
title('range and azimuth compressed')
axis([2600,2900,1000,1200])