自制合成孔径雷达(5) SAR代码解读

在lecture 4的ppt里，有SAR的处理框图。SAR的matlab代码跟这个基本相符。

matlab代码分为两个文件，SBAND_RMA_opendata是读取数据，读取完了调用SBAND_RMA_IFP处理数据。

有一点要注意，默认do_all_plots是0，很多绘制中间数据的代码没有真的用上。

下面是SBAND_RMA_opendata.m

%-------------------------------------------%
%Process raw data here
clear all;
close all;

%read the raw data .wave file here
[Y,FS,NBITS] = wavread('towardswarehouse.wav');

%constants
c = 3E8; %(m/s) speed of light

%radar parameters
Tp = 20E-3; %(s) pulse time
Trp = 0.25; %(s) min range profile time duration
N = Tp*FS; %# of samples per pulse
fstart = 2260E6; %(Hz) LFM start frequency
fstop = 2590E6; %(Hz) LFM stop frequency
%fstart = 2402E6; %(Hz) LFM start frequency for ISM band
%fstop = 2495E6; %(Hz) LFM stop frequency for ISM band
BW = fstop-fstart; %(Hz) transmti bandwidth
f = linspace(fstart, fstop, N/2); %instantaneous transmit frequency

%the input appears to be inverted
trig = -1*Y(:,1);
s = -1*Y(:,2);
clear Y;

%parse data here by position (silence between recorded data)
rpstart = abs(trig)>mean(abs(trig));
count = 0;
Nrp = Trp*FS; %min # samples between range profiles

for ii = Nrp+1:size(rpstart,1)-Nrp
    if rpstart(ii) == 1 & sum(rpstart(ii-Nrp:ii-1)) == 0
        count = count + 1;
        RP(count,:) = s(ii:ii+Nrp-1);
        RPtrig(count,:) = trig(ii:ii+Nrp-1);
    end
end

%parse data by pulse
count = 0;
thresh = 0.08;
clear ii;
for jj = 1:size(RP,1)
    %clear SIF;
    SIF = zeros(N,1);
    start = (RPtrig(jj,:)> thresh);
    count = 0;
    jj
    for ii = 12:(size(start,2)-2*N)
        [Y I] =  max(RPtrig(jj,ii:ii+2*N));
        if mean(start(ii-10:ii-2)) == 0 & I == 1
            count = count + 1;
            SIF = RP(jj,ii:ii+N-1)' + SIF;
        end
    end
    %hilbert transform
    q = ifft(SIF/count);
    sif(jj,:) = fft(q(size(q,1)/2+1:size(q,1)));
end
sif(find(isnan(sif))) = 1E-30; %set all Nan values to 0

%SAR data should be ready here
clear s;
s = sif;
save routsidewarehouse2 s; %for image data

%-------------------------------------------%
%load additional varaibles and setup constants for radar here
clear all;
c = 3E8; %(m/s) speed of light

%load IQ converted data here
load routsidewarehouse2 s; %load variable sif %for image data

for ii = 1:size(s,1)
    s(ii,:) = s(ii,:) - mean(s,1);
end

%sif = s-sif_sub; %perform coherent background subtraction
%sif = sif_sub; %image just the background
sif = s; %image without background subtraction
clear s;
clear sif_sub;

%***********************************************************************
%radar parameters
fc = (2590E6 - 2260E6)/2 + 2260E6; %(Hz) center radar frequency
B = (2590E6 - 2260E6); %(hz) bandwidth
cr = B/20E-3; %(Hz/sec) chirp rate
Tp = 20E-3; %(sec) pulse width
%VERY IMPORTANT, change Rs to distance to cal target
%Rs = (12+9/12)*.3048; %(m) y coordinate to scene center (down range), make this value equal to distance to cal target
Rs = 0;
Xa = 0; %(m) beginning of new aperture length
delta_x = 2*(1/12)*0.3048; %(m) 2 inch antenna spacing
L = delta_x*(size(sif,1)); %(m) aperture length
Xa = linspace(-L/2, L/2, (L/delta_x)); %(m) cross range position of radar on aperture L
Za = 0;
Ya = Rs; %THIS IS VERY IMPORTANT, SEE GEOMETRY FIGURE 10.6
t = linspace(0, Tp, size(sif,2)); %(s) fast time, CHECK SAMPLE RATE
Kr = linspace(((4*pi/c)*(fc - B/2)), ((4*pi/c)*(fc + B/2)), (size(t,2)));

