CN102946288A - 基于自相关矩阵重构的压缩频谱感知方法 - Google Patents

Abstract

本发明提出了一种基于自相关矩阵重构的压缩频谱感知方法，主要解决现有感知算法采样速率高和计算开销过大的问题。其实现方法是：次级用户对频谱环境利用压缩感知获得观测序列，并利用观测序列的自相关矩阵重构出奈奎斯特采样序列的自相关向量，以此得到奈奎斯特采样的自相关矩阵估计值；之后采用多信号分类MUSIC算法，根据自相关矩阵估计值的特征值得到被占用信道数目的估计值；根据特征值和占用信道数目的估计值构造特征谱，并根据各个信道上特征对应谱幅值相加所得和值，判断出各个被占用信道的标号。本发明能够降低次级用户接收机的采样速率，重构端算法复杂度低，可用于认知无线电系统中快速判断频谱占用情况。

Description

基于自相关矩阵重构的压缩频谱感知方法

技术领域：

本发明属于通信技术领域，涉及一种频谱感知方法，更进一步涉及一种基于自相关矩阵重构的压缩频谱感知方法，可用于认知无线电系统。

背景技术：

近几年来，无线通信的快速增长导致授权频带和非授权频带上的无线服务需求都急剧增长。然而，目前的固定频谱分配策略使得频谱使用效率很低。为了解决频谱资源短缺和频谱利用率不高的问题，人们提出了认知无线电CR的概念。CR技术能够有效地利用未占用频带，进而提高频谱利用率。在CR系统中，每个次级用户SU最重要的任务之一就是频谱感知，即对射频环境进行检测，以判断在感兴趣的信道或频段上是否有授权用户LU在进行通信。频谱感知的主要目的在于高效地识别未利用的频段或物理信道，即频谱空洞，以供CR使用，从而提高整个认知系统的吞吐量和服务质量并保护LU的通信不受干扰。

一些成熟的、得到广泛应用的频谱感知技术，如能量检测，匹配滤波器检测和循环平稳特征检测等技术，均是以窄带感知为目的提出的方法。而未来的认知无线电应该有扫描宽带频段的能力，而这个带宽可能宽至数GHz。SU需要对整个带宽上LU的功率谱进行估计，并自适应地改变工作频段以进行传输。这就需要SU的射频前端使用一系列可调谐的窄带带通滤波器，针对每个子带逐次进行检测；或是一次性地，利用高采样率的数模转换器ADC和高性能的数字信号处理器DSP对整个带宽进行检测。这样两种检测策略均有实现上的困难存在：前者需要大量的射频参数，而且每个滤波器的调谐范围必须预先设定；后者则对ADC和DSP的处理速率，以及存储设备的容量有很高的要求。

压缩感知是一种以低于奈奎斯特采样速率获取稀疏性信号的采样方法，这一技术在频谱感知领域也得到了广泛的应用。由于目前固定频率分配策略导致的低效频谱利用率，大部分场景内SU探测到的信号在频率域上只有部分频段功率值很大，而其余的大部分频段功率较低，即在频域上是稀疏的，这种性质允许接收机对信号以低于奈奎斯特速率的欠采样方式获取采样点，之后通过一定的信号重构算法恢复原始信号。在利用压缩感知进行频谱感知时，SU对宽带信号进行压缩采样，并利用重构算法，将频域稀疏的原始信号从压缩采样值中还原出来，并以此为基础对信号进行判决，从而判断观测频段是否被LU占用。这样的方式虽然能有效降低采样速率和存储容量，但是当前的技术在重构模块仍有挑战：压缩重构需要从少量的采样点中求解出更多的数值，这相当于求解一个未知数远多于方程数的欠定线性方程组。现有的重构算法多是用迭代的方式求解一个凸优化问题，这种循环迭代的方式在采样点数较多时计算复杂度很高，很难适应快速变化的频谱环境。

