% ========================================================================
% Extended Chirp Scaling Algorithm
% Maurice Rueegg
% April 2005, May 2009
% ========================================================================
% Implementation of the basic algorithm for Extended Chirp Scaling as
% described in the paper:
% "Extended Chirp Scaling Algorithm for Air- and Spaceborne SAR Data
% Processing in Stripmap and ScanSAR Imaging Modes". A. Moreira,
% J. Mittermayer, R. Scheiber. IEEE Trans Geoscience and Remote Sensing,
% 34(5): 1123-1136, 1996
% Names of variables are according to the paper.
% ========================================================================
% Copyright (c) Maurice Rueegg (2005,2009)
% Permission to use these scripts for SAR processing and related research
% is granted. Please acknowledge the author.
% The author disclaims any responsibility for or liability related to this
% Software. Any questions or claims related to the Software can be
% directed to info@znuerb.com 
% ========================================================================

clear all; close all;
doTimePlots = 1; % 1 to draw plots, 0 to skip all but the inital and final plots
doSpectralPlots = 0; % 1 to draw spectral plots, 0 to skip them

% ========================================================================
% set parameters, read parameters
% ========================================================================
% example raw data point target from file 'pointtarget.raw'
% the file contains a simulated point target with Envisat ASAR
% system parameters:
c = 299792458; % speed of light in vacuum
f_c = 5331004416; % carrier frequency
f_s = 19207680; % data sampling rate
echoes = 2206; % number of radar echoes in data file
samples = 672; % number of samples per radar echo
t_near = 0.00533035424126647; % near range fast time (2 x range)
tau = t_near:1/f_s:t_near+samples/f_s-1/f_s; % fast time
f_dc = 0; % Doppler centroid (squint angle)
v = 7500; % SAR satellite platform velocity
PRF = 2067.120103315; % pulse repetition frequency
t_p = 2.160594095e-005; % chirp pulse duration
B = 16000000; % chirp bandwidth

% additional paramters, see paper
r_ref = (tau(1)+samples/2/f_s)/2*c;
alpha = 1;
f_a = -PRF/2+f_dc:PRF/echoes:f_dc+PRF/2-PRF/echoes;
f_r = -f_s/2:f_s/samples:f_s/2-f_s/samples;
lambda = c/f_c;

% read the raw data
data=readMatrix('pointtarget.raw',1,samples,1,echoes);

G=abs(data);
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
figure;
colormap(gray(256))
image(256-cg*(G-ng));
axis('image');axis('xy')
xlabel('range')
ylabel('azimuth')
title('SAR raw data')

% ========================================================================
% azimuth fft
% ========================================================================
data = ftx(data);

if (doSpectralPlots==1)
    G=angle(data);
    xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
    figure;
    colormap(jet(256))
    image(256-cg*(G-ng));
    axis('image');axis('xy')
    xlabel('range')
    ylabel('Doppler frequency phase')
    title('range signal/Doppler domain')
end;

% ========================================================================
% chirp scaling, range scaling: H1
% ========================================================================
beta = (1 - (f_a*lambda/2/v).^2).^0.5;
a = 1./beta - 1;
R = r_ref./beta;
a_scl = a + (1-alpha).*(1+a)./alpha;
k_r = -B./t_p;
% WARNING: in the work of Moreira et al, 1996 k_r is defined negative
% because they assume a "down-chirp", whereas we assume an "up-chirp" and
% must explicitely set k_r as negative for the following transfer functions
% to be equal to the theoretical formulations.

k_inv = 1./k_r - (2.*lambda.*r_ref.*(beta.^2-1))./(c^2.*beta.^3);
k = 1./k_inv;
%k = (k_r*c.^2.*beta.^3) ./ (c.^2.*beta.^3 - ...
%k_r.*2.*lambda.*r_ref.*(beta.^2-1)); % same thing as above

x = k.*a_scl;
X = x(:)*ones(1,samples);
Tau = ones(echoes,1)*tau;
y = 2.*R./c;
Y = y(:)*ones(1,samples);
Z = (Tau-Y).^2;
H1 = exp(-j*pi*X.*Z);

data = data .* H1;

if (doSpectralPlots==1)
    G=angle(data);
    xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
    figure;
    colormap(jet(256))
    image(256-cg*(G-ng));
    axis('image');axis('xy')
    xlabel('range')
    ylabel('Doppler frequency phase')
    title('chirp scaling and range scaling')
