Frequency Locator Polynomial Based Fast Algorithm for Sparse Aliased Spectrum Recovery

1.   引言

随着无线通信设备和用户的剧增，无线频谱感知技术研究需求增长，特别是宽带频谱感知(Wideband Spectrum Sensing, WSS)技术的研究在近几年成为热点。在实践中，WSS给采样设备带来巨大挑战，若采用传统的奈奎斯特采样率捕获宽频信号，比如GHz频谱，需要至少2 GS/s采样率，该类采集卡不仅价格昂贵，且分辨率低。在此种背景下，研究者们寻求一种可以采用欠采样实现的宽频感知技术，此类技术主要包括有压缩感知(Compressed Sensing)技术^[1–5]、互素采样(Co-prime Sampling)技术^[6–8]以及延时错位采样(shifted sampling)技术^[9–11]。

压缩感知技术假设信号频谱是稀疏的，以低于奈奎斯特采样率采集信号，目前该类研究在理论上已经相对较完善。值得注意的是，以色列的Eldar教授团队提出了一种实时压缩采样技术实现方案^[1,2]。该方案的核心技术是一种所谓的宽带调制转换器(Modulated Wideband Converter, MWC)，MWC的主要功能是通过一定规则将频谱混叠到基带，然后使用对应的解混叠算法恢复出原始信号。然而，与其他压缩感知算法相同，该方法的计算以及硬件复杂度依然较高。

互素采样技术主要利用中国剩余定理，即在使用采样率互素的采集通道的前提下，原始频点在一路欠采样频谱中混叠，则在另一路中是不混叠的。文献[6]结合互素采样和子采样(MultiCoset Sampling, MCS)，实现一种多节点联合感知，可见实现成本较高。互素采样算法在计算复杂度上相比压缩感知有较大优势。但对时钟精度要求高，比如感知GHz的频谱，则时钟精度需要达到至少 $1 \times {10^{ - 9}}$ s，而且需具备产生多个互素的时钟频率，此外，还需要不同类型的采集卡，因此硬件成本高，不易实现。

延时错位采样技术相对上述两种实现，具有硬件复杂度低的优势。文献[9]提供了一种基于稀疏傅里叶变换的不同延时的多通道宽带频谱感知方案-BigBand，采用3块50 MS/s的ADC捕获和解译0.9 GHz的频谱，实验结果表明该算法具有较好的宽带稀疏频谱感知能力。但BigBand对噪声很敏感，算法的鲁棒性较低，同时算法复杂度相对较高。

本文针对BigBand进行改进，通过重新设计混叠频谱解算算法，使得计算复杂度有效降低，且混叠频谱恢复的准确率有明显提升，仿真结果也验证了该算法的可行性。

2.   问题描述

首先考虑无噪声污染的复信号情形。已知时域信号的欠采样会导致频域的混叠，设感知信号 $x$ 的频谱 $X$ 带宽为 $B$ ，欠采样倍数为 $\alpha$ ，则欠采样后的时域信号样本 ${\hat x_i} = {x_{i + mB/\alpha }}$ , $i \in 0,1,·\!·\!·,B/\alpha - 1,$ $m = 0,1,·\!·\!·,\alpha - 1$ ，则其DFT (Discrete Fourier Transform) $\hat X$ 满足： ${\hat X_i} = \displaystyle\sum\nolimits_{m = 0}^{\alpha - 1} {{X_{{f_{i,m}}}} }$ ，其中 ${X_{{f_{i,m}}}}$ 表示原始频谱中混叠在第 $i$ 个频点位置上的第 $m$ 个频点谱值。假设采样通道个数为 $C$ ，每个通道的延时为 ${\tau _s},s = 0,1,·\!·\!·,C - 1$ ，则

其中， ${w^{{f_{i,l}}{\tau _s}}} = {{\rm e}^{{\rm j}2{{π}} {f_{i,l}}{\tau _s}}}$ , $\hat X_i^{({\tau _s})}$ 表示延时为 ${\tau _s}$ 的通道的第 $i$ 个频点的谱值。

考虑频谱是稀疏的，设每个混叠频点上混叠有 $r(r < C)$ 个谱值非零的频点，则式(1)可以写为