在Signal Processing,IEEE Transactions,2012《Compressive Sampling for PowerSpectrum Estimation》一文中，作者Geert Leus提出一种通过重构信号自相关向量的方式实现的信号功率谱重构方法。该方法在利用感知矩阵对原始的模拟信号进行欠采样后，在重构端并不完全重构完整的原始信号，而是着眼于通过还原原始信号自相关向量这一信号统计特征的方式，来间接地重构原信号的功率谱。文章提出，对于宽平稳信号来说，其自相关矩阵是一个Toeplitz矩阵，这意味着这个自相关矩阵中的一列就已经包含了整个自相关矩阵的所有信息；同时，欠采样序列并不是一个宽平稳的过程，因此其自相关矩阵不具有Toeplitz矩阵的性质。文章根据这一特点，利用欠采样序列的自相关矩阵的所有列的信息，通过简单的最小二乘法，估计出原始信号自相关矩阵的其中一列，即奈奎斯特采样序列的自相关向量，并据此得到原始信号的功率谱估计。但是仍然存在一定的问题：在检测带宽非常大时，重构整个频段的功率谱仍然会带来计算开销过大的问题，同时，功率谱的估计误差在进行检测时也会带来检测性能的损失。

发明内容

本发明针对上述现有基于功率谱重构的压缩频谱感知方法的不足，提出一种基于自相关矩阵重构的压缩频谱感知方法，以在不降低检测概率的条件下，减小宽带压缩频谱感知的计算开销。

实现本发明的技术思路在于，利用欠采样序列的自相关矩阵重构出奈奎斯特采样序列的自相关向量，并根据以此得到奈奎斯特采样的自相关矩阵的估计值，之后采用多信号分类MUSIC频率估计算法，对整个频段进行分析，判断出各个被占用信道出现的位置。具体实现步骤包括如下：

（1）次级用户对频域稀疏的宽带模拟信号进行压缩采样：

（1a）次级用户利用数模转换器ADC对模拟信号进行奈奎斯特采样，得到N×1维采样序列x[k]，其中k取正整数，N表示模拟信号频段上的信道总数；

（1b）生成M×N维观测矩阵Φ，并利用观测矩阵Φ对采样序列x[k]进行压缩采样，得到M×1维观测序列y[k]，其中观测矩阵Φ是一个高斯随机矩阵， $\sqrt{2N-1}<M<N;$

（2）估计观测序列y[k]的自相关矩阵R_y，并对该自相关矩阵R_y做向量化处理，得到M²×1维向量vec(R_y)；

（3）利用步骤（1）中的观测矩阵Φ，构造重构矩阵Θ；

（4）估计采样序列x[k]的自相关矩阵R_x：

（4a）根据步骤（2）中的向量vec(R_y)和步骤（3）中的重构矩阵Θ，通过最小二乘法求解关于r_x的方程：vec(R_y)＝Θr_x，得到最小二乘解

作为自相关向量r_x的估计值，其中，r_x表示由采样序列x[k]的自相关矩阵R_x第一行和第一列的元素构成的(2N-1)×1维向量；

（4b）将上述

排列成Toeplitz矩阵得到

将作为自相关矩阵R_x的估计值；

（5）根据自相关矩阵R_x的估计值

判定授权用户占用的信道数目和标号：

（5a）对所述

进行特征值分解，得到降序排列的特征值λ₁≥λ₂≥...≥λ_N和对应的特征向量u₁,u₂,...,u_N；

（5b）根据上述N个特征值，找到其中明显大于其他特征值的K个特征值，将K作为授权用户占用信道数目的判定结果；

（5c）利用上述最小的N-K个特征值对应的特征向量，构造噪声子空间特征矩阵G＝[u_K+1,u_K+2,...,u_N]；

（5d）构造频率矩阵A:

其中L表示特征谱的分辨率，L＝QN，其中Q表示在每个信道上特征谱的幅值数目，

l＝1,2,...,L是归一化的角频率；

（5e）利用上述频率矩阵A和噪声子空间特征矩阵G，计算采样序列x[k]的特征谱；

（5f）将N个信道中每个信道对应的特征谱的幅值相加求和，选择和值中较大的K个，将这K个值对应的标号作为占用信道的标号。

本发明与现有技术相比的主要优点：

（1）由于本发明采用了压缩采样，相比传统的宽带CR频谱检测，本发明降低了采样速率，节省了存储空间。

（2）本发明避免了传统压缩感知求解凸优化问题的重构方法，利用欠采样序列的自相关矩阵和奈奎斯特采样序列的自相关向量之间的线性关系，通过求解一个最小二乘解进行重构，仅需完成矩阵求逆及乘法运算，降低了运算复杂度，从而降低了频谱感知的计算开销。

