自制合成孔径雷达(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 profilesfor ii = Nrp+1:size(rpstart,1)-Nrpif rpstart(ii) == 1 & sum(rpstart(ii-Nrp:ii-1)) == 0count = 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;jjfor 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 == 1count = count + 1;SIF = RP(jj,ii:ii+N-1)' + SIF;endend%hilbert transformq = 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 datafor 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)-Nrpif rpstart(ii) == 1 & sum(rpstart(ii-Nrp:ii-1)) == 0count = 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 == 1count = count + 1;SIF = RP(jj,ii:ii+N-1)' + SIF;endend

然后做希尔伯特变换

    %hilbert transformq = 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:NH(ii) = 0.5 + 0.5*cos(2*pi*(ii-N/2)/N);
endfor 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);endfigcount = 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);endfigcount = figcount + 1;
endif 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);endfigcount = figcount + 1;
endif 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);endfigcount = 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_iffor 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 purposesKxx(jj,ii) = Kx(jj); %generate 2d Kx for plotting purposesend
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);endfigcount = figcount + 1;
endclear 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);endfigcount = 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 dataclear 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);endfigcount = figcount + 1;
end%apply hanning window to data, cleans up data ALOT
N = size(Ky_even,2);
for ii = 1:NH(ii) = 0.5 + 0.5*cos(2*pi*(ii-N/2)/N);
endfor 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/2S_image = fliplr(rot90(S_image));cr1 = -80; %(ft)cr2 = 80; %(ft)dr1 = 1 + Rs/.3048; %(ft)dr2 = 350 + Rs/.3048; %(ft)%data truncationdr_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 againclear ii;for ii = 1:size(trunc_image,2)trunc_image(:,ii) = (trunc_image(:,ii)').*(abs(downrange*.3048)).^(3/2);endtrunc_image = dbv(trunc_image); %added to scale to d^3/2imagesc(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);endfigcount = 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_iffor 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 purposesKxx(jj,ii) = Kx(jj); %generate 2d Kx for plotting purposesend
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 dataclear 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/2S_image = fliplr(rot90(S_image));cr1 = -80; %(ft)cr2 = 80; %(ft)dr1 = 1 + Rs/.3048; %(ft)dr2 = 350 + Rs/.3048; %(ft)%data truncationdr_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 againclear ii;for ii = 1:size(trunc_image,2)trunc_image(:,ii) = (trunc_image(:,ii)').*(abs(downrange*.3048)).^(3/2);endtrunc_image = dbv(trunc_image); %added to scale to d^3/2imagesc(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);endfigcount = figcount + 1;

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

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

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

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

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

自制合成孔径雷达(3) doppler代码解读
上一篇帖子,看完了基于SDR的多普勒雷达,就可以看看硬件雷达的多普勒测速的DSP代码了. 先看一下这个图: 我们需要的多普勒频移的测量结果是从混频器(Multiply Conjugate)输出的,也就 ...
SAR成像系列：【8】合成孔径雷达（SAR）成像算法-压缩感知（Compressed Sensing，CS）成像算法（附Matlab代码）
压缩感知(Compressed Sensing,CS)该理论指出:对于满足约束等距条件(Restricted Isometry Property,RIP)的稀疏或可压缩信号,通过低于(甚至远低于)Ny ...
SAR成像系列：【9】合成孔径雷达（SAR）成像算法-波数域（omega-K）成像算法[也叫距离徙动（RM）算法]（附Matlab代码）
波数域()成像算法作为本系列的最后一种成像算法介绍.关于SAR成像的其他的各种改进算法就不一一列举了.在实际成像中,万变不离其踪,最主要的是关注成像的几何模型,再根据指标选择不同的基础成像算法,然后进 ...
SAR成像系列：【1】合成孔径雷达（SAR）成像概述
本系列主要介绍合成孔径雷达(SAR)成像的关键技术,帮助入门者更好的理解雷达成像原理及算法. (1)雷达原理雷达的英文式 Radar ,源于 Radio Detection and Ranging ...
SAR成像系列：【4】合成孔径雷达（SAR）距离历程分析
(1)SAR的距离历程-正侧视合成孔径雷达(SAR)成像中,所有算法的基础是对距离历程(Range History)的分析."合成孔径"从表面上看是合成了雷达阵列,从信号上看是获 ...
SAR成像系列：【5】合成孔径雷达（SAR）成像算法-距离多普勒（RD）算法（附Matlab代码）
完整的距离多普勒算法主要包括距离压缩.距离徙动矫正(矫正距离走动和距离弯曲).方位压缩等步骤.其中距离走动矫正即可在时域进行也可在频域进行,而距离弯曲矫正一般在多普勒域进行.在距离多普勒域叫作RCMC ...
SAR成像系列：【3】合成孔径雷达（SAR）的二维回波信号与简单距离多普勒（RD）算法（附matlab代码）
合成孔径雷达发射信号以线性调频信号(LFM)为基础,目前大部分合成孔径雷达都是LFM体制,为了减轻雷达重量也采用线性调频连续波(FMCW)体制:为了获得大带宽亦采用线性调频步进频(FMSF)体制. ( ...
【合成孔径雷达】SAR 距离多普勒算法（RD）+ Chirp Scaling算法（CS）+ InSAR干涉SAR人造场景仿真和干涉处理 + PolSAR极化定标【MATLAB】
最近发现在CSDN等很多平台有抄袭本文内容,或本人GitHub源码的行为. 在此特别声明如下:欢迎转载,欢迎分享,这是对我工作的最大肯定.但是请务必标明出处和来源.严禁任何形式的.无说明的.或抄袭成个 ...
什么是合成孔径雷达（SAR）
"合成孔径雷达",什么是"合成",什么是"孔径",什么是"雷达"? 合成孔径雷达(Synthetic Aperture ...

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

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

最新文章

热门文章