在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