因此，从欠采样得到的混叠频谱中恢复原始频谱的问题就转化为求频点位置 ${f} = {[{f_{i,1}},{f_{i,2}},·\!·\!·,{f_{i,r}}]^{\rm T}}$ 及其对应谱值 ${X} = {[{X_{{f_{i,1}}}},{X_{{f_{i,2}}}},·\!·\!·,{X_{{f_{i,r}}}}]^{\rm T}}$ 。BigBand中使用最小均方差的方法首先确定 ${f}$ ，此时式(2)就变成关于 ${X}$ 的一个超定线性方程组。注意到，问题式(2)的未知数个数为 $2(C - 1)$ ，而已知的样本数量为 $C$ ，因此综合来说，该问题是个欠定的非线性问题，也因此导致了该算法对噪声的敏感度较高。下面介绍本文提出的一种通过构建频率位置多项式的方法来解决文献[9]中的问题。

3.   基于频率位置多项式的混叠频谱解算算法

3.1   频率位置多项式

首先给出以下定义和定理。对第 $i$ 个混叠频点，定义： ${B_i} = \{ f|0 \le f < B,f \equiv i\,{\rm mod} B/\alpha \}$ 。

定义 1　令 ${L_i} = \{ {f_{i,1}},{f_{i,2}}, ·\!·\!· ,{f_{i,r}}\} = \{ f \in {B_i}|$ ${X_f} \ne 0\}$ 为 ${B_i}$ 中混叠有非零频点的位置。定义

为频率位置多项式(Frequency Locator Polynomial, FLP)，其中 $\tau$ 为延时。

定理 1　令 ${\tau _i} = i/B,i = 0,1,·\!·\!·,C - 1$ 表示通道 $i$ 的延时。假设混叠频谱中每个频点最多混叠有 $r,r < C$ 个谱值非零的频点，并令 $G(z)$ 为频率位置多项式，使得

则有如式(4)的线性方程成立：

其中， $i = 0,1,·\!·\!·,B/\alpha - 1$ , ${a_0} = 1$ 。

证明　令

令 ${\omega _s} \!=\! {{\rm e}^{{\rm j}2{{π}} {f_{i,s}}\tau }},s \!=\! 1,2, ·\!·\!· ,r$ ，则 ${\omega _s},s = 1,2, ·\!·\!· ,r$ 为多项式 $G(z)$ 的根。考虑 $r = C - 1$ 情形，式(2)可写为

令 ${{a}}^{(t)} \!\!=\!\! \left[ {\overbrace {0, ·\!\!·\!\!· ,0}^t,{a_0}, ·\!\!·\!\!· ,{a_r},\overbrace {0, ·\!\!·\!\!· ,0}^{C - 1 - t - r}} \right],t \!=\! 0,1,·\!\!·\!\!·,$ $C - r - 1$ 。则有

因此，

证毕

由式(5)可知，原始频谱中的谱值非零的频点位置信息包含在多项式 $G(z)$ 中，因此如果可以确定 $G(z)$ ，即FLP，则可以确定 ${f}$ ，进而利用式(2)求解得到 ${X}$ 。而事实上，由于噪声存在，所得到的FLP系数与真实值存在偏差，因此需要对 $G(z)$ 进行估计。下面提供以下估计方法。

3.2   频率位置多项式的估计

令各信号采样通道的延时满足 ${\tau _i} = i/B,i = 0,$ $1,·\!·\!·,C \!- \!1$ ，假设非零频点的位置在连续 $d(d \ge C \!-\! 1)$ 个时间窗内保持不变，且 $d(C - r) \ge r$ 。令 $\hat X_{i,t}^{({\tau _s})}$ 表示第 $t$ 个时间窗的 $\hat X_i^{({\tau _s})}$ 值，则式(2)可以扩展成如式(9)方程组。

其中， ${{c}} = {[{c_1},{c_2}, ·\!·\!·, {c_r} ]^{\rm{T}}}$ 为FLP的系数 ${[{a_1},{a_2}, ·\!·\!·, {a_r} ]^{\rm{T}}}$ 的估计值。式(9)的构造过程如图1所示，图中的时间轴单位是 $N/{f_{\rm nyq}}$ ， $N$ 为FFT点数， ${f_{\rm nyq}}$ 为Nyquist采样率。

可见式(9)是一个关于 ${c}$ 超定线性方程组，取 $d$ 使得 $d(C - r) \ge r$ 且 ${{ {X}}_2}$ 列满秩，于是有：

(1) 无噪声的情形下，利用类似BCH译码中搜索错位多项式根的方法得到精确的FLP：

当 $d = 1$ , $C - r \ge r$ ，则利用Berlekamp-Massey算法可得到式(9)唯一的非平凡解，该解唯一确定FLP。

图  1  频率位置多项式系数估计方程组构造过程

下载: 全尺寸图片 幻灯片