%Save background subtracted and callibrated data
save sif sif delta_x Rs Kr Xa;
%clear all;

%run IFP
SBAND_RMA_IFP;

 

它所做的与上一个测距类似，读取每一段方波打开时的数据（脉冲）。

for ii = Nrp+1:size(rpstart,1)-Nrp
    if rpstart(ii) == 1 & sum(rpstart(ii-Nrp:ii-1)) == 0
        count = count + 1;
        RP(count,:) = s(ii:ii+Nrp-1);
        RPtrig(count,:) = trig(ii:ii+Nrp-1);
    end
end

读完以后还把这些脉冲数据加在一起

for ii = 12:(size(start,2)-2*N)
        [Y I] =  max(RPtrig(jj,ii:ii+2*N));
        if mean(start(ii-10:ii-2)) == 0 & I == 1
            count = count + 1;
            SIF = RP(jj,ii:ii+N-1)' + SIF;
        end
    end

然后做希尔伯特变换

    %hilbert transform
    q = ifft(SIF/count);
    sif(jj,:) = fft(q(size(q,1)/2+1:size(q,1)));

后续算法都在SBAND_RMA_IFP.m中完成。

%***********************************************************************
%a note on formatting, our convention is sif(Xa,t)

clear all;
load sif;

%apply hanning window to data first
N = size(sif,2);
for ii = 1:N
    H(ii) = 0.5 + 0.5*cos(2*pi*(ii-N/2)/N);
end

for ii = 1:size(sif,1)
    sif_h(ii,:) = sif(ii,:).*H;
end
sif = sif_h;

figcount = 1;
close_as_you_go = 0;
do_all_plots = 0;

set(0,'defaultaxesfontsize',13); %set font size on plots so we can see it in the dissertation

% NOTE: the function 'dbv.m' is just dataout = 20*log10(abs(datain));
%***********************************************************************
if do_all_plots == 1,
    figure(figcount);
    S_image = angle(sif);
    imagesc(Kr, Xa, S_image);
    colormap(gray);
    title('Phase Before Along Track FFT');
    xlabel('K_r (rad/m)');
    ylabel('Synthetic Aperture Position, Xa (m)');
    cbar = colorbar;
    set(get(cbar, 'Title'), 'String', 'radians','fontsize',13); 
    print(gcf, '-djpeg100', 'phase_before_along_track_fft.jpg');
    if close_as_you_go == 1,
        close(figcount);
    end
    figcount = figcount + 1;
end

%along track FFT (in the slow time domain)
%first, symetrically cross range zero pad so that the radar can squint
zpad = 2048; %cross range symetrical zero pad
szeros = zeros(zpad, size(sif,2));
for ii = 1:size(sif,2)
    index = round((zpad - size(sif,1))/2);
    szeros(index+1:(index + size(sif,1)),ii) = sif(:,ii); %symetrical zero pad
end
sif = szeros;
clear ii index szeros;


S = fftshift(fft(sif, [], 1), 1);
clear sif;
Kx = linspace((-pi/delta_x), (pi/delta_x), (size(S,1)));

if do_all_plots == 1,
    figure(figcount);
    S_image = dbv(S);
    imagesc(Kr, Kx, S_image, [max(max(S_image))-40, max(max(S_image))]);
    colormap(gray);
    title('Magnitude After Along Track FFT');
    xlabel('K_r (rad/m)');
    ylabel('K_x (rad/m)');
    cbar = colorbar;
    set(get(cbar, 'Title'), 'String', 'dB','fontsize',13); 
    print(gcf, '-djpeg100', 'mag_after_along_track_fft.jpg');
    if close_as_you_go == 1,
        close(figcount);
    end
    figcount = figcount + 1;
end

if do_all_plots == 1,
    figure(figcount);
    S_image = angle(S);
    imagesc(Kr, Kx, S_image);
    colormap(gray);
    title('Phase After Along Track FFT');
    xlabel('K_r (rad/m)');
    ylabel('K_x (rad/m)');
    cbar = colorbar;
    set(get(cbar, 'Title'), 'String', 'radians','fontsize',13); 
    print(gcf, '-djpeg100', 'phase_after_along_track_fft.jpg');
    if close_as_you_go == 1,
        close(figcount);
    end
    figcount = figcount + 1;
