SAR成像学习（四）距离方向成像matlab代码解析 2

如果发射信号是线性调频信号，上一次讲的距离成像算法流程（匹配滤波方法）依然可以用，但那个流程要求T x =4X 0 c >T p  。如果T x <T p  ，即幅宽相对较小的情况，上一讲中的流程会带来一个问题，解决这个问题的办法是pulse compression。本文将会讨论这个puse compression的原理和实现。

1 what is pulse compression

对于线性调频信号：p(t)=a(t)exp(jβt+jαt 2 ) ，信号持续时间为T p  ，瞬时频率为β+aαt 带宽为：

±ω 0 =±αT p  

puse compression就是讲其对应的回波信号s(t) 与参考信号exp(jβt+jαt 2 ) 共轭相乘。

s(t)=∑ n σ n a(t−t n )exp(jβ(t−t n ))exp(jα(t−t n ) 2 ) 

其中t n =2x n /c 
puse compression的结果：

s c (t)=s ∗ (t)exp(jβt+jαt 2 )=∑ n σ n a ∗ (t−t n )exp(jβt n −jαt 2 n )exp(j2αtt n ) 

注意其中的频率成分2αt n  包含了目标的位置信息！这也正是puse compression能够进行成像处理的本后原因（代码在上一讲中）。
对上式进行傅里叶变换：

S c (ω) =∑ n σ n exp(jβt n −jαt 2 n −jωt n )psf ω (ω−2αt n )=∑ n σ n exp(jβt n −jαt 2 n −jωt n )psf ω (ω−4αx n /c)  

其中的exp(jβt n −jαt 2 n −jωt n ) 是相位函数。
其中psf ω  为频域的点扩散函数：

psf ω =F (t) [a ∗ (t)] 

很明显这是一个sinc函数。注意sinc函数峰值位置对应了（注意对应关系）了目标的位置。这个对应关系为：

ω=4αxc  

这就很好地解释了pulse compression方法进行成像的原理。现在我们来看一看s c (t) （时域）压缩信号要求的采样间隔。因为时域信号的范围是x∈[X c −X 0 ,X c +X 0 ] ，所以对应的有ω∈[4α(X c −X 0 )c ,4α(X c +X 0 )c ] 。现在我们找到了频域的带宽24αX 0 c =2αT x  ，其中T x =4X 0 c  。
那么时域压缩信号的采样间隔为：

Δ tc <=παT x   

注意上面的T x  对应成像的幅宽。

2 为什么要在匹配滤波成像过程中进行pulse compression

上一讲中说到利用匹配滤波进行成像处理的话，信号的采样间隔是Δ t <=πB 0  =παT p   。
可见两种成像方法（匹配滤波和时域压缩）所要求的采样间隔是不一样的。
注意硬件上希望采样间隔较大。
现在考虑一种情况：T x <T p  ，此时压缩信号所要求的采样间隔较大，如果我用时域压缩的方法进行成像处理，没有问题。但是假如我偏偏想用匹配滤波的方法进行成像处理，而且硬件上又必须要用Δ tc  进行采样，结果是什么？结果是回波信号亚采样，频谱混叠，成像失败！
有没有解决办法？有，就是在匹配滤波成像过程中进行pulse compression，利用没有混叠的压缩信号频谱得到没有混叠的回波信号，然后就是按照匹配滤波的流程进行成像。
也就是说因为亚采样（为了迁就硬件上的要求），我们不能直接利用回波信号，而是要用后处理得到的回波信号。
T x >T p  和T x <T p  两种情况下的匹配滤波成像处理流程别如下：

T x >T p   匹配滤波成像处理流程

T x <T p   匹配滤波成像处理流程

这两个过程在matlab代码中体现得很清楚了，不再展开。

3 为什么两种方法的处理结果不同呢

这是上一讲最后的一个疑问：case T x >T p   （由下面的代码其实上一讲也有代码）两种方法恢复出来的结果不一样！具体说是时域压缩方法结果是错误的。想了好久之后终于找到了答案。
时域压缩方法结果确实是错误的。因为时域压缩信号亚采样了，因为代码中用的是回波信号的采样间隔，而T x >T p   。
那么怎样才能得到正确的结果呢？书上说可以通过没有混叠的回波信号得到没有混叠的压缩信号。就像可以通过没有混叠的压缩信号得到没有混叠的回波信号一样。
这样绕来绕去的原因是：硬件上希望采样间隔较大。
所以我们可以总结：任何情况下T x >T p  或者T x <T p  ，如果必须满足硬件上的要求（采样间隔大一些），那么回波信号和压缩信号必然有一个是亚采样的，直接利用匹配滤波方法或者时域压缩方法的话，必然有一个会失败。避免失败的方法就是通过没有亚采样的那个信号间接得到本来亚采样的信号。
这个间接的方法说白了就是内插，升采样。

4 代码及注释