当 $d > 1$ ，上述问题就变成了多序列综合问题。由定理1以及列满秩的条件可得， ${[{a_1},{a_2}, ·\!·\!· ,{a_r} ]^{\rm{T}}}$ 是式(9)的唯一解，即 ${c} = {[{a_1} ,{a_2}, ·\!·\!·,{a_r} ]^{\rm{T}}}$ ，该解可以通过文献[12]中的算法快速得到。

(2) 有噪声情形下，注意到 ${{{X}}_2}$ 列满秩，可利用最小二乘法求解得到 ${c}$ ，即 ${c} = {({{X}}_2^{\rm{H}}{{X}}_2})^{ - 1}{{X}}_2^{\rm{H}}{{{X}}_1}$ 。

需要注意以下几点：

(a) 在 $d(C - r) \ge r$ 条件下， ${{{X}}_2}$ 列满秩的概率很大。

(b) 由于在实际情况下，可使用的采集通道数很有限，而文献[9]中又证明了3个欠采样通道就足够，即 $C = 3$ ，因此定理中假设非零频点位置保持 $T\,$ 个时间窗不变在实际频谱感知中是合理的。

(c) $c(z)$ 与 $G(z)$ 的区别：在无噪声情况下，二者是相等的；而在有噪声的情况下，由复数多项式幅值的连续性可知， $G(z)$ 的每个根 ${z_r}$ 落在 $c(z)$ 的某一个根的很小的领域内，也就是说 $\left| {c({z_r})} \right| \approx 0$ 。在下一节，将介绍一种快速搜索 $G(z)$ 多项式根的方法。

综上，确定原始频谱中谱值非零的频点位置就可以转换成求解一系列的线性方程组得到FLP和对应的根。

3.3   频率位置多项式根快速搜索

本节给出一种简单而又快速的FLP根的确定方法。由3.1节和3.2节可知，多项式FLP的根是与非零频点的位置唯一对应的，并且已知该位置的可能取值，即 $\{ {f_{i,1}},{f_{i,2}}, ·\!·\!·,{f_{i,r}}\} \!\subseteq\! {B_i},i = 0,1,·\!·\!·,B/\alpha - 1$ 。称 $\{ {w^{{f_{i,s}}\tau }}|{f_{i,s}} = i + (s - 1)B,s = 1,2, ·\!·\!· ,\alpha \}$ 为第 $i$ 个混叠频点对应的FLP的侯选根集合， ${w^{{f_{i,s}}\tau }},$ $s = 1,2, ·\!·\!· ,\alpha$ 为 $c(z)$ 的候选根。

无噪声条件下，可直接通过精确计算来确定FLP的根。

有噪声条件下，选取候选根集合中 $r$ 个使得 $\left| {c(z)} \right|$ 最小的元素作为FLP的根。这样处理的依据在3.2节的注意里面已经提到。

通过上述方法，只需要将候选的数个根代入估计得到的FLP，从而确定非零频点的位置。

3.4   稀疏混叠频谱恢复算法

不难发现，本文算法可以分为两步完成。

步骤 1　确定非零频点 ${f}$ 。由定理1可知，频率位置多项式是存在的，并且是由非零频点位置唯一确定的；由3.2节的方法又可得到该多项式的估计；又由3.3节给出的方法搜索FLP的根，从而确定非零频点位置。此过程可认为是解式(9)得到FLP，然后从侯选根集合中找出FLP的根，进而换算对应频点位置。

步骤 2　计算谱值 ${X}$ 。步骤1中确定了频点位置，则式(6)就变成关于 ${X}$ 的超定线性方程组，因此可以利用最小二乘法解出 ${X} = {ω} \backslash \hat {{X}}$ 。

3.5   算法复杂度分析

BigBand提供了两种采样模式：一种是互素采样(co-prime sampling)，一种是延时错位采样(shifted sampling)。在算法实现过程中使用的是后者，本文从实用角度出发，此处只针对后者做比较。第2种方法在频谱恢复中采用的是一种遍历的搜索方法，对每个混叠频点求解式(6)次数为 $C_\alpha ^{C - 1}$ ，由于求解式(6)的计算复杂度为 $O({C^3})$ ，因此文献算法总的计算复杂度为 $O(C_\alpha ^{C - 1}{C^3}B/\alpha )$ 。

对于本文算法，需要分别求解式(6)和式(9)各一次，因此总的耗时为 $O(2{C^3}B/\alpha )$ 。