end;

if (doTimePlots==1)
    figure;
    colormap(jet(256))
    mesh(abs(iftx(data)));
    xlabel('range')
    ylabel('azimuth')
    zlabel('intensity')
    title('after H1')
    axis('tight');
end;

% ========================================================================
% range fft
% ========================================================================
data = fty(data);

if (doSpectralPlots==1)
    G=angle(data);
    xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
    figure;
    colormap(jet(256))
    image(256-cg*(G-ng));
    axis('image');axis('xy')
    xlabel('signal frequency phase')
    ylabel('Doppler frequency phase')
    title('2D spectral domain')
end;

% ========================================================================
% bulk rcmc, range compression: H2
% ========================================================================
x = 1./(k.*(1+a_scl));
X = x(:)*ones(1,samples);
y = f_r.^2;
Y = ones(echoes,1)*y;
z = f_r;
Z = ones(echoes,1)*z;
A = a(:)*ones(1,samples);
Z = Z .* A;
H2 = exp(-j*pi*X.*Y) .* exp(j*4*pi*r_ref/c.*Z);
data = data .* H2;

if (doSpectralPlots==1)
    G=angle(data);
    xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
    figure;
    colormap(jet(256))
    image(256-cg*(G-ng));
    axis('image');axis('xy')
    xlabel('signal frequency phase')
    ylabel('Doppler frequency phase')
    title('bulk rcmc and range compression')
end;

if (doTimePlots==1)
    figure;
    colormap(jet(256))
    mesh(abs(iftx(ifty(data))));
    xlabel('range')
    ylabel('azimuth')
    zlabel('intensity')
    title('after H2')
    axis('tight');
end;

% ========================================================================
% range ifft
% ========================================================================
data = ifty(data);

if (doSpectralPlots==1)
    G=angle(data);
    xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
    figure;
    colormap(jet(256))
    image(256-cg*(G-ng));
    axis('image');axis('xy')
    xlabel('range')
    ylabel('Doppler frequency phase')
    title('range Doppler Domain after range compression')
end;

% ========================================================================
% angle correction: H3
% ========================================================================
r_0 = tau/2*c;
x = k.*a_scl.*(1+a).^2 ./ (c^2.*(1+a_scl));
X = x(:)*ones(1,samples);
z = (r_0-r_ref).^2;
Z = ones(echoes,1)*z;
dphi = 4*pi*X.*Z;
H3 = exp(j*dphi);

data = data .* H3;

if (doSpectralPlots==1)
    G=angle(data);
    xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
    figure;
    colormap(jet(256))
    image(256-cg*(G-ng));
    axis('image');axis('xy')
    xlabel('range')
    ylabel('Doppler frequency phase')
    title('range Doppler domain after angle correction')
end;

if (doTimePlots==1)
    figure;
    colormap(jet(256))
    mesh(abs(iftx(data)));
    xlabel('range')
    ylabel('azimuth')
    zlabel('intensity')
    title('after H3')
    axis('tight');
end;

% ========================================================================
% azimuth compression: H4
% ========================================================================
r_0_scl = r_ref + (r_0-r_ref)/alpha;
X = ones(echoes,1)*r_0_scl;
Z = (beta(:)-1)*ones(1,samples);
H4 = exp(j*4*pi/lambda*X.*Z);

data = data .* H4;

if (doSpectralPlots==1)
    G=angle(data);
    xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
    figure;
    colormap(jet(256))
    image(256-cg*(G-ng));
    axis('image');axis('xy')
    xlabel('range')
    ylabel('Doppler frequency phase')
    title('azimuth compression')
end;

if (doTimePlots==1)
    figure;
    colormap(jet(256))
    mesh(abs(iftx(data)));
    xlabel('range')
    ylabel('azimuth')
    zlabel('intensity')
    title('focused (after H4)')
    axis('tight');
end;
% ========================================================================
% azimuth ifft
% ========================================================================
data = iftx(data);

% ========================================================================
% write result
% ========================================================================

G=abs(data);
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
figure;
colormap(gray(256))
image(256-cg*(G-ng));
axis('image');axis('xy')
xlabel('range')
ylabel('azimuth')
title('focused data')

% = end ==================================================================