（3）本发明在检测过程中使用MUSIC频率估计算法，无需对信噪比进行估计，避免了先验信息误差带来的性能损失。

附图说明:

图1是本发明的实现流程图；

图2是本发明仿真实验中使用的宽带模拟信号的功率谱示意图；

图3是本发明在固定采样长度时不同采样比下的检测性能比较；

图4是本发明在不同采样长度条件下同一采样比的检测性能比较；

图5是本发明在不同采样长度下重构的自相关矩阵特征值和特征谱的均方误差图。

具体实施方式:

以下对本发明的原理及技术方案做进一步的描述：

参照图1，本发明的实现流程图如下：

步骤1，次级用户利用观测矩阵Φ对模拟信号x(t)进行压缩采样。

1.1）假设模拟信号x(t)的频率范围为[f_min,f_max]，这段频率上共有N个互不重叠的物理信道供授权用户LU使用，其标号分别为1,2,...,N，每个信道带宽为B，并且假设次级用户SU拥有各子带位置的先验信息。

将x(t)通过模数转换器ADC进行奈奎斯特采样，取其输出序列x[n]中的N个点构成N×1维采样序列x[k]：

x[k]＝[x[kN],x[kN+1],...,x[kN+N-1]]^T

其中T表示矩阵转置，k取正整数，N表示频段上的信道总数。

1.2）生成一个M×N维观测矩阵Φ，并利用观测矩阵Φ对采样序列x[k]进行压缩采样，得到观测序列y[k]：

其中

观测矩阵Φ是一个高斯随机矩阵，其每个元素服从均值为0、方差为1的高斯随机分布，

是观测矩阵Φ的第i行的行向量。

步骤2，估计观测序列y[k]的M×M维自相关矩阵R_y：

R_y＝E[y[k]y^H[k]]

其中H表示矩阵共轭转置；并将自相关矩阵R_y矩阵向量化，得到向量化结果vec(R_y)，即将R_y中的每一个M×1维的列向量，依次放入一个M²×1维的向量vec(R_y)中去。

步骤3，构造重构矩阵Θ：

3.1）对步骤1中观测矩阵

的各行向量进行相关运算，得到相关系数

其中

分别表示观测矩阵Φ第i行和第j行的行向量，

表示行向量

的第n+s个元素，

表示行向量

的第n个元素，s＝1-N，...0,...,N-1，N表示模拟信号频段上的信道总数，T表示转置运算，*表示复共轭运算；

3.2）利用上述相关系数求得相关向量

i,j＝0,1,...,M-1

3.3）利用上述相关向量构成重构矩阵Θ：

步骤4，利用步骤2中的向量vec(R_y)和步骤3中的重构矩阵Θ恢复采样序列x[k]的自相关矩阵R_x：

4.1）利用最小二乘方法，求解下面关于r_x的方程：

vec(R_y)＝Θr_x

得到最小二乘解

作为r_x的估计值，其中r_x表示R_x第一行和第一列的元素构成的(2N-1)×1维的向量，即：

r_x＝[[R_x]_1,N,...,[R_x]_1,2,[R_x]_1,1,[R_x]_2,1,...,[R_x]_N,1]^T

＝[r_x(1-N),...,r_x(-1),r_x(0),r_x(1),...,r_x(N-1)]^T

这里[R_x]_ij＝r_x(i-j)＝r_x ^*(j-i)＝E[x[n+i-j]x^*[n]]，表示矩阵R_x的第i行、第j列的元素。

4.2）将上述

排列成一个Toeplitz矩阵得到R_x，将R_x作为

的估计值，即：

其中，

表示

的第i个元素。

由于宽平稳信号的自相关矩阵R_x的元素满足[R_x]_ij＝r_x(i-j)＝r_x ^*(j-i)，即r_x包含了R_x中的所有信息，因此将排成Toeplitz矩阵即可得到

步骤5，根据上述

判定LU占用的信道数目和标号：

5.1）对

进行特征值分解，得到降序排列的特征值λ₁≥λ≥₂...≥λ_N，以及每个特征值对应的特征向量u₁,u₂,...,u_N。

5.2）利用信息论准则分析特征值，判断出明显大于其他特征值的K个特征值，即信号子空间对应的特征值，将K作为LU占用的信道数目的判定结果。

