0 概述

近年来，随着无线通信技术的不断发展^[1-3]，用于用户通信的授权频谱资源日趋紧缺，若无线通信系统能通过非授权频谱进行用户通信，则其可用的频谱资源将大幅增加。为支持长期演进增强（Long Term Evolution-Advanced，LTE-A）系统在非授权频谱上通信，长期演进（Long Term Evolution，LTE）系统引入许可辅助接入（Licensed⁃Assisted Access，LAA）^[4]、非授权长期演进（Long Term Evolution-Unlicensed，LTE-U）^[5]以及MulteFire^[6]等关键技术。当前5G非授权无线频谱独立新空口（5G New Radio in Unlicensed Spectrum Stand Alone，5G NR-USA）技术主要用于实现通信系统在非授权频谱上通信。5G NR-USA与LAA、LTE-U的最大区别在于其通过非授权频谱独立进行通信，而不使用授权频谱作为锚点传输控制信令，因此，需要用户终端独立地采用非授权频谱进行频谱感知，这要求终端使用的频谱感知方法复杂度较低。

传统频谱感知方法^[7]主要包括能量检测法、匹配滤波器检测法以及循环平稳检测法^[8]等。能量检测法按照接收信号能量大小进行检测，具有复杂度低与实时性高的优点，但其在信噪比低时鲁棒性较差，无法进行信号类型识别。匹配滤波器检测法要求接收端使用与发送端一致的滤波器，且每种信号均对应一个滤波器，这在实际应用中难以实现。循环平稳检测法根据信号的二阶循环统计量进行判断，在低信噪比时仍具有良好的检测性能，且通过信号的循环谱可得到具体特征^[9]，但是该方法计算复杂度较高，不适合在终端应用。

针对上述问题，本文提出一种基于循环平稳特性的循环谱切面检测方法，通过分析切面中频率的约束关系，建立循环谱单一切面检测信号的循环平稳特征。由于在非授权频段中5 250 MHz~5 350 MHz和5 470 MHz~5 725 MHz为雷达可用频段，在该频段进行通信时需进行频谱感知，因此本文针对最常见的线性调频（Linear Frequency Modulation，LFM）雷达信号^[10]和正交频分复用（Orthogonal Frequency Division Multiplexing，OFDM）信号^[11-12]进行循环平稳特征检测。

1 循环平稳特性

循环平稳特性^[13]指信号的统计特性按照一定周期平稳变化，该特性通常由通信系统中传输信号的采样、调制、添加循环前缀等人为操作引入。由于噪声通常不具备这一循环平稳特性，因此在低信噪比下根据循环平稳特征检测可很好地区分信号与噪声。

如果1个非平稳信号的N阶统计量随时间呈周期性变化，则称该信号具有N阶循环平稳特性。对于具有二阶循环平稳特性的信号x(t)，当采样点的数量N趋于无穷大时，可用时间平均值表示统计平均值，该信号的自相关函数表示为：

$\begin{array}{*{20}{l}} {{R_x}\left( {t, \tau } \right) = \mathop {{\rm{lim}}}\limits_{T \to \infty } {{\hat R}_x}{{\left( {t, \tau } \right)}_T} = \mathop {{\rm{lim}}}\limits_{N \to \infty } \frac{1}{{2N + 1}} \cdot }\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathop \sum \limits_{n = - N}^N x\left( {t + nT + \frac{\tau }{2}} \right){x^{\rm{*}}}\left( {t + nT - \frac{\tau }{2}} \right)} \end{array}$

(1)

其中，T为采样周期。由于x(t)具有二阶循环平稳特性，即x(t)的二阶统计特性呈周期性变化，因此R_x（t, τ）是周期为T的函数，可用傅里叶级数表示如下：

$ {R_x}\left( {t, \tau } \right) = \mathop \sum \limits_{m = - \infty }^{ + \infty } R_x^{\frac{m}{T}}\left( \tau \right){{\rm{e}}^{{\rm{j}}2{\rm{ \mathsf{ π} }}t\frac{m}{T}}}, m = 0, \pm 1, \cdots , \pm n $