由上述分析可知，本文算法计算速度与BigBand的比值是 $C_\alpha \!\!^{C - 1}/2$ 。也就是说，当 $C \!=\! 3,\alpha \!=\! 18$ (文献[9]实验中使用)，这个比值为85.5；而当 $C = 4,\alpha = 18$ ，比值为484.5。可见，本文算法在计算复杂度上具有较大优势。

3.6   同类算法比较

(1)与压缩感知比较：本文的求解框架和压缩感知整体上是一致的，即Support Recovery(对应本文非零频点定位)和谱值计算(两类算法没有区别)。谱支撑域恢复是一个NP-hard问题，计算负担较大，一般压缩感知使用的是一种贪婪算法。其中OMP(Orthogonal Matching Pursuit)使用较多^[13]，计算复杂度为 $O\left( {{K^2}B\ln B} \right)$ ，其中 $K$ 为频谱稀疏度，针对本文的情形，可支持的稀疏度最大值为 $K = (C - 1)B/\alpha$ 。

注意到本文的多子带信号模型，稀疏字典具有傅里叶变换结构；另外由于采样率大于或等于子信号最大带宽，且子信号个数最大值( $C - 1$ )已知，则频域里每个混叠频点上最多混叠有 $C - 1$ 个信号频率。利用以上两点信息，可以使用结构化的稀疏表示模型来解决本文的问题。事实上，文献[9]中的BigBand算法使用的就是这一方法来确定信号频点位置。

此外，在硬件实现上，如文献[9]中叙述，本文算法(和BigBand)可以在易得商用的标准采集卡上实现；而压缩感知则需要额外的模拟混频电路，不能使用常用的AD设备，因此硬件复杂度相对较高。

(2)与Prony方法比较：不难发现，本文算法与Prony模型类似。Prony方法通过求解特征方程确定频谱中非零频点位置，继而求得谱值。但两种方法有以下几个方面的区别：

第一，Prony方法对噪声很敏感，原始的Prony方法是不适用于处理含噪信号的，经过数百年的发展，适用于有噪声情形的信号频率成分分析的改进Prony方法被提出。另一方面，本文算法则适用于有噪声的频谱恢复。

第二，二者的应用场合不同。Prony方法适合处理频率成分较简单的信号感知，比如多音信号，而对信号带宽较大的通信信号不适用。本文算法是为宽带频谱感知设计，支持更多类的信号，包括多音信号。

第三，Prony模型一般直接对时域信号进行处理，虽然只需要求解一个特征方程，但随着正弦波个数的增多，在有噪声的情况下，性能会大为降低。而本文算法，在频域进行处理，针对每个混叠频点分别进行解算，有效地提高了频率分辨能力，即提高了算法所支持的频谱稀疏度。

考虑到本文算法模型与Prony模型类似，也可以将本文算法看作是Prony模型在欠采样信号频域处理的一种推广应用。

图  2  本文原始频谱、欠采样频谱及恢复的频谱

下载: 全尺寸图片 幻灯片

4.   实验与结果分析

设置 $B$ =100 MHz，采样通道数为 $C = 4$ , $d = r = 3$ ，欠采样倍数为 $\alpha = 10$ ，即每个通道的采样率为 $B/\alpha = 10\; {\rm{MS/s}}$ 。实验中采用一块四通道采集卡(分辨率16 bit，采用率100 MS/s，模拟输入带宽100 M)，输入信号如表1所示，信号分别用3台Agilent E4438C生成，并经由合路器与功分器输入至采集卡。

表  1  原始频谱中的信号

信号	调制类型	载频(MHz)	波特率(MHz)
A	QPSK	23	2
B	QPSK	33	2
C	QPSK	43	2

下载: 导出CSV 

| 显示表格

4.1   频谱恢复

图2和图3分别给出了本文算法和BigBand^[9]对经过10倍欠采样后的混叠频谱进行恢复的结果图，其中图2(c)和图3的频谱恢复错误分别为 $6.9 \times$ ${10^{ - 6}}$ 和0.1162，此处的错误率的定义是：(虚判频点数+错判频点数)/(总带宽)。可见本文算法的准确度较文献算法具有较大的优势，在 $C = 4$ , $r = 3$ 的条件下，文献[9]算法性能基本不能接受。通过一定的仿真发现，对 $r = 3$ ，文献[9]算法需要的采样通道数至少为 $C = 6$ 才能达到本文算法的同样准确度，这说明，在相同采集通道数条件下，本文算法性能要高于文献算法。