5.3）构造噪声子空间特征矩阵G＝[u_K+1,u_K+2,...,u_N]，这里u_K+1,u_K+2,...,u_N为最小的N-K个特征值对应的特征向量，即噪声子空间对应的特征向量。

5.4）构造频率矩阵A：

其中，L表示特征谱的分辨率，L越大，特征谱的频率分辨率越高，L可取值为信道数目N的整数倍，即：L＝QN，其中Q表示在每个信道上特征谱的幅值数目，

l＝1,2,...L是归一化的角频率。

5.5）根据上述频率矩阵A和噪声子空间特征矩阵G计算特征谱P_music(ω_l)：

$P_{\mathrm{music}}\left({\mathrm{\&omega;}}_l\right)=\frac{1}{a^H\left({\mathrm{\&omega;}}_l\right){\mathrm{GG}}^{H}a\left({\mathrm{\&omega;}}_l\right)}=\frac{1}{\underset{i=K+1}{\overset{N}{\mathrm{\&Sigma;}}}{\left|a^H\left({\mathrm{\&omega;}}_l\right)u_i\right|}^2},$ l＝1，2,...,L

其中H表示共轭转置，a(ω_l)表示频率矩阵A第l列的列向量。

5.6）根据上述特征谱P_music(ω_l)，利用下式对每个信道对应的特征谱的幅值相加求和，得到和值P(n)：

$P\left(n\right)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{P}_{\mathrm{music}}\left(\right[(n-1)Q+i\left]\frac{2\mathrm{\&pi;}}{L}\right),$ n＝1,2,...,N

并选择较大的K个P(n)值，将对应的标号n作为占用信道的标号。

本发明的效果可通过仿真进一步说明：

1）仿真参数：

仿真场景为检测LU用户在DVB-T 2K模式下，利用OFDM调制传输的信号，LU将随机使用20个物理信道中的10个，这些信道中心频率的分布范围为590~765MHz。其中每个物理信道的带宽为7.61MHz，子载波数为1705，码元周期为224μs。一个典型的功率谱实例如图2所示。

2）仿真内容与结果：

仿真1，在固定采样长度为1024点的条件下，对本发明与未采用压缩感知的检测方法的检测概率P_d和虚警概率P_f随信噪比变化进行1000次蒙特卡洛实验，得到对比仿真结果如图3所示。图3中八条曲线分别表示采样比为M/N＝0.4、MN＝0.5、MN＝0.6时本发明的检测性能以及M/N＝1即未压缩采样时的检测概率和虚警概率。从图3可以看出，本方法在信噪比高于-2dB时检测概率几乎为1，虚警概率始终保持在0.05以下，欠采样条件下的检测性能并未显著下降。

仿真2，在固定采样比M/N＝0.5的条件下，对不同的采样序列长度下本发明的检测性能随信噪比变化进行蒙特卡洛仿真，得到仿真结果如图4所示。图4中的六条曲线分别表示采样序列长度为1024点、2048点及4096点时本发明的检测概率和虚警概率。从图4可以看出，在同一欠采样比条件下，奈奎斯特采样序列长度越长，本发明的检测概率越高；因此对比图3、图4可以看出，影响本发明检测性能的主要因素是奈奎斯特采样点数，而采样比对方法性能的影响很小。因此，只要保证了原奈奎斯特采样的长度，可以在满足M²≥2N-1的条件下尽可能的降低采样比，以减少计算开销。

仿真3，在固定信噪比为0dB的条件下，对不同采样长度下本发明估计的特征谱和特征值的均方误差MSE随采样比的变化进行蒙特卡洛仿真，得到仿真图5，如图5所示。图5中横轴表示采样比，其变化范围为0.45~0.95，纵轴表示均方误差值。从图5可以看出，特征值重构的均方误差保持在10^-2以下，特征谱P_music(ω_l)的均方误差保持在10^-4以下，较低的均方误差保证了检测方法的性能。

Claims (4)

1.一种基于自相关矩阵重构的压缩频谱感知方法，包括如下步骤：

（1）次级用户对频域稀疏的宽带模拟信号进行压缩采样：

（1a）次级用户利用数模转换器ADC对模拟信号进行奈奎斯特采样，得到N×1维采样序列x[k]，其中k取正整数，N表示模拟信号频段上的信道总数；