(2)

$ R_x^{\frac{m}{T}}\left( \tau \right) = \mathop {{\rm{lim}}}\limits_{T \to \infty } \frac{1}{T}\int_{ - \frac{T}{2}}^{\frac{T}{2}} {x\left( {t + \frac{\tau }{2}} \right){x^{\rm{*}}}\left( {t - \frac{\tau }{2}} \right){{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\alpha t}}{\rm{d}}t} $

(3)

其中，$R_x^{\frac{m}{T}}$为傅里叶系数，α=m/T为该信号的循环频率，$R_x^\alpha \left( \tau \right)$为信号x(t)的循环自相关。令u(t)=x(t)·e^-jπαt，v(t)=x(t)·e^jπαt，使用换元法将$R_x^\alpha \left( \tau \right)$表示为互相关形式：

$ R_x^\alpha \left( \tau \right) = {R_{uv}}\left( \tau \right) = \frac{1}{T}\int_{ - \frac{T}{2}}^{\frac{T}{2}} {u\left( {t + \frac{\tau }{2}} \right){v^{\rm{*}}}\left( {t - \frac{\tau }{2}} \right){\rm{d}}t} $

(4)

根据互相关函数与功率谱密度互为傅里叶变换对，得到x(t)的循环谱密度函数$表达式为：

$ \begin{array}{*{20}{l}} {S_x^\alpha \left( f \right) = {S_{UV}}\left( f \right) = \int_{ - \infty }^{ + \infty } {{R_{uv}}\left( \tau \right){{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}f\tau }}{\rm{d}}\tau } = }\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\int_{ - \infty }^{ + \infty } {R_x^\alpha \left( \tau \right){{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}f\tau }}{\rm{d}}\tau } } \end{array} $

(5)

由式（2）可知，信号x(t)的循环自相关与循环谱密度函数仅在离散点的周期整数倍位置有数值，而其在离散点其他位置的数值恒等于0，即循环谱密度函数在循环频率轴上以α=m/T的间隔出现峰值。图 1为码元速率为2 000波特的LFM雷达信号循环谱在频率f=f_c处的切面以及加性高斯白噪声（Additive White Gaussian Noise，AWGN）循环谱在f=0处的切面。可以看出，雷达信号在α=0和α=m/T处存在峰值，峰值间隔为信号的码元速率1/T，因为数字信号处理时采用较多窗函数，造成频谱存在泄露和混叠的情况，所以其在α≠m/T处数值不为0。加性高斯白噪声不具有二阶循环平稳特性，当α≠0时其循环谱密度函数${S_x^\alpha \left( f \right) = 0}$，同样由于数字信号处理中使用较多窗函数，因此其在α≠m/T处不为0。

2 FAM的优化

随着对循环平稳特性研究的不断深入，研究人员提出频域平滑方法（FSM）^[14]、时域平滑方法、快速傅里叶变换累加方法（Fast Fourier Transformation Accumulation Method，FAM）^[15]、分段谱相关（SSCA）方法^[16]等信号循环谱检测方法。其中：频域平滑方法检测性能最好，但其观测时间长且计算复杂度较高，应用实时性较差；时域平滑方法通过加窗截断待测信号减少数据长度，但是由于对采样数据进行加窗与在时域上进行平均化处理这两个过程不能并行，不适用于实际工程；FAM和SSCA方法均基于时域平滑方法进行优化，其中，SSCA方法使用频域和时域混合计算，造成其检测性能较时域平滑方法要差，FAM因为使用两次FFT简化计算，所以其检测性能受到的影响较小。此外，目前关于FFT模块的研究已较成熟，因此，FAM更适用于DSP或FPGA工程开发。本文以FAM作为基础方法，进一步优化其在工程应用中的计算复杂度。

2.1 传统FAM

FAM^[15]是基于时域平滑方法利用快速傅里叶变换（Fast Fourier Transformation，FFT）简化计算的方法。时域平滑方法中循环谱计算公式^[17]如下：