end

if do_all_plots == 1,
    figure(figcount);
    S_image = dbv(fftshift(fft(S, [], 2), 2));
    imagesc(linspace(-0.5, 0.5, size(S, 2)), Kx, S_image, [max(max(S_image))-40, max(max(S_image))]);
    colormap(gray);
    title('Magnitude of 2-D FFT of Input Data');
    xlabel('R_{relative} (dimensionless)');
    ylabel('K_x (rad/m)');
    cbar = colorbar;
    set(get(cbar, 'Title'), 'String', 'dB','fontsize',13); 
    print(gcf, '-djpeg100', 'mag_after_2D_fft.jpg');
    if close_as_you_go == 1,
        close(figcount);
    end
    figcount = figcount + 1;
end


%**********************************************************************
%matched filter

%create the matched filter eq 10.8
for ii = 1:size(S,2) %step thru each time step row to find phi_if
    for jj = 1:size(S,1) %step through each cross range in the current time step row
        %phi_mf(jj,ii) = -Rs*Kr(ii) + Rs*sqrt((Kr(ii))^2 - (Kx(jj))^2);
        phi_mf(jj,ii) = Rs*sqrt((Kr(ii))^2 - (Kx(jj))^2);
        Krr(jj,ii) = Kr(ii); %generate 2d Kr for plotting purposes
        Kxx(jj,ii) = Kx(jj); %generate 2d Kx for plotting purposes
    end
end
smf = exp(j*phi_mf); %%%%%%%%%%%%

%note, we are in the Kx and Kr domain, thus our convention is S_mf(Kx,Kr)

%appsly matched filter to S
S_mf = S.*smf;

clear smf phi_mf;

if do_all_plots == 1,
    figure(figcount);
    S_image = angle(S);
    imagesc(Kr, Kx, S_image);
    colormap(gray);
    title('Phase After Matched Filter');
    xlabel('K_r (rad/m)');
    ylabel('K_x (rad/m)');
    cbar = colorbar;
    set(get(cbar, 'Title'), 'String', 'radians','fontsize',13); 
    print(gcf, '-djpeg100', 'phase_after_matched_filter.jpg');
    if close_as_you_go == 1,
        close(figcount);
    end
    figcount = figcount + 1;
end

clear S;

if do_all_plots == 1,
    figure(figcount);
    S_image = dbv(fftshift(fft(S_mf, [], 2), 2));
    imagesc(linspace(-0.5, 0.5, size(S_mf, 2)), Kx, S_image, [max(max(S_image))-40, max(max(S_image))]);
    colormap(gray);
    title('Magnitude of 2-D FFT of Matched Filtered Data');
    xlabel('R_{relative} (dimensionless)');
    ylabel('K_x (rad/m)');
    cbar = colorbar;
    set(get(cbar, 'Title'), 'String', 'dB','fontsize',13); 
    print(gcf, '-djpeg100', 'mag_after_downrange_fft_of_matched_filtered_data.jpg');
    if close_as_you_go == 1,
        close(figcount);
    end
    figcount = figcount + 1;
end



%**********************************************************************
%perform the Stolt interpolation

%FOR DATA ANALYSIS
kstart =73;
kstop = 108.5;

%kstart = 95;
%kstop = 102;


Ky_even = linspace(kstart, kstop, 1024); %create evenly spaced Ky for interp for real data

clear Ky S_St;
%for ii = 1:size(Kx,2)
count = 0;
for ii = 1:zpad;
%for ii = round(.2*zpad):round((1-.2)*zpad)
    count = count + 1;
    Ky(count,:) = sqrt(Kr.^2 - Kx(ii)^2);
    %S_st(ii,:) = (interp1(Ky(ii,:), S_mf(ii,:), Ky_even)).*H;
    S_st(count,:) = (interp1(Ky(count,:), S_mf(ii,:), Ky_even));
end
S_st(find(isnan(S_st))) = 1E-30; %set all Nan values to 0
clear S_mf ii Ky;