（1b）生成M×N维观测矩阵Φ，并利用观测矩阵Φ对采样序列x[k]进行压缩采样，得到M×1维观测序列y[k]，其中观测矩阵Φ是一个高斯随机矩阵， $\sqrt{2N-1}<M<N;$

（2）估计观测序列y[k]的自相关矩阵R_y，并对该自相关矩阵R_y做向量化处理，得到M²×1维向量vec(R_y)；

（3）利用步骤（1）中的观测矩阵Φ，构造重构矩阵Θ；

（4）估计采样序列x[k]的自相关矩阵R_x：

（4a）根据步骤（2）中的向量vec(R_y)和步骤（3）中的重构矩阵Θ，通过最小二乘法求解关于r_x的方程：vec(R_y)＝Θr_x，得到最小二乘解

作为自相关向量r_x的估计值，其中，r_x表示由采样序列x[k]的自相关矩阵R_x第一行和第一列的元素构成的(2N-1)×1维向量；

（4b）将上述

排列成Toeplitz矩阵得到

将

作为自相关矩阵R_x的估计值；

（5）根据自相关矩阵R_x的估计值

判定授权用户占用的信道数目和标号：

（5a）对所述

进行特征值分解，得到降序排列的特征值λ₁≥λ₂≥...≥λ_N和对应的特征向量u₁,u₂,...,u_N；

（5b）根据上述N个特征值，找到其中明显大于其他特征值的K个特征值，将K作为授权用户占用信道数目的判定结果；

（5c）利用上述最小的N-K个特征值对应的特征向量，构造噪声子空间特征矩阵G＝[u_K+1,u_K+2,...,u_N]；

（5d）构造频率矩阵A:

其中L表示特征谱的分辨率，L＝QN，其中Q表示在每个信道上特征谱的幅值数目，

l＝1,2,...,L是归一化的角频率；

（5e）利用上述频率矩阵A和噪声子空间特征矩阵G，计算采样序列x[k]的特征谱；

（5f）将N个信道中每个信道对应的特征谱的幅值相加求和，选择和值中较大的K个，将这K个值对应的标号作为占用信道的标号。

2.根据权利要求1中所述的基于自相关矩阵重构的压缩频谱感知算法，其中步骤（2）所述的对观测序列自相关矩阵R_y做向量化处理，是指将观测序列自相关矩阵R_y中的每一个M×1维的列向量，依次放入一个M²×1维的向量vec(R_y)中去。

3.根据权利要求1中所述的基于自相关矩阵重构的压缩频谱感知算法，其中步骤（3）所述的利用步骤（1）中的观测矩阵Φ构造重构矩阵Θ，按如下步骤进行：

（3a）对步骤（1）中观测矩阵

的各行向量进行相关运算，得到相关系数

其中

分别表示观测矩阵Φ第i行和第j行的行向量，表示行向量

的第n+s个元素，

表示行向量

的第n个元素，s＝1-N,...0,...,N-1，N表示模拟信号频段上的信道总数，T表示转置运算，*表示复共轭运算；

（3b）利用上述相关系数求得相关向量

i,j＝0,1,...,M-1

（3c）利用上述相关向量构成重构矩阵Θ：

4.根据权利要求1中所述的基于自相关矩阵重构的压缩频谱感知算法，其中步骤（5e）所述的利用频率矩阵A和噪声子空间特征矩阵G，计算采样序列x[k]的特征谱P_music(ω₁)，按如下公式进行：

$P_{\mathrm{music}}\left({\mathrm{\&omega;}}_l\right)=\frac{1}{a^H\left({\mathrm{\&omega;}}_l\right){\mathrm{GG}}^{H}a\left({\mathrm{\&omega;}}_l\right)}=\frac{1}{\underset{i=K+1}{\overset{N}{\mathrm{\&Sigma;}}}{\left|a^H\left({\mathrm{\&omega;}}_l\right)u_i\right|}^2},$     l＝1，2,...,L

其中H表示共轭转置，a(ω1)表示频率矩阵A第l列的列向量，u_i表示第i个特征值对应的特征向量，N表示模拟信号频段上的信道总数，

l＝1,2,...,L是归一化的角频率，L表示特征谱的分辨率，K是授权用户占用信道数目的判定结果。

Images