$ S_{{x_T}}^\alpha {\left( {n, f} \right)_{{\rm{\Delta }}t}} = \frac{1}{T}{\left\langle {{X_T}\left( {n, f + \alpha /2} \right)X_T^{\rm{*}}\left( {n, f - \alpha /2} \right)} \right\rangle _{{\rm{\Delta }}t}} $

(6)

其中，〈·〉_∆t表示在时间∆t内取平均值，X_T(n, f)表示信号x(n)的复包络。由式（6）可得到信号x(n)的频率分量X_T(n, f+α/2)和X_T(n, f-α/2)在时间∆t内的平均相关程度，使用FFT计算得到：

$ {X_T}\left( {n, f} \right) = \mathop \sum \limits_{r = - \frac{{N'}}{2}}^{\frac{{N'}}{2}} a\left( r \right)x\left( {n - r} \right){{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}f\left( {n - r} \right){T_s}}} $

(7)

其中，T_s为采样间隔，a(r)为长度T=N'·T_s的数据衰减窗。将S^α_xT(n, f)_△t沿α轴频移ε，使用FFT代替时域平均化过程^[15]，计算公式如下：

$S_{{x_T}}^{{\alpha _i} + \varepsilon }{\left( {n, {f_j}} \right)_{{\rm{\Delta }}t}} = \mathop \sum \limits_r {X_T}\left( {n, {f_k}} \right)X_T^{\rm{*}}\left( {n, {f_l}} \right){g_c}\left( {n - r} \right){{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\varepsilon r{T_s}}}$

(8)

其中，ε为频移因子。为简化计算，在数据序列上进行时域抽取以降低需计算的采样点数，计算公式如下：

$ S_{{x_T}}^{{\alpha _i} + \varepsilon }{\left( {rL, {f_j}} \right)_{{\rm{\Delta }}t}} = \mathop \sum \limits_r {X_T}\left( {rL, {f_k}} \right)X_T^{\rm{*}}\left( {rL, {f_l}} \right){g_c}\left( {n - r} \right){{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\varepsilon r{T_s}}} $

(9)

将频移因子ε离散化，令ε=q∆α，代入式（9）得到FAM计算公式为：

$ \begin{array}{l} S_{{x_T}}^{{\alpha _i} + q{\rm{\Delta }}\alpha }{\left( {rL, {f_j}} \right)_{{\rm{\Delta }}t}} = \mathop \sum \limits_r {X_T}\left( {rL, {f_k}} \right)·\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;X_T^{\rm{*}}\left( {rL, {f_l}} \right){g_c}\left( {n - r} \right){{\rm{e}}^{\frac{{ - {\rm{j}}2{\rm{ \mathsf{ π} }}rq}}{P}}} \end{array} $

(10)

其中，L为时域抽取的采样点数，P=N/L为采样点数的实际计算数目，∆α=1/∆t为FAM的循环频率分辨率，α_i和f_j的计算公式如下：