这里的代码包括了上一讲和这一讲的所有内容。
参考数据及代码来源

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%基本参数clear clccolormap(gray(256))cj=sqrt(-1);pi2=2*pi;%c=3e8;                   % Propagation speedB0=100e6;                % Baseband bandwidth is plus/minus B0w0=pi2*B0;fc=1e9;                  % Carrier frequencywc=pi2*fc;Xc=2.e3;                 % Range distance to center of target areaX0=50;                   % target area in range is within [Xc-X0,Xc+X0]%% Case 1: Tp=.1e-6;               % Chirp pulse duration% % % Case 2: Tp=10e-6;               % Chirp pulse durationalpha=w0/Tp;             % Chirp ratewcm=wc-alpha*Tp;         % Modified chirp carrier 即是式子中的beta%dx=c/(4*B0);             % Range resolution%采样间隔dt=pi/(2*alpha*Tp);      % Time domain sampling (gurad band plus minus                         % 50 per) or use dt=1/(2*B0) for a general                         % radar signal                         % 可能是因为信号是实信号 所以是采样频率是信号带宽的2倍Tx=(4*X0)/c;             % Range swath echo time period%压缩信号的采样间隔dtc=pi/(2*alpha*Tx);     % Time domain sampling for compressed signal                         % (gurad band plus minus 50 per)%Ts=(2*(Xc-X0))/c;        % Start time of samplingTf=(2*(Xc+X0))/c+Tp;     % End time of sampling% If Tx < Tp, choose compressed signal parameters for measurement% 因为硬件倾向于信号有着较小的带宽，所以如果压缩后的信号可以满足这个，首先会选择使用这个压缩后的信号% 可想而知，这样子回波信号（未压缩的）采样的间隔就过大了，所以需要进行升采样。flag=0;                  % flag=0 indicates that Tx > Tpif Tx < Tp, flag=1;                 % flag=1 indicates that Tx < TP dt_temp=dt;             % Store dt dt=dtc;                 % Choose dtc (dtc > dt) for data acquisition %就是说用dtc去采样回波信号（它应该是用dt采样才不会混叠）end;% Measurement parameters%n=2*ceil((.5*(Tf-Ts))/dt);   % Number of time samplest=Ts+(0:n-1)*dt;             % Time array for data acquisition%注意这里用了dt 后面以自变量的形式去求的回波信号序列也就是采样了%如果它本质上是dtc，那么就是亚采样了。%这三个参数也是取决于dt的 如果Tx < Tp后面会重新计算 因为这些参数是属于未压缩的回波信号的%如果Tx > Tp 这里的设置当然就是对的啦dw=pi2/(n*dt);               % Frequency domain samplingw=wc+dw*(-n/2:n/2-1);        % Frequency array (centered at carrier)x=Xc+.5*c*dt*(-n/2:n/2-1);   % range bins (array); reference signal is                             % for target at x=Xc.kx=(2*w)/c;                  % Spatial (range) frequency array%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ntarget=4;                        % number of targets%%%%%%%%%%%%% Targets' parameters  %%%%%%%%%%%%%%%%%%%% xn: range;               fn: reflectivity 发射系数%xn(1)=0;                   fn(1)=1;xn(2)=.7*X0;               fn(2)=.8;xn(3)=xn(2)+2*dx;          fn(3)=1.;xn(4)=-.5*X0;              fn(4)=.8;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SIMULATION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%s=zeros(1,n);              % Initialize echoed signal array%调节na可以对比不同情况下的匹配滤波和时间域压缩两种方法的成像结果的差异%结果是随着na的增大 时间域压缩方法变得不好na=0;                      % Number of harmonics in random phase               ar=rand(1,na);             % Amplitude of harmonicster=2*pi*rand(1,na);       % Phase of harmonicsfor i=1:ntarget;   td=t-2*(Xc+xn(i))/c;   pha=wcm*td+alpha*(td.^2);         % Chirp (LFM) phase   for j=1:na;                       % Loop for CPM harmonics 相位调制！出于ECCM的考虑 amplitude modulating an LFM chirp      pha=pha+ar(j)*cos((w(n/2+1+j)-wc)*td+ter(j));   end;   s=s+fn(i)*exp(cj*pha).*(td >= 0 & td <= Tp);%回波信号 采样后的end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 第一种成像方法 匹配滤波% If flag=1, i.e., Tx < Tp, perform upsmapling%if flag == 1, td0=t-2*(Xc+0)/c; pha0=wcm*td0+alpha*(td0.^2);        % Reference chirp phase scb=conj(s).*exp(cj*pha0);          % Baseband compressed signal                                     % (This is done by hardware) scb=[scb,scb(n:-1:1)];  % Append mirror image in time to reduce wrap                         % around errors in interpolation (upsampling)                         %DFT前进行了一个镜像扩展 fscb=fty(scb);          % F.T. of compressed signal w.r.t. time% dt=dt_temp;                     % Time sampling for echoed signal 注意这个是理论上的回波信号采样间隔 之前存下的 n_up=2*ceil((.5*(Tf-Ts))/dt);   % Number of time samples for upsampling %理论上的采样个数 nz=n_up-n;                      % number of zeros for upsmapling is 2*nz %缺少的采样个数 fscb=(n_up/n)*[zeros(1,nz),fscb,zeros(1,nz)]; %注意镜像扩展 系数是？% scb=ifty(fscb); scb=scb(1:n_up);            % Remove mirror image in time%% Upsampled parameters%未压缩的回波信号的一些参数 n=n_up; t=Ts+(0:n-1)*dt;             % Time array for data acquisition dw=pi2/(n*dt);               % Frequency domain sampling w=wc+dw*(-n/2:n/2-1);        % Frequency array (centered at carrier) x=Xc+.5*c*dt*(-n/2:n/2-1);   % range bins (array); reference signal is                              % for target at x=Xc. kx=(2*w)/c;                  % Spatial (range) frequency array% td0=t-2*(Xc+0)/c; s=conj(scb).*exp(cj*wcm*td0+cj*alpha*(td0.^2));  % Decompressionend;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Reference echoed signal%td0=t-2*(Xc+0)/c;pha0=wcm*td0+alpha*(td0.^2);         % Chirp (LFM) phasefor j=1:na;                          % Loop for CPM harmonics   pha0=pha0+ar(j)*cos((w(n/2+1+j)-wc)*td0+ter(j));end;s0=exp(cj*pha0).*(td0 >= 0 & td0 <= Tp);%% Baseband conversion%sb=s.*exp(-cj*wc*t);sb0=s0.*exp(-cj*wc*t);%figureplot(t,real(sb))xlabel('Time, sec')ylabel('Real Part')title('Baseband Echoed Signal')axis('square')axis([Ts Tf 1.1*min(real(sb)) 1.1*max(real(sb))])print P1.1a.pspause(1)%figureplot(t,real(sb0))xlabel('Time, sec')ylabel('Real Part')title('Baseband Reference Echoed Signal')axis('square')axis([Ts Tf 1.1*min(real(sb0)) 1.1*max(real(sb0))])print P1.2a.pspause(1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RECONSTRUCTION%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fourier Transform%fsb=fty(sb);fsb0=fty(sb0);% Power equalization%mag=abs(fsb0);amp_max=1/sqrt(2);     % Maximum amplitude for equalizationafsb0=abs(fsb0);P_max=max(afsb0);I=find(afsb0 > amp_max*P_max);nI=length(I);fsb0(I)=((amp_max*(P_max^2)*ones(1,nI))./afsb0(I)).* ...        exp(cj*angle(fsb0(I)));%% Apply a window (e.g., power window) on fsb0 here%E=sum(mag.*abs(fsb0));%figureplot((w-wc)/pi2,abs(fsb)) %注意自变量的范围 单位xlabel('Frequency, Hertz')ylabel('Magnitude')title('Baseband Echoed Signal Spectrum')axis('square')print P1.3a.pspause(1)%figureplot((w-wc)/pi2,abs(fsb0))xlabel('Frequency, Hertz')ylabel('Magnitude')title('Baseband Reference Echoed Signal Spectrum')axis('square')print P1.4a.pspause(1)%% Matched Filtering%fsmb=fsb.*conj(fsb0);%figureplot((w-wc)/pi2,abs(fsmb))xlabel('Frequency, Hertz')ylabel('Magnitude')title('Baseband Matched Filtered Signal Spectrum')axis('square')print P1.5a.pspause(1)%% Inverse Fourier Transform%smb=ifty(fsmb);       % Matched filtered signal (range reconstruction)%% Display%figureplot(x,(n/E)*abs(smb)) %注意自变量的对应关系xlabel('Range, meters')ylabel('Magnitude')title('Range Reconstruction Via Matched Filtering')axis([Xc-X0 Xc+X0 0 1.1]); axis('square')print P1.6a.pspause(1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Time domain Compression 第二种成像方法% 当Tp<Tx (case 1)总是得不到想要的形式？这种方法恢复出来的结果看起来也不对% 这是因为 这个时候压缩信号的采样间隔太大了 需要通果不混叠的未压缩信号恢复出不混叠的压缩信号 在进行成像处理！% 当然如果Tp>Tx (case 2) 两种方法恢复出来的结果基本是一致的。%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%td0=t-2*(Xc+0)/c;scb=conj(s).* ...     exp(cj*wcm*td0+cj*alpha*(td0.^2));  % Baseband compressed signal%figuresubplot(1,2,1)plot(t,real(s))subplot(1,2,2)plot(t,real(scb))xlabel('Time, sec')ylabel('Real Part')title('Time Domain Compressed Signal')axis('square')print P1.7a.pspause(1)fscb=fty(scb);X=(c*(w-wc))/(4*alpha);     % Range array for time domain compression 注意这个转换关系figuresubplot(1,2,1)plot(w,(dt/Tp)*abs(fscb))subplot(1,2,2)plot(X+Xc,(dt/Tp)*abs(fscb))xlabel('Range, meters')ylabel('Magnitude')title('Range Reconstruction Via Time Domain Compression')axis([Xc-X0 Xc+X0 0 1.1]); axis('square')%axis([Xc-X0 Xc+X0 0 1.1*max(abs(fscb))]); axis('square')print P1.8a.pspause(1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%