造成这种区别的主要原因在于：BigBand在每个混叠频点解算中可用的样本数只有 $C$ 个，而待定参数 ${f}$ 与 ${X}$ 总个数是 $2(C - 1)$ ，因此该解算过程是一种欠定非线性问题。而本文算法利用信号频谱位置在短时间内不变的，即频率位置多项式是不变的，通过构造式(9)，对每个混叠频点分别解算，可使用的样本数为 $(d + 1)C \ge {C^2}$ ( $d \ge C - 1$ )。并且将非线性问题式(6)通过构造频率位置多项式转换成两个超定的线性问题，这样，就导致本文算法较BigBand性能要提高较多。

4.2   性能分析

由3.2节可以看出，本文算法的性能与候选根的FLP赋值有密切关系。首先对 $\left| {c({w^{{f_{i,s}}\tau }})} \right|,$ $s \!=\! 1,2, ·\!·\!· ,\alpha$ 进行升序排列得到 $[{v_1}, {v_2}, ·\!·\!· , {v_\alpha }]$ ，显然 ${v_1} , ·\!·\!· , {v_r}$ 所对应的是FLP的根，即 ${v_1} \approx ·\!·\!· \approx {v_r} \approx 0 < < {v_{r + 1}} ·\!··$ 。定义算法对候选根的区分能力： ${v_{r + 1}}/{v_r}$ ，可见该值越大，说明估计得到的FLP越接近真实值。也就是说， ${v_{r + 1}}/{v_r}$ 表示了本文所提出的根搜索方法从候选根集合中区分真实根和非真实根的能力。

图  3  文献[9]恢复的频谱

下载: 全尺寸图片 幻灯片

图  4  侯选根的多项式幅值分布

下载: 全尺寸图片 幻灯片

为验证本文给出的根搜索方法的有效性，在4.1节中所搭建的实验平台上，通过多次试验，设置不同的信号频谱位置，记录每一次的升序排列后的侯选根的FLP赋值，然后取平均值(100次)，结果如图4所示。可以看出， ${v_3}$ 与 ${v_4}$ 之间有明显的差值，这是由于实验中 $r = 3$ 。注意到随着信噪比的降低， ${v_4}/{v_3}$ 的值在降低，但即使在信噪比为10 dB的条件下，这个值仍然足以将真实根从侯选根集合中区分出来。上述结果也从侧面说明了本文算法的可行性。

图5给出了不同信噪比条件下的混叠频谱恢复错误率(虚假谱线出现概率)，从图中可以看出，本文算法对噪声的鲁棒性明显强于文献[9]算法(BigBand)。在信噪比较高时(30 dB)平均错误率为 $1 \times {10^{ - 6}}$ 量级，即使在SNR=10 dB，该值也能保持在0.1量级。

为进一步测试算法的重构性能，图6中给出了解调误码率随信噪比的变化曲线。实验中使用的信号与表1一致，信噪比的定义是：信号总功率/宽带噪声功率。由前面的结果可以看出，针对本文所处理的3个信号混叠，四通道采样下的BigBand性能恶化严重，因而此处没有进一步列出其误码率，而只是给出了Nyquist采样下的解调误码率和本文算法在每个通道以10倍欠采样下的误码率性能比较。由图6中结果可以看出，本文算法相比Nyquist采样下的误码率有接近10 dB的信噪比差距，该结果和图2中给出的结果是吻合的，图2中恢复的信号相对原始信号有10 dB的信噪比损失，这也是下一步工作中需要重点关注和解决的问题。

图  5  不同信噪比下算法错误率比较

下载: 全尺寸图片 幻灯片

图  6  不同信噪比下的解调误码率

下载: 全尺寸图片 幻灯片

5.   结束语

本文在BigBand的基础上，通过重新设计稀疏频谱恢复算法，不仅使得原算法的计算复杂度得到有效降低，也使得算法在频谱恢复的准确度得到有力提升。本文的贡献的主要在于：提出了一种通过建立类似的错位多项式的频率位置多项式来快速定位谱值非零的频点的算法，从而准确计算其对应谱值；将稀疏频谱恢复的非线性问题转换成两种线性问题，通过求解一系列的线性方程组使问题得到解决。实验结果表明本文算法的可行性，所提出的频率位置多项式可以准确地表征非零频点的位置。因此，可以相信本文算法在实际的宽带频谱感知领域将有广阔的应用前景。

还有以下几点需要后续研究：首先，如何补偿欠采样所带来的信噪比损失；其次，关于鲁棒性的具体的数值分析需要研究；此外，算法可处理的信噪比理论性能界限需要推导。