if do_all_plots == 1,
    figure(figcount);
    S_image = angle(S_st);
    imagesc(Ky_even, Kx, S_image);
    %imagesc(S_image);
    colormap(gray);
    title('Phase After Stolt Interpolation');
    xlabel('K_y (rad/m)');
    ylabel('K_x (rad/m)');
    cbar = colorbar;
    set(get(cbar, 'Title'), 'String', 'radians','fontsize',13); 
    print(gcf, '-djpeg100', 'phase_after_stolt_interpolation.jpg');
    if close_as_you_go == 1,
        close(figcount);
    end
    figcount = figcount + 1;
end

%apply hanning window to data, cleans up data ALOT
N = size(Ky_even,2);
for ii = 1:N
    H(ii) = 0.5 + 0.5*cos(2*pi*(ii-N/2)/N);
end

for ii = 1:size(S_st,1)
    S_sth(ii,:) = S_st(ii,:).*H;
end
%S_st = S_sth;

%*********************************************************************
%perform the inverse FFT's
%new notation:  v(x,y), where x is crossrange
%first in the range dimmension
clear v Kr Krr Kxx Ky_even;
v = ifft2(S_st,(size(S_st,1)*4),(size(S_st,2)*4));

%bw = (3E8/(4*pi))*(max(xx)-min(xx));
bw = 3E8*(kstop-kstart)/(4*pi);
max_range = (3E8*size(S_st,2)/(2*bw))*1/.3048;
    figure(figcount);
    S_image = v; %edited to scale range to d^3/2
    S_image = fliplr(rot90(S_image));
    cr1 = -80; %(ft)
    cr2 = 80; %(ft)
    dr1 = 1 + Rs/.3048; %(ft)
    dr2 = 350 + Rs/.3048; %(ft)
    %data truncation
    dr_index1 = round((dr1/max_range)*size(S_image,1));
    dr_index2 = round((dr2/max_range)*size(S_image,1));
    cr_index1 = round(( (cr1+zpad*delta_x/(2*.3048)) /(zpad*delta_x/.3048))*size(S_image,2));
    cr_index2 = round(( (cr2+zpad*delta_x/(2*.3048))/(zpad*delta_x/.3048))*size(S_image,2));
    trunc_image = S_image(dr_index1:dr_index2,cr_index1:cr_index2);
    downrange = linspace(-1*dr1,-1*dr2, size(trunc_image,1)) + Rs/.3048;
    crossrange = linspace(cr1, cr2, size(trunc_image, 2));
    %scale down range columns by range^(3/2), delete to make like
    %dissertation again
    clear ii;
    for ii = 1:size(trunc_image,2)
        trunc_image(:,ii) = (trunc_image(:,ii)').*(abs(downrange*.3048)).^(3/2);
    end
    trunc_image = dbv(trunc_image); %added to scale to d^3/2
    imagesc(crossrange, downrange, trunc_image, [max(max(trunc_image))-40, max(max(trunc_image))-0]);  
    colormap('default');
    title('Final Image');
    ylabel('Downrange (ft)');
    xlabel('Crossrange (ft)');
    axis equal;
    cbar = colorbar;
    set(get(cbar, 'Title'), 'String', 'dB','fontsize',13); 
    print(gcf, '-djpeg100', 'final_image.jpg');
    if close_as_you_go == 1,
        close(figcount);
    end
    figcount = figcount + 1;

先做了fft，不知道是不是就是ppt里说的calibration matrix。

%along track FFT (in the slow time domain)
%first, symetrically cross range zero pad so that the radar can squint
zpad = 2048; %cross range symetrical zero pad
szeros = zeros(zpad, size(sif,2));
for ii = 1:size(sif,2)
    index = round((zpad - size(sif,1))/2);
    szeros(index+1:(index + size(sif,1)),ii) = sif(:,ii); %symetrical zero pad
end
sif = szeros;
clear ii index szeros;


S = fftshift(fft(sif, [], 1), 1);
clear sif;
Kx = linspace((-pi/delta_x), (pi/delta_x), (size(S,1)));

然后就是matched filter

%matched filter

%create the matched filter eq 10.8
for ii = 1:size(S,2) %step thru each time step row to find phi_if
    for jj = 1:size(S,1) %step through each cross range in the current time step row
        %phi_mf(jj,ii) = -Rs*Kr(ii) + Rs*sqrt((Kr(ii))^2 - (Kx(jj))^2);
        phi_mf(jj,ii) = Rs*sqrt((Kr(ii))^2 - (Kx(jj))^2);
        Krr(jj,ii) = Kr(ii); %generate 2d Kr for plotting purposes
        Kxx(jj,ii) = Kx(jj); %generate 2d Kx for plotting purposes
    end
end
smf = exp(j*phi_mf); %%%%%%%%%%%%

%note, we are in the Kx and Kr domain, thus our convention is S_mf(Kx,Kr)

%appsly matched filter to S
S_mf = S.*smf;

clear smf phi_mf;

然后是stolt interpolation也就是ppt里的stolt transform 

%perform the Stolt interpolation

%FOR DATA ANALYSIS
kstart =73;
kstop = 108.5;

%kstart = 95;
%kstop = 102;


Ky_even = linspace(kstart, kstop, 1024); %create evenly spaced Ky for interp for real data

clear Ky S_St;
%for ii = 1:size(Kx,2)
count = 0;
for ii = 1:zpad;
%for ii = round(.2*zpad):round((1-.2)*zpad)
    count = count + 1;
    Ky(count,:) = sqrt(Kr.^2 - Kx(ii)^2);
    %S_st(ii,:) = (interp1(Ky(ii,:), S_mf(ii,:), Ky_even)).*H;
    S_st(count,:) = (interp1(Ky(count,:), S_mf(ii,:), Ky_even));
end
S_st(find(isnan(S_st))) = 1E-30; %set all Nan values to 0
clear S_mf ii Ky;

最后是IFFT，也就是ppt上的2D IDFT

%perform the inverse FFT's
%new notation:  v(x,y), where x is crossrange
%first in the range dimmension
clear v Kr Krr Kxx Ky_even;
v = ifft2(S_st,(size(S_st,1)*4),(size(S_st,2)*4));

然后就可以绘图了 

bw = 3E8*(kstop-kstart)/(4*pi);
max_range = (3E8*size(S_st,2)/(2*bw))*1/.3048;
    figure(figcount);
    S_image = v; %edited to scale range to d^3/2
    S_image = fliplr(rot90(S_image));
    cr1 = -80; %(ft)
    cr2 = 80; %(ft)
    dr1 = 1 + Rs/.3048; %(ft)
    dr2 = 350 + Rs/.3048; %(ft)
    %data truncation
    dr_index1 = round((dr1/max_range)*size(S_image,1));
    dr_index2 = round((dr2/max_range)*size(S_image,1));
    cr_index1 = round(( (cr1+zpad*delta_x/(2*.3048)) /(zpad*delta_x/.3048))*size(S_image,2));
    cr_index2 = round(( (cr2+zpad*delta_x/(2*.3048))/(zpad*delta_x/.3048))*size(S_image,2));
    trunc_image = S_image(dr_index1:dr_index2,cr_index1:cr_index2);
    downrange = linspace(-1*dr1,-1*dr2, size(trunc_image,1)) + Rs/.3048;
    crossrange = linspace(cr1, cr2, size(trunc_image, 2));
    %scale down range columns by range^(3/2), delete to make like
    %dissertation again
    clear ii;
    for ii = 1:size(trunc_image,2)
        trunc_image(:,ii) = (trunc_image(:,ii)').*(abs(downrange*.3048)).^(3/2);
    end
    trunc_image = dbv(trunc_image); %added to scale to d^3/2
    imagesc(crossrange, downrange, trunc_image, [max(max(trunc_image))-40, max(max(trunc_image))-0]);  
    colormap('default');
    title('Final Image');
    ylabel('Downrange (ft)');
    xlabel('Crossrange (ft)');
    axis equal;
    cbar = colorbar;
    set(get(cbar, 'Title'), 'String', 'dB','fontsize',13); 
    print(gcf, '-djpeg100', 'final_image.jpg');
    if close_as_you_go == 1,
        close(figcount);
    end
    figcount = figcount + 1;
    

除了calibration matrix外，基本上ppt和matlab代码都能对上。

至于为啥按照ppt上的模块就能得到sar图像，以及具体matlab代码如何实现了那些模块的算法，还要后续仔细看论文才能知道。

 

这有个文章是对应ppt的中文翻译

http://www.360doc.com/content/14/1125/00/12979957_427822702.shtml