$ {\alpha _i} = {f_k} - {f_l} = \left( {k - l} \right)\left( {\frac{{{f_s}}}{{N'}}} \right), k, l \in \left[ { - \frac{{N'}}{2}, \frac{{N'}}{2} - 1} \right] $

(11)

$ {f_j} = \frac{{{f_k} + {f_l}}}{2} = \frac{{k + l}}{2}\left( {\frac{{{f_s}}}{{N'}}} \right), k, l \in \left[ { - \frac{{N'}}{2}, \frac{{N'}}{2} - 1} \right] $

(12)

为最小化频移误差带来的影响，通常只估计在[-∆α/2, ∆α/2]范围内循环谱频移的数值，其中∆α=f_s/N'，循环谱频移范围值q∆α∈ [-∆α/2, ∆α/2]。

FAM的具体流程如图 2所示，使用P点的FFT代替累加与低通滤波操作。FAM使用FFT代替时域平滑方法中每一段的平滑过程，可缩短窗口平滑时间，并充分利用FFT并行化优势，进一步提升计算速度。因此，FAM是循环平稳检测的主流方法。

2.2 f-切面法

使用FAM构建的雷达信号循环谱包括该信号全部信息，由式（10）可知，计算X_T(n, f_j)的相关程度时需遍历所有频率f_j，对每个频率f_j需遍历计算全部循环频率α_i，导致计算复杂度较高。由图 1可知，符合循环平稳特性的信号在α=m/T处出现峰值，因此根据该特性得出：无需计算待测信号全部循环谱，仅计算某个f_j切面的循环谱即可进行信号检测。为进一步提高效率，只计算f_j=0切面的循环谱进行信号检测。

由式（12）可知，在计算f_j=0时，对于每个k(k∈[-N'/2, 0])，都有唯一的l与其对应并满足l=-k，因此，仅需N'组的采样数据就可得到所有f_j=0对应的循环谱密度函数。图 3为所有k和l组合的矩阵，其中方框内的组合均满足k+l=0。

图 4为传统FAM和f-切面优化FAM（以下称为f-切面法）的X_T(n, f_j)选取方案（向上和向下的箭头表示选取顺序）。在得到每段信号的复包络X_T(n, f_j)后，图 4中标号为1到N'的采样数据所对应k、l的取值范围分别为[-N'/2, N'/2-1]，传统FAM对于每个k均遍历N'个l进行共轭相乘，而f-切面法结合限制条件k+l=0后，对于每个k仅选择1个l进行共轭相乘。由k，l∈[-N'/2, N'/2-1]得到k的实际取值范围为[-N'/2+1, N'/2-1]。

对于$\forall {f_k} \in \left[ { - N'/2 + 1, 0} \right]$，存在如下特性：

${X_T}\left( {rL, {f_k}} \right)X_T^{\rm{*}}\left( {rL, {f_{ - k}}} \right) = {\left( {{X_T}\left( {rL, {f_{ - k}}} \right)X_T^{\rm{*}}\left( {rL, {f_k}} \right)} \right)^{\rm{*}}}$

(13)

通过使用式（13）中共轭运算可减少式（10）的计算量。在实际设计时，为增加FFT效率，N'设置为2的整数次幂。

循环频率坐标α₀的计算公式为：

${\alpha _0} = {\alpha _i} + q{\rm{\Delta }}\alpha $

(14)

将式（11）代入式（14）得到α₀，再将式（10）计算得到的$S_{{x_T}}^{{\alpha _i} + q{\rm{\Delta }}\alpha }{\left( {nL, {f_0}} \right)_{{\rm{\Delta }}t}}$以α₀为x轴绘制得到信号的二维循环谱。

2.3 α-切面法

使用f-切面法构建f=0切面虽然可以有效减少计算量，然而由于不同信号具有不同循环平稳特征，因此OFDM等部分信号在f=0切面不会出现明显的谱峰。无循环前缀的OFDM信号^[18]循环谱密度函数为：

$ S_x^\alpha \left( f \right) = \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{4T}}{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\alpha {t_0}}}\left\{ {\mathop \sum \limits_{k = - N/2}^{N/2} {{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}k{\rm{\Delta }}f{t_0}}} \cdot } \right.}\\ {[Q\left( {f + {f_c} + k{\rm{\Delta }}f + \alpha /2} \right){Q^{\rm{*}}}\left( {f + {f_c} - k{\rm{\Delta }}f - \alpha /2} \right) + }\\ {\left. {Q\left( {f - {f_c} + k{\rm{\Delta }}f + \alpha /2} \right){Q^{\rm{*}}}\left( {f - {f_c} - k{\rm{\Delta }}f - \alpha /2} \right)]} \right\}, }\\ {\alpha = m/{T_s}}\\ {0, \alpha \ne m/{T_s}} \end{array}} \right. $

(15)

在f=0切面上，由于多个子载波造成频谱混叠而失去原始形状，因此在α=2f_c±m/T处不存在清晰谱峰。

单频信号中QPSK由于其正交分量与同向分量平衡造成在f=0切面无明显的谱峰，其循环谱密度函数^[19]为：

$ \begin{array}{l} S_x^\alpha \left( f \right) = \frac{1}{{2T}}\left[ {Q\left( {f + \frac{a}{2} + {f_c}} \right){Q^{\rm{*}}}\left( {f - \frac{a}{2} + {f_c}} \right)S_c^\alpha \left( {f + {f_c}} \right) + } \right.\\ \left. {\;\;\;\;\;\;\;\;\;\;\;\;\;Q\left( {f + \frac{a}{2} - {f_c}} \right){Q^{\rm{*}}}\left( {f - \frac{a}{2} - {f_c}} \right)S_c^\alpha \left( {f - {f_c}} \right)} \right]{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\alpha {t_0}}} \end{array} $

(16)

因此，使用f-切面法对上述两种信号的检测性较差。但依据循环平稳特征的定义，在f=±f_c切面上，循环平稳特征信号在α≠0处存在峰值。

根据循环谱的定义，信号的循环谱密度函数在α=0时将退化为信号的功率谱密度函数，并在对应的f=±f_c处出现峰值，这适用于所有具有循环平稳特性的信号，由此可估计出待测信号的载波频率f_c，进而采用f-切面法构建得到待测信号循环谱的f=±f_c切面，同时通过该方法检测类似OFDM或QPSK等在f=0切面无明显循环平稳特性的信号。

通过f-切面法构建待测信号循环谱的切面，即令α=0，则由式（11）和式（14）得到：

${\alpha _0} = \left( {{f_k} - {f_l}} \right) + q{\rm{\Delta }}\alpha = 0$

(17)

由2.1节中循环谱频移范围、Δa值和Δα值可得到参数q∈[-P/8, P/8)，式（17）转化为：

$ {\alpha _0} = \left[ {k - l + q \cdot \left( {\frac{4}{P}} \right)} \right] \cdot fs/N' = 0 $

(18)

由于k、s、q均为整数，因此只有当q=0时式（18）才成立，同时k=l。

由f-切面法可知，${X_T}\left( {n, f + \alpha /2} \right)X_T^{\rm{*}}\left( {n, f - \alpha /2} \right)$表示每段信号在频谱上间隔为α的相关程度。当α=0时，${X_T}\left( {n, f + \alpha /2} \right)X_T^{\rm{*}}\left( {n, f - \alpha /2} \right)$转换为|X_T (n, f)|²，即仅需计算每个k对应的X_T(n, f_k)模值的平方，再将其与q=0同时代入式（10）得到该信号循环谱的α=0切面。由图 3可知，满足k=l共有N'组，因此，通过计算N'次模值的平方可得到α=0切面的全部f_j组合。

图 5为OFDM信号在信噪比为0 dB时循环谱的α=0切面。可以看出，在f=±f_c处有2个峰值，由此估计出该OFDM信号的载波频率。由于OFDM信号是多载波信号，在f=±f_c处会偏移k∆f，因此只能估计得到不精确的载波频率f_c。

令式（12）中f_j等于估计的载波频率，即：

${f_j} = \frac{{{f_k} + {f_l}}}{2} = \frac{{k + l}}{2}\left( {\frac{{{f_s}}}{{N'}}} \right) = {f_c} = c\left( {\frac{{{f_s}}}{{N'}}} \right)$

(19)

由式（19）得到c=(k+l)/2，其中c为载波频率f_c在全部频率f_j中的位置。根据循环谱对称性可知，如果c > N'/2，令c=1+(N'-c)，则有2c < N'。

对f-切面法进行扩展，挑选出2c-1组满足约束条件k=2c-l的数据，可得到待测信号在f=f_c切面的循环谱。如图 6所示，对于每个k∈[1, 2c-l]只有1个l与之对应，因此，仅需进行2c-l组共轭相乘就可得到f=f_c切面的循环谱。

图 7为使用α-切面法构建的OFDM信号在信噪比为0 dB时循环谱的f=f_c切面。可以看出多个子载波相互混叠导致谱线不明显，但在α=0两侧各有1个明显的谱峰，谱峰位置与α=0的间隔大小数值上等于OFDM信号的码元速率1/T，这与文献[18]使用传统FAM得到的仿真结果一致。

在实际工程应用中，f-切面法用来检测大部分单频信号及常见的LFM雷达信号，若未知待测信号类型或需要检测OFDM等多载波信号，则可使用α-切面法构建信号循环谱的f=±f_c切面，并估计出信号的载波频率和码元速率。

2.4 信号检测

由循环平稳特性分析可知，在使用f-切面法和α-切面法所得信号循环谱的f=0或f=f_c切面上，信号会在α≠0处出现峰值，而噪声只会在α=0处出现峰值，因此，可通过统计超过门限的峰值数目来检测是否存在信号，若满足检测统计量$S_x^{{\alpha _0}}\left( {rL, {f_0}} \right) > \gamma $（γ为检测门限）的α₀超过1个，则认为该频段内存在信号。对于加性高斯白噪声，$S_x^{{\alpha _0}}\left( {rL, {f_0}} \right)$服从自由度为2的卡方分布，可根据其累积分布概率选择满足给定虚警率时的门限γ，信号检测流程如图 8所示。

综上所述，使用f-切面法和α-切面法可直接得到待测信号循环谱的切面作为$S_x^{{\alpha _0}}\left( {rL, {f_0}} \right)$进行信号的循环特征检测，而无需计算全部循环谱，从而避免建立不必要的切面，可加快计算速度，保证频谱检测的实时性。

3 仿真与结果分析

使用MATLAB软件搭建仿真链路，分别使用f-切面法和α-切面法构建f_j=0和f_j=f_c的循环谱切面进行信号检测，仿真参数设置如表 1所示。

表 2为f-切面法、α-切面法与传统FAM在式（10）中计算复杂度的对比。可以看出：与传统FAM相比，f-切面法在共轭相乘和第二次FFT中乘法运算次数均减少到1/N'，使用上述参数进行性能仿真，共减少3×10⁹次乘法运算；由于α-切面法需多构建1次α=0切面，因此其计算量是f-切面法的两倍，但是当载波频率与N'相差较大（2c﹤﹤N'）时，α-切面法的计算量与f-切面法接近。与传统FAM相比，f-切面法和α-切面法的计算复杂度均大幅降低。

为分析不同信号检测方法之间的性能差异，分别将f-切面法、α-切面法、传统FAM^[15]和使用恒虚警率（Constant False Alarm Rate，CFAR）自适应门限的能量检测法^[20]（CFAR-ED）在不同信噪比和噪声估计误差（Z）下得到的检测概率进行对比，恒虚警率为0.05，结果如图 9所示。可以看出：随着信噪比逐渐升高，f-切面法和α-切面法的检测概率与传统FAM逐渐接近；当信噪比为－18 dB时，f-切面法的检测概率达到传统FAM的95%以上，且高于α-切面法的检测概率；当信噪比低于－16 dB时，f-切面法和α-切面法的检测概率低于传统FAM，这是因为当信噪比过低时，单一切面上信号特征会被噪声淹没；当信噪比达到－16 dB后，f-切面法和α-切面法的检测概率均与传统FAM一致；α-切面法较f-切面法检测性能略差，这是因为噪声在循环谱上集中在α=0切面，会降低所估计信号载波频率的精确度。由于循环平稳检测与噪声估计无关，因此其不受噪声估计误差的影响，当噪声估计误差大于0.017 dB时，能量检测法的检测性能劣于f-切面法和α-切面法。

综上所述，与传统FAM相比，本文提出的f-切面法和α-切面法在不影响检测性能的情况下，可大幅降低循环平稳检测法的计算复杂度。

4 结束语

在非授权频谱检测中，循环平稳检测法因计算复杂度高而难以在实际工程中推广应用。针对该问题，本文提出一种利用改进FAM的循环谱切面检测方法，采用构建循环谱单一切面的方式设计f-切面和α-切面两种FAM优化方案，对不同类型信号进行循环平稳特征检测。实验结果表明，与传统FAM相比，该方法的检测性能在低信噪比时略有下降，但计算复杂度大幅降低，检测性能较能量检测方法有大幅提升。下一步将具体分析不同雷达信号的循环平稳特征，充分利用待测信号的循环谱信息，在工程应用中实现信道内多种信号的盲识别与参数估计。