伪随机编码冲雷达信号处理仿真
目录
1.题目论述
仿真伪随机相位编码脉冲雷达的信号处理。设码频为各学生学号末两位数（22），单位为MHz，伪码周期内码长为127，占空比10%，雷达载频为10GHz，输入噪声为高斯白噪声。
目标模拟分单目标和双目标两种情况，目标回波输入信噪比可变（-35dB～10dB），目标速度可变（0～1000m/s），目标幅度可变（1～100），目标距离可变（0～10000m），相干积累总时宽不大于10ms。
单目标时，给出回波视频表达式；脉压和FFT 后的表达式；仿真m序列的双值电平循环自相关函数，给出脉压后和FFT 后的输出图形；通过仿真说明各级处理的增益，与各级时宽和带宽的关系；仿真说明脉压时多卜勒敏感现象和多卜勒容限及其性能损失（脉压主旁比与多卜勒的曲线）。
双目标时，仿真出大目标旁瓣盖掩盖小目标的情况；仿真出距离分辨和速度分辨的情况。
2.回波信号及脉压-FFT后的表达式
2.1视频表达式：
伪随机二相编码脉冲雷达的发射信号复数形式可以表示为：
\[{s_s}(t) = c(t)\exp (j2\pi {f_0}t)\] [1]
其中为载频，为相位调制因子，即m序列的二值函数，取1或-1.
则接受到的回波信号可写为：
\[{s_r}(t) = Ac(t - \tau )\exp [j2\pi {f_0}(t - \tau ) + j2\pi {f_d}t]\] [2]
下混频后得到：
\[s(t) = Ac(t - \tau )\exp (j2\pi {f_d}t - j2\pi {f_0}\tau )\][3]
其中，
$$\tau = 2d/c = d \times 6.667e - 9$$ [4]
$${f_d} = 2{V_d}/\lambda = {v_d} \times \frac{{200}}{3}$$ [5]
2.2脉压后的表达式
信号处理从回波信号的脉压开始。为之前所得混频后的回波信号，则脉压可以通过匹配滤波来实现。脉压后的表达式为：
\[{s_{MF}}(k) = h(t) \otimes s(t) = \int\limits_{ - \infty }^{ + \infty } {h(t - \tau )} \cdot s(\tau )d\tau \][6]
其中，是匹配滤波器的传递函数，仿真时用一个周期的m序列的反褶作为脉压器。
2.3 FFT后的表达式
在有限相干时间内经过相干积累得到的信号称为一帧。由于采用数字信号的处理方法，首先需要采样，也即AD变换。设定采样频率为，相干时间等于N个m序列周期，序列码长为p，则脉压后的样点数为N*p个，记为.
由于发射的是10%占空比的脉冲信号，即每10个序列周期发射一个完整序列。脉压后的每隔10p个样点就会出现一个峰。
按照每个距离门i对每隔10p进行采样，就可以得到INT(N/10)组数据(INT为向下取值)。第i组数据的表达式为
\[{s_{RC}}[i,n] = {s_{MF}}[10pn + i],n = 0,1,2, \cdot \cdot \cdot \] [7]
实际上也就是对多普勒信号以频率采样后的信号。对作FFT就可以提取到多普勒频率，根据式[5]就可以转换为目标速度.此外，对多普勒信号波形的采样频率也决定了雷达的测速范围。这将在后面再次提及。
至此，可以写出作NFFT点数FFT后的表达式：\[\begin{array}{l}
X[i,k] = FFT({s_{RC}}[i,n]) = \sum\limits_{n = 0}^{NFFT - 1} {{s_{RC}}[i,n]\exp ( - j2\pi kn/NFFT)} \\
{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{n = 0}^{NFFT - 1} {{s_{MF}}(10pn + i)\exp ( - j2\pi kn/NFFT)}
\end{array}\]
[8]
3 m序列的性质与仿真
m序列有多个性质，包括均衡特性、游程分布和自相关特性等。
均衡性是指：m序列中，1的个数比-1的个数多一个，所以.\[\sum\limits_{i = 0}^{p - 1} {c[i]} = - 1\]
自相关性可以通过自相关函数来体现。
周期自相关函数：由周期函数的自相关函数定义：
$$R(\tau ) = \frac{1}{{{T_0}}}\int_{ - \frac{{{T_o}}}{2}}^{\frac{{{T_o}}}{2}} {s(t)s(t + \tau )dt} $$
从而可知m序列的周期自相关函数为：
\[R(\tau ) = \left\{ {\begin{array}{*{20}{c}}
{1 - \frac{{p + 1}}{{{T_0}}}}&{0 \le \left| {\tau - i{T_0}} \right| \le \frac{{{T_0}}}{p},}&{i = 0,1,2, \cdot \cdot \cdot }\\
{ - \frac{1}{p}}&{Others}&{}
\end{array}} \right.\]
需要注意，上式中可以在连续域上取值。若取m码的码元离散值，自相关函数为：
\[R(x) = \left\{ {\begin{array}{*{20}{c}}
1&{x = 0}&{}\\
{ - \frac{1}{p}}&{x = 1,2,3, \cdot \cdot \cdot ,p - 1}&{}
\end{array}} \right.\] [9]
利用MATLAB仿真时，使用idinput()函数可以直接生成m序列。计算其循环自相关函数，仿真的结果如下图，可见与理论式[7]相符。
Fig. 1 m序列的循环自相关函数
非周期自相关函数：根据自相关函数表达式\({R_x}(\tau ) = E[X(t)X(t + \tau )]\)，对单个周期m序列求自相关，结果如下图：
Fig. 2 127点m序列非周期自相关函数
非周期自相关函数与循环自相关函数一样，可以看到在$\tau = 0$处有最大自相关性，但非周期自相关函数在$\tau \ne 0$处的自相关性相比于周期自相关函数，更加杂乱。m序列调相连续波雷达工作利用了m序列的周期自相关性，而脉冲雷达脉冲内发射单个m序列调相的波形，因而工作在非周期自相关特性下。
4 噪声与回波信号
雷达系统接收到的回波信号受信道特性影响，含有加性高斯白噪声。高斯白噪声指的是：1.噪声幅值符合高斯分布；2.噪声的功率谱密度为常数。
4.1 噪声的生成
在MATLAB中，可以利用randn()函数产生标准差为1，均值为0的高斯信号，但并不能保证产生的信号具有均匀的频谱密度。wgn(m,n,p) 可以产生一个m行n列的高斯白噪声的矩阵，噪声功率是p dBW.而awgn(x,SNR)可以直接在信号x中加入信噪比为SNR的高斯白噪声。
接下来验证这三个函数产生噪声的性质。
高斯分布特性可以通过统计幅值分布，作概率密度直方图来简单估计。根据经验，统计的区间数和样本数之间有最优关系：
$$K = 1.87{(N - 1)^{\frac{2}{5}}}$$
[10]
若有1000点噪声序列，根据公式，应划分29个区间。
噪声的"白"特性可以由功率谱密度来验证，功率谱密度和自相关函数为傅立叶变换对关系，因而首先求得自相关函数，作傅立叶变换可得功率谱密度。
下图给出了三个函数产生的噪声序列，绘出了概率分布直方图；并分别给出了时域波形、自相关函数和自相关函数的傅立叶变换——功率谱密度。
Fig. 3 randn(),wgn(),awgn()产生序列的直方图统计
Fig. 4 randn(),wgn(),awgn()产生序列的自相关函数和功率谱密度
其中randn函数和wgn函数直接产生1000点噪声序列，而awgn是对1000点幅值为1的信号添加5dB信噪比的噪声，作图时减去了直流分量1，得到均值为0的噪声再作自相关和功率谱密度。
结果可见，三种噪声序列幅值都符合高斯分布特性。虽然randn函数的描述中并未提及功率谱密度均匀，但三种函数产生的噪声都属于高斯白噪声，即幅值呈高斯分布，功率谱密度为常数。
高斯过程通过线性系统仍旧为高斯过程。所以对于含高斯白噪声的回波信号，混频到视频后仍为高斯白噪声。仿真时可以在视频信号中加入高斯白噪声，模拟实际系统的噪声性能。
式[3]的回波信号中加入特定信噪比的噪声，与发射信号对比如下图。
|
|
|
|
Fig. 5 不同信噪比的回波信号与发射信号对比
上图给出了发射的脉冲基带信号(左上图)，模拟目标距离2000m和速度 20m/s情况下的回波。回波信号中可见时延和多普勒频率调制。当信噪比逐渐减小时，噪声幅度增大。图中给出了信噪比为50dB，10dB，0dB，-10dB，-30dB情况下的回波波形。当信噪比达到-30dB时，有用信号已经完全淹没在噪声中。
4.2验证awgn()函数输出结果的信噪比
awgn(x,snr,'measured')函数据帮助文件描述，首先会计算信号序列x的功率，然后根据所给的信噪比SNR调整高斯白噪声的幅值，加到信号上。
下图是利用awgn()函数对发射脉冲信号增加20dB噪声后的结果。
Fig. 6 利用awgn函数对发射信号增加20dB噪声
其中标识的区域1是纯噪声，而区域2是信号叠加了噪声。
记原信号为$x[i]$，添加了噪声后信号为$x'[i]$。20dB信噪比是指一个周期内信号与噪声的平均功率之比是20dB。由于发射信号$x[i]$是127点长度的脉冲，占空比10%，即周期为1270，且${x^2}[i] = 1$。计算一个周期内的平均信号功率为：
$${S_{av}} = \frac{{\sum\limits_{i = 0}^{i = 126} {{x^2}[i]} }}{{1270}} = 0.1$$. [11]
噪声的平均功率可以通过计算$x'[i]$序列中区域1中样点的能量，除以点数，得到平均功率。即噪声平均功率：
$${N_{av}} = \frac{{\sum\limits_{i = 127}^{i = 1269} {x{'^2}[i]} }}{{127*9}}$$. [12]
这一项无法直接写出结果，但可以利用MATLAB对数据进行统计计算，得到${N_{av}}=9.4182e-4$.
所以信噪比$SNR = \frac{{{S_{av}}}}{{{N_{av}}}} = \frac{{0.1}}{{9.4182e - 4}} = 106.177 = 20.26dB$.符合设定参数
生成-5dB信噪比的信号再次进行验证，给出输入输出信号对比图如下：
Fig. 7利用awgn函数对发射信号增加-5dB噪声
根据式[11]、[12]，${S_{av}} = 0.1$不变，统计${N_{av}} = \frac{{\sum\limits_{i = 127}^{i = 1269} {x{'^2}[i]} }}{{127*9}}=0.3198$.得：$$SNR = \frac{{{S_{av}}}}{{{N_{av}}}} = \frac{{0.1}}{{0.3198}} = 0.3127 = - 5.045dB$$，
符合设定的参数。
4.3小结
信号处理的仿真从回波正式开始，目的就是从噪声中检测出回波所含的延时信息和多普勒信息。极低的信噪比给信号处理带来挑战，同时也将体现现代信号处理方法的强大能力。
5 各级处理增益和时宽带宽的仿真
对回波信号的处理典型方法是首先通过脉压，提取距离信息。依据不同发射信号类型，此步骤可以得到不同的增益。然后，按照距离门排列，作FFT，得到多普勒频率信息，进而确定目标速度。FFT又可以得到一定的增益。通过这两个步骤，可以发掘出深埋于噪声下的回波信号所含的信息。
5.1 脉压增益
理论增益：根据式[6]，以127点m序列作为本地匹配信号，对回波进行相关积累，则脉压的理论增益$D = BT$，其中是B带宽，对于二相编码序列，\(B = \frac{1}{\tau }\)；D为时宽。所以\(D = \frac{1}{\tau }{\kern 1pt} {\kern 1pt} p\tau = p\). 仿真时使用127点m序列，因而理论压缩比为127(21dB)，并且可知，只与m序列长度p有关。
仿真：
1. 10dB信噪比回波情况（大信噪比）对目标速度为0，距离为0，信噪比为10dB的回波按照式[6]所述进行脉压处理，得到序列，并去除所得结果的前127个点，使距离对齐，得到${s_{MF}}[k]$。则在第0个样点处将会有脉压产生的峰，这点也是信噪比最高的点。回波信号与脉压后的结果如下图所示：
Fig. 8 10dB信噪比回波脉压后的图像
为了便于观察，对脉压信号取对数坐标，结果如下图：
Fig. 9 10dB信噪比（20dB瞬时信噪比）脉压后取对数
从图上粗略判断，可见脉压后的信噪比大约在40dB左右，即增益为20dB左右，接下来精确计算脉压后的信噪比。
脉冲压缩是指将长度127点的脉冲序列经过匹配滤波等方法，转换成了时宽很窄的峰，而幅值增大。脉压的增益D就是脉压后的最大信号功率与噪声功率之比相对于脉压前瞬时信噪比的增益。
记脉压前信号为$x'[i]$，平均信噪比为$SNR{I_{av}}$（也是awgn函数参数设置的信噪比），脉压输出信号为${s_{MF}}[k]$。则由于输入信号是10%占空比的脉冲，从平均信噪比换算到信号噪声的瞬时功率比，输入信号瞬时信噪比\(SNR{I_{ins}}\)是\(SNR{I_{av}}\)的十倍，即
$$SNR{I_{ins}} = 10SNR{I_{av}}$$[13]
脉压后的信号功率即k=0时的峰值高度，而噪声功率可以取fig,7所框出的纯噪声信号区域的平均功率。
所以脉压后峰值信噪比可以定义为:$$SNR{O_{ins}} = \frac{{{s_{MF}}^2[0]}}{{\frac{{\sum\nolimits_{200}^{799} {{s_{MF}}^2[k]} }}{{600}}}}$$
[14]
当输入信号平均信噪比为10dB时，根据式[13], $SNR{I_{ins}}=20dB$。脉压后的信噪比根据式[14]，在MATLAB中统计噪声平均功率，得到
\(SNR{O_{ins}}=125.8^2/1.1829=13378.68=41.26dB\)，故增益为
D=\(SNR{O_{ins}}\) / \(SNR{O_{av}}\)=41.26dB-20dB=21.26dB，与理论的21dB相差1.2%.
2. -15dB信噪比回波情况（小信噪比）改变回波信号的平均信噪比\(SNR{O_{av}}\)到-15dB，再次验证脉压后的增益。按照之前同样的方法，作回波与脉压后的图，并对脉压信号取对数，结果如下：
Fig. 10 -15dB信噪比回波脉压后的图像
Fig. 11 -15dB信噪比（-5dB瞬时信噪比）脉压后取对数
127点m序列应该得到21dB增益，回波瞬时信噪比为-5dB，因而脉压后信噪比应为21+(-5)=16dB。从fig.10来看，信噪比确实在16dB左右。
按照式[14]，精确计算脉压后的信噪比，得
\(SNR{O_{ins}}\)=(116.5)^2/375.7187=36.1234=15.58dB.
增益D=\(SNR{O_{ins}}\)/ \(SNR{O_{av}}\)=15.58dB-（-5）dB=20.58dB，与理论的21dB相差2%.
从以上两例可见，仿真结果都接近理论值。小信噪比时仿真得到的增益误差略大于大信噪比时情况，原因在于统计的信号其实是叠加了噪声的信号，大信噪比时噪声对信号幅值的影响较小，因而误差较小。
5.2脉压后的时宽和带宽
在不加入噪声和多普勒频移的情况下，对回波进行脉压，考察时宽和带宽的变化。由于回波信号的采样频率等于m序列码频，脉压输出只有一个样点宽度的峰值。因而需要对脉压后的信号作内插，在每两个样点间插入20个差值样点，结果如下，脉压后的时宽变为1个子码的时宽，是原来的1/127.
Fig. 12 脉压信号的时宽
对脉压信号作FFT，变换到频域，得到下图，可见脉压之后信号的带宽仍为22MHz.
Fig. 13 脉压信号的带宽
5.3 FFT的增益
5.1已经验证，127点m序列脉冲雷达的回波信号经过脉压处理，可以得到21dB的增益。之前仿真时的条件是距离、速度都为0。当目标距离不为零时，脉压的峰值会延迟出现，延迟时间就是回波延迟时间，从而可以测得距离信息；当目标速度不为0时，回波将会被多普勒频率调制。此时经过脉压，得到的窄脉冲峰值也会含有多普勒频率的幅度包络。Fig.4中已经展示过这种情况。对这一振幅包络作FFT就可以得到多普勒频率，进而得到速度信息。
在作FFT之前，需要先对脉压信号按照距离门进行排列，便于最后速度和距离的同时显示。题目中规定相干时间不大于10ms，因而最多对173个脉冲周期积累，每个周期内有1270个样点。把每个样点作为一个距离门，一共有1270个距离门；为了数据对齐，每个距离门取171个样作FFT。
对171个点作FFT，理论将产生增益$$G = 10\log (171) = 22.3dB$$
FFT输入信噪比：仿真时取距离1000m，速度2m/s，回波信噪比0dB。将脉压后的信号按照距离门排列，得到171*1270维矩阵C[m,n]。结果如下。
Fig. 14 按距离门排列后的脉压图像和取对数之后
在图上找出最大幅值处的信号，位于第148个距离门的第88个样点上，最大信号幅度为128.8. 按照式[14]的形式在矩阵第88行上计算出噪声平均功率，算出信噪比。此处为了避开距离旁瓣，只统计噪声功率，样点从300取到1199，共900个点，再取平均。数学式为：
$$SNRi = \frac{{{C^2}{{[m,n]}_{\max }}}}{{{{\sum\nolimits_{k = 300}^{1199} {{C^2}[88,k]} } \mathord{\left/
{\vphantom {{\sum\nolimits_{k = 300}^{1199} {{C^2}[88,k]} } {900}}} \right.
\kern-\nulldelimiterspace} {900}}}}$$[15]
计算得$$SNRi = {(128.8)^2}/6.8243 = 2430.94 = 33dB$$.
FFT输出信噪比: 对C[m,n]的每行C[1:end，n]作256点FFT，得到256*1270维的矩阵FFT[j,k].作图如下：
Fig. 15 作FFT之后和对数坐标显示
从图中看到，最大信号幅值达到了2.169e4, 对应的位置在第148个距离门，第3个FFT样点处。按照同样的方法统计噪声功率，计算信噪比：
$$SNR0 = \frac{{FF{T^2}{{[j,k]}_{\max }}}}{{{{\sum\nolimits_{k = 300}^{1199} {{C^2}[3,k]} } \mathord{\left/
{\vphantom {{\sum\nolimits_{k = 300}^{1199} {{C^2}[3,k]} } {900}}} \right.
\kern-\nulldelimiterspace} {900}}}} = {(2.169e4)^2}/1.12e3 = 420050 = 56dB$$
所以FFT环节的增益G=\({{SNR0} \mathord{\left/
{\vphantom {{SNR0} {SNRi}}} \right.
\kern-\nulldelimiterspace} {SNRi}}\)=56dB-33dB=23dB.与理论增益22.3dB相比，误差3%.
5.4 FFT处理后的带宽
在无噪声情况下，设置目标参数距离为0m，速度10m/s, 脉压并作FFT后从二维矩阵中分离出含有FFT峰值的向量，取对数后得到下图。从图中读出FFT后的-3dB带宽为102Hz.
Fig. 16 FFT后的带宽
6 多普勒敏感现象与多普勒容限
6.1 多普勒敏感现象
理论分析：
式[3]表明了当目标具有径向速度时，反射的回波频率会加上多普勒频移，经过下混频去除载频后，残留的多普勒频移就表现为对基带信号进行了余弦包络的调幅，余弦波的频率就是多普勒频率${f_d}$。
当${f_d}$为0时，匹配滤波后得到最大增益；当${f_d}$增大，多普勒周期可以和随机序列周期相比拟时，被多普勒调幅后的回波信号与本地相干积累信号失配，主瓣降低，旁瓣升高。
设多普勒频率为${f_d}$，码元宽度为$\tau $，相关器长度为p，采样率为$\frac{1}{\tau }$，可从理论推导出多普勒敏感引起的脉压增益损失为：
$${L_p} = \left| {\frac{{\sin (\pi {f_d}p\tau )}}{{p\sin (\pi {f_d}\tau )}}} \right|$$ [16]
在本次仿真中$p = 127$，\(\tau = \frac{1}{{22e6}} = 4.545e - 8\).
根据表达式[16]绘出理论增益损失与多普勒频率、目标速度的关系曲线如下：
Fig. 17 理论增益损失与多普勒频率、目标速度的关系曲线
图中还标定了脉压增益下降-3dB时的多普勒频率为76.7KHz，对应速度1152m/s.
之前已经提到，以为间隔对混频后的信号进行采样率为的采样，按照距离门再次采样后得到多普勒波形的采样，采样率降低十倍，即
$${f_{s2}} = {{{f_{s2}}} \mathord{\left/
{\vphantom {{{f_{s2}}} {10p}}} \right.
\kern-\nulldelimiterspace} {10p}} = \frac{1}{{10p\tau }}$$ [17]
根据奈奎斯特采样定理，若要通过FFT正确分离出多普勒频率，采样率和多普勒频率之间需要满足：$${f_{s2}} \geqslant 2{f_d}$$.根据式[9], 得fd<=8661, 再由式[5]换算到速度，得到该雷达测速范围${V_d} \leqslant 129.86$m/s.
所以，尽管相位编码雷达对多普勒敏感，但由于此次仿真题目数据设计欠妥，不能在测速范围内仿真出多普勒敏感现象。倘若忽略超过测速范围而产生的盲速效应，依旧能够看到多普勒频移带来的增益损失现象，解释相位编码雷达多普勒容限低的特点。
仿真：
在大信噪比(100dB)下，设置目标距离d=0m，速度v可变。在速度v=0时得到最大脉压峰值，而峰值高度随着多普勒频率的增加而降低。下图展示了在不同目标速度情况下的回波信号和脉压峰值高度。
Fig. 18 不同目标速度下的回波信号和脉压峰值
从图中可以看到，当多普勒频率增大时，回波开始失配，脉压后的信号高度减小；当多普勒周期正好是m序列周期时，完全失配，脉压后峰值最小。
脉压的增益D等于脉压后的瞬时信噪比与回波瞬时信噪比的比值。回波瞬时信噪比由式[13]已经证明：$$SNR{I_{ins}} = 10SNR{I_{av}}$$；脉压后的信噪比可以定义为最大信号功率与平均噪声功率之比，即$$SNR{O_{ins}} = \frac{{{s_{\left. {MF} \right|\max }}^2}}{{{N_{av}}}}$$，所以增益D=$SNR{O_{ins}}$/$SNR{O_{av}}$. 仿真时只改变多普勒频率，因而$SNR{O_{ins}}$不变，平均噪声功率N_av不变。在fd=0时得到最大脉压增益D0，并以此归一化其他多普勒频率下的增益，表示为.
$$\eqalign{
& D({f_d}) = {{\frac{{SNR{O_{ins}}({f_d})}}{{SNR{I_{ins}}({f_d})}}} \mathord{\left/
{\vphantom {{\frac{{SNR{O_{ins}}({f_d})}}{{SNR{I_{ins}}({f_d})}}} {{D_0}}}} \right.
\kern-\nulldelimiterspace} {{D_0}}} = \frac{{{{\frac{{{s_{\left. {MF} \right|\max }}{{({f_d})}^2}}}{{{N_{av}}}}} \mathord{\left/
{\vphantom {{\frac{{{s_{\left. {MF} \right|\max }}{{({f_d})}^2}}}{{{N_{av}}}}} {10SNR{I_{av}}}}} \right.
\kern-\nulldelimiterspace} {10SNR{I_{av}}}}}}{{{D_0}}} = \frac{{{{\frac{{{s_{\left. {MF} \right|\max }}{{({f_d})}^2}}}{{{N_{av}}}}} \mathord{\left/
{\vphantom {{\frac{{{s_{\left. {MF} \right|\max }}{{({f_d})}^2}}}{{{N_{av}}}}} {10SNR{I_{av}}}}} \right.
\kern-\nulldelimiterspace} {10SNR{I_{av}}}}}}{{{{\frac{{{s_{\left. {MF} \right|\max }}{{(0)}^2}}}{{{N_{av}}}}} \mathord{\left/
{\vphantom {{\frac{{{s_{\left. {MF} \right|\max }}{{(0)}^2}}}{{{N_{av}}}}} {10SNR{I_{av}}}}} \right.
\kern-\nulldelimiterspace} {10SNR{I_{av}}}}}} \cr
& {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{{{s_{\left. {MF} \right|\max }}{{({f_d})}^2}}}{{{s_{\left. {MF} \right|\max }}{{(0)}^2}}} \cr} $$[18]
上式表明，计算多普勒敏感时的脉压增益损失比时，不需要统计基底噪声，归一化后的脉压增益即是脉压信号功率与无多普勒时脉压信号功率之比。
据此改变多普勒频率，得到200组对应的相对增益数据，将其与理论曲线绘在一起（完整程序见附录-程序3），得到的结果如下图：
Fig. 19 多普勒敏感造成增益损失的仿真
仿真结果与理论曲线相符，并且可以看到，在f_d= 173KHz时，增益损失到0，已经完全没有脉压增益。
6.2多普勒容限
理论：多普勒效应除了使得脉压峰值下降外，还会造成旁瓣提高，恶化信号性能。多普勒容限定义为一个区间，在此区间内，虽然由于多普勒敏感效应而不能产生最大输出信噪比，但相关积累后仍然可以用于信号检测。
研究多普勒容限的关键是得到不同多普勒频率下的主瓣旁瓣比。设脉压后的序列为$${s_{MF}}[k,{f_d}]$$，其中fd是多普勒频率参数。脉压主旁瓣比与目标距离无关，故可以假设在${s_{MF}}[k,{f_d}]$处有最大匹配峰值。所以主瓣高度是${s_{MF}}[{k_o},{f_d}]$，旁瓣存在于主瓣之后126个样点中，即旁瓣的范围为：$${s_{MF}}[k,{f_d}]{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({k_0} < k < {k_0} + 127)$$，主旁比MSR定义为主瓣能量与峰值旁瓣能量之比，即
$$\eqalign{
& MSR = 20\log \left( {\frac{{MainLobe({f_d})}}{{SideLob{e_{\max }}({f_d})}}} \right) \cr
& {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = 20\log \left( {\frac{{\left| {{s_{MF}}[{k_0},{f_d}]} \right|}}{{\max (\left| {{s_{MF}}[k,{f_d}]} \right|,{\kern 1pt} {k_0} < k < {k_0} + 127)}}} \right) \cr} $$[19]
由于MSR与目标距离无关，仿真时设目标距离为0,则主瓣出现在${s_{MF}}[0,{f_d}]$处。
根据式[19]，统计多普勒频率${s_{MF}}[0,{f_d}]$（也即Fig.16中脉压增益损失到0的第一个点的多普勒频率）时的MSR, 得到结果表示在下图中，统计所用程序见附录-程序4。
主瓣幅值按照Fig.16中的规律因为多普勒敏感而下降，峰值旁瓣起初保持在9作用，在频率大于40KHz后开始增大。MSR在53.6KHz时下降3dB。
7 距离分辨力和速度分辨力
7.1 距离分辩
距离分辨力由指同一方向上两个大小相等点目标之间的最小可分辨距离，分辨力d0取决于脉压后的脉冲宽度${\tau _0}$，而二项编码信号脉压后的时宽${\tau _0}$已经验证就等于子码宽度\({\tau _0}\)，所以$$d = \frac{{c{\tau _0}}}{2} = \frac{c}{{2B}}$$，B是脉压输出信号的带宽。代入数据计算得到距离分辨力$${d_0} = \frac{{c\tau }}{2} = \frac{{3e8 \times \frac{1}{{22e6}}}}{2} = 6.82m$$，所以理论最小距离分辨能力是d_0=6.82m. 仿真时设置
仿真：仿真时两个目标速度都设为10m/s; 一号目标距离设定为200m恒定，二号目标距离与一号目标依次相差0.5倍、1倍、1.5倍、2倍最小分辨距离，即依次设距离参数为：200+0.5d_0=203.41m，200+d_0=206.82m，200+1.5d_0=210.23m和200+2d_0=213.64m。仿真结果展示如下：
两目标相差0.5距离 两目标相差距离
两目标相差1.5距离 两目标相差2距离
Fig. 20 距离分辨率的仿真
仿真说明当两目标距离相差大于1.5倍最小分辨距离时可以分辨。
7.2 速度分辨
5.2节中已经论述过，所用m序列长度和占空比决定了多普勒频率幅度包络的采样率，根据奈奎斯特定理可以得出多普勒频率最大测量范围，结果是\(\left| {{V_d}} \right| \leqslant 129.86\)m/s. FFT是DFT的快速算法，而DFT就是对DTFT的采样，采样点数就是作FFT的点数。FFT点数决定了速度分辨率。本次仿真中，对171点FFT通道的数据末尾补零后作2^8=256点FFT，数据末尾补零并不能提高分辨力，但可以提高判断的精度。对171点数据作FFT，得到171个测速区间，平分2*129.86的测速范围，所以测速分辨率$${V_{d0}} = \frac{{129.86 \times 2}}{{171}} = 1.52m/s$$.
测速分辨率也可以通过公式:\({V_{d0}} = \frac{\lambda }{{2T}}\) 得到，其中\(\lambda \)是载波波长，T是相干时间。此处\(\lambda \)=0.03m，171个周期信号积累，即T=171*(1270*1/22e6) =9.87ms, 代入公式中，得到\({V_{d0}} = \frac{\lambda }{{2T}} = 0.03/(2*9.87e - 3) = 1.519\)m/s, 结果是相同的。
与仿真距离分辨时的做法类似，控制两目标距离都为200m恒定，目标一速度分别为80，仿真目标二的速度分别为80+0.5\({V_{d0}}\)=80.76m/s、80+\({V_{d0}}\)=81.52m/s、80+1.5\({V_{d0}}\)=82.28m/s、80+2\({V_{d0}}\)=83.04m/s的情况。结果如下：
两目标相差0.5速度 两目标相差1速度
两目标相差.1.5速度 两目标相差2速度
Fig. 21 速度分辨率的仿真
仿真说明，两目标速度相差达到1.5时可以较好地分辨。
8. 大目标掩盖小目标
8.1 旁瓣掩盖的仿真
对回波脉压会产生距离旁瓣；作FFT后，由于频率泄露也会产生速度旁瓣。当同时接收到两个目标的回波信号，且它们的信号幅值相差过大时，大目标的旁瓣会掩盖小目标的主瓣，从而影响小目标的识别。
仿真时两个目标的幅度A，距离d(m)，速度v(m/s)信息表示成参数向量(A,d,v), 以下给出了大目标被速度旁瓣和距离旁瓣掩盖的情况。
速度旁瓣掩盖：
距离旁瓣掩盖：
Fig. 22 大目标掩盖小目标情况仿真
8.2 抑制速度旁瓣的措施-加窗
在8.1节大目标掩盖小目标的仿真中，发现速度旁瓣幅值高、衰减慢，严重影响小目标的识别能力。
对无限长正弦波做傅立叶变换理应得到纯净的频谱冲击函数。但在实际数字信号处理中，只能取有限长的样点数，相当于用矩形窗截取了无限长序列。矩形窗具有很宽的频谱分量，频域特性是sinc函数包络。因此在本次信号处理中对171点信号作FFT会造成频谱泄露，产生旁瓣。对FFT之前的信号加窗可以抑制旁瓣。
海明(Hamming)窗是一种升余弦窗函数，表示式为：
$$w[n] = 0.54 + 0.46\cos \left( {\frac{{\pi n}}{M}} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - M \leqslant n \leqslant M$$
海明窗可以有效降低旁瓣，而主瓣扩展和增益损失相对较小，适用于当前场景。仿真条件设为10dB信噪比，两目标的幅度、距离、速度参数分别设为(25,200,80)和(1,200,95)，这也是之前仿真的小目标被大目标距离旁瓣掩盖的情况。
比较加窗前后对距离重叠目标的处理结果，如下所示：
Fig. 23 加窗前后频谱对比
可见，加窗后作FFT可以有效抑制旁瓣，这有利于小目标的检测；但同时主瓣宽度增大，信号幅值减小，这又将降低速度分辨率。
采用与式[15]相同的方法，统计输出信噪比
$$SNR = \frac{{{C^2}{{[m,n]}_{\max }}}}{{\sum\nolimits_{k = 50}^{999} {{C^2}[120,k]/950} }}$$
统计加窗前后大目标的信噪比，考察加窗后对增益的影响。计算得：
加窗前的输出信噪比SNR1= (5.29e5)^2 / 1.37e5=2.04e6=63dB;
加窗后的输出信噪比SNR2=(2.88e5)^2 / 5.17e4=1.60e6=62dB.
从仿真结果可见，加窗之后虽然信号幅值降低了，但噪声也被压制，输出信噪比并没有过多损失。采用海明窗来消除速度旁瓣是切实有效的，
9 Appendix
附全部程序：
程序1：主要程序
%-------2018年6月1日---------------
% Ver1.1
%修改说明：修改了几处硬编码，初始时设置的参数可以应用作用到整个程序，
%现在处理部分与配置参数不再耦合
clc;
clear all;
close all;
tic;
%%--设置两个目标的参数------------
A1=1;%回波强度
d1=1570;%距离m
v1=60;%速度m/s
A2=1;
d2=200;
v2=22;
snr= 10;%回波信噪比(可变)
%%--设置两个目标的参数--------------------
%%--根据所设参数计算相应的频偏
fd1=200/3*v1;%对应10GHz下的多普勒频偏,用于模拟生成回波
fd2=200/3*v2;
%-------雷达信号可调参数设置----------
fm=22e6;%码频
fc=10e9;%载频
t=10e-3;%相干积累总时宽
%-------雷达信号可调参数设置------------------
%------由雷达信号参数计算待用的参数
fs=fm*1270;%多普勒余弦作ft时的采样率
T=1/fm*127;%周期
N=round(10e-3/(127*1/22E6)); %积累时间不大于10ms，计算得到127长度M序列的延拓次数(N=173)
distancePerCodeword= 1/fm * 3e8 /2 ;%每个采样点代表的距离
%--生成基带信号
%----生成127点m序列
n = 7; % 阶次
p = 2^n -1; % 循环周期
ms = idinput(p, 'prbs')';% Gives a p-length pseudorandom sequence
msfull=repmat(ms,1,N);%将127长度的M序列周期延拓N次
figure
plot(autocorr(msfull,1024),'r');%●作M序列循环自相关函数,最大长度1024点
title('127-M序列的自相关函数'),xlabel('样点n'),ylabel('Rx')
%至此已经产生了127点M序列ms，并周期延拓N=1732次，存储在msfull
%以下产生占空比方波，与ms1点乘
y=(1+square((0:2*pi/1270:2*pi*N/10),10))/2;%产生占空比方波存入y
y(N*127+1:end)=[ ];%删除多余的点，使y序列的长度与msfull的长度相等
%产生了1270点周期，10%占空比的方波
% figure
% stairs(y)
% ylim([-2 2])
%——————**单目标情况————————————————
%-*1.生成回波信号
i=linspace(0,N*127-1,N*127);
c1=A1*exp(j*(2*pi*fd1/fm*i+2*pi*fc*2*d1/3e8));%模拟由于距离和速度产生的多普勒频移和相移
send=msfull.*y;%发射信号
temp1=c1.*msfull.*y;%加上频移后的信号存入temp1
%clear y;
%echo signal
echo1=zeros(1,N*127);
samplePointsPerUnitDistance=2 / 3e8 *fm; %距离d1带来延时，再换算到码字(样点)数
echo1(1,1+round(samplePointsPerUnitDistance*d1):N*127)=temp1(1,1:(N*127-round(samplePointsPerUnitDistance*d1)));%加入时延
echoSingle=awgn(echo1,snr,'measured');%adds white Gaussian noise, 得到单目标回波
clear temp1;
figure;%发射和接收信号作图
subplot(211);
stairs(send);title('发射的基带信号'),xlabel('样点数');
xlim([0 10000]),ylim([-2 2]);
subplot(212);
stairs(echoSingle);title('下混频后的回波信号(时移,频移,噪声)'),,xlabel('样点数');
xlim([0 10000]);
% plot(fftshift(abs(fft(echo))));
% xlim([0 10000])
%-*2.脉压，使用127点m序列作为匹配波形，得到匹配滤波后的数据存在MF1---
local=fliplr(repmat(ms,1,1));%反褶
MF1=conv(local,echoSingle);
MF1=MF1(127*1:1:end);%去除前1*127点数据
figure;%●作脉压后的图像
subplot(211);
x=linspace( 0,length(MF1)*distancePerCodeword,length(MF1) );%生成对应的距离坐标以绘图
plot(x,MF1);title('脉压后的图像'),xlabel('距离'),ylabel('幅值');
xlim([0 10000])
subplot(212);
%MF1=MF1./(max(max(MF1)));%归一化
plot(x,20.*log10(MF1) );title('脉压后的图像(对数坐标)'),xlabel('距离'),ylabel('db');
% semilogy(x,MF1);
% ylim([-10 500])
xlim([0 10000])
% figure;%对脉压后的图像作fft
% plot(abs(fft(MF1)));
% xlim([0 1000])
%-*3.距离门重排---数据存入RangeCellA[]，1270*173的矩阵，1270个距离门，173点幅值数据
%RangeCellA=zeros(1270,173);
for i=1:1270
RangeCellA(: ,i)=MF1(i:1270:1270*170+i);
end
figure;mesh(1:distancePerCodeword:1270*distancePerCodeword,1:171,real(RangeCellA)),xlabel('距离'),ylabel('样点序号');
%pause
%-*4.FFT---对同一距离门的数据作fft，得到速度量-------
NFFT=2^nextpow2(length(RangeCellA(:,1)));
%生成Hamming窗函数
W1=0.54+0.46*cos(pi.*linspace(-1,1,171))';%----海明窗
W1=repmat(W1,1,1270);
FFTA=fft(W1.*RangeCellA,NFFT);%加窗后对每列作FFT,点数为NFFT=256
figure;
subplot(211)
surfl(1:1:1270,1:NFFT,abs(FFTA)),title('作FFT'),xlabel('距离门序号'),ylabel('FFT样点');
shading flat;
%FFTAspeed=FFTA(1:NFFT/2+1,:);
subplot(212)
fs2=1/(10*127*1/fm);%距离门重排后实际上采样率降低
maxVelocity=fs2/2 * 3e8 /(2*fc);
speed=[2*maxVelocity/NFFT*(0:1:NFFT/2-1) 2*maxVelocity/NFFT*(-NFFT/2+1:1:0) ];%FFT后的对应速度维
surfl(1:distancePerCodeword:1270*distancePerCodeword,speed,abs(FFTA) );
shading flat;
title('换算到速度维'),xlabel('距离 m'),ylabel('速度 m/s');
sound(sin(1:1:300));
disp('单目标时的处理完成,press any key to continue..');
toc;
pause;
%---------------单目标处理到此结束。---
%-------------**开始双目标信号处理-----
tic
%--*1.生成第二个目标的回波,-----
i=linspace(0,N*127-1,N*127);
c2=A2*exp(j*2*pi*fd2/fm*i+j*2*pi*fc*2*d2/3e8);%多普勒频移和相移
temp2=c2.*msfull.*y;
clear y;
echo2=zeros(1,N*127);
echo2(1,1+round(samplePointsPerUnitDistance*d2):N*127)=temp2(1,1:(N*127-round(samplePointsPerUnitDistance*d2)));%加入时延
echoDouble=awgn(echo2+echo1,snr,'measured');%adds white Gaussian noise,得到了双目标的回波，存入echoDouble
figure;%双目标时回波信号作图
subplot(211);
stairs(send);title('发射的基带信号'),xlabel('样点数'),ylabel('幅值');
xlim([0 10000]),ylim([-2 2]);
subplot(212);
stairs(echoDouble),title('双目标的回波'),xlabel('样点数'),ylabel('幅值');
xlim([0 10000]);
%-*2.对回波脉压,匹配滤波后的数据存入MF2---------
MF2=conv(local,echoDouble);
MF2=MF2(127*1:1:end);%去除前1*127点数据
clear echoDouble;
clear echoSingle;
figure;%●作脉压后的图像
x=linspace( 0,length(MF2)*75/11,length(MF2) );%生成对应的距离坐标
plot(x,MF2);title('脉压后的图像'),xlabel('距离');
xlim([0 10000]);
%-*3.距离门重排---数据存入RangeCellA[]，1270*173的矩阵，1270个距离门，173点幅值数据
RangeCellB=zeros(171,1270);%双目标距离门矩阵
for i=1:1270
RangeCellB(: ,i)=MF2(i:1270:1270*170+i);
end
figure;mesh(1:distancePerCodeword:1270*distancePerCodeword,1:171,real(RangeCellB) ),title('按距离门排列'),xlabel('距离'),ylabel('样点序号');
%pause
%-*4.FFT---对同一距离门的数据作fft，得到速度量-------
NFFT=2^nextpow2(length(RangeCellB(:,1)));
FFTB=fft(W1.*RangeCellB,NFFT);%加Hamming窗后对每列作FFT,点数为NFFT=256
figure;
subplot(211)
surfl(1:1:1270,1:NFFT,abs(FFTB)),title('作FFT后'),xlabel('距离门序号'),ylabel('FFT样点');
shading flat;
%FFTBspeed=FFTB(1:NFFT/2+1,:);%取前一半FFT的值|改为复数形式后可以测到负速度
subplot(212)%换算到速度维
surfl(1:distancePerCodeword:1270*distancePerCodeword,speed,abs(FFTB) );
shading flat;
title('换算到速度维'),xlabel('距离 m'),ylabel('速度 m/s');
sound(sin(1:1:200));
toc
10 Reference
[1] 曹丽娜，通信原理（第七版）[M] , 国防工业出版社，2015: 390-395
[2] 顾红, 刘国岁. 随机二相编码连续波雷达的研究[J]. 电子学报, 1995(12):71-74.
[3] 季鹏，解决多普勒敏感问题的PRC-CW雷达信号处理的研究[D] 2009: 13-14
推荐阅读