CN106230441B - 一种基于m序列的可变维度的压缩感知观测矩阵构造方法 - Google Patents

Abstract

一种基于m序列的可变维度的压缩感知观测矩阵构造方法，本发明涉及基于m序列的可变维度的压缩感知观测矩阵构造方法。本发明的目的是为了解决现有压缩感知观测矩阵硬件存储空间大、应用范围受限、其信号重构性能和重构概率低以及硬件复杂度高的缺点。具体过程为：根据待采样信号长度N选择子矩阵维数；得到m序列优选对；分别产生第一组m序列和第二组m序列；得到A和B；合并为Φ₁；选取Φ₁中前N列构成Φ₂；根据实际所需要的观测量，随机生成S个子集Γ_i；计算互相关；得到最优子集；根据最优子集中元素，选取Φ₂对应的行序号，构成本发明的观测矩阵。本发明用于信号处理和通信技术领域。

Description

一种基于m序列的可变维度的压缩感知观测矩阵构造方法

技术领域

本发明涉及基于m序列的可变维度的压缩感知观测矩阵构造方法。

背景技术

信号处理和通信技术领域必然会涉及到采样过程，采样理论中公认的基础理论是奈奎斯特采样定理，其要求采样频率必须大于等于被采样信号最高频率的两倍。随着信号处理以及通信技术的不断发展，信号包含的载波频率越来越高，已经远远超过当前模数转换器(ADC)所具有的采样速率。而且，需要大量的硬件资源来存储和处理奈奎斯特采样得到的大量的采样数据。这使得宽带高频的模拟信号采样成为信号处理发展严峻的瓶颈。因此，对于包含超高频率的信号，当前的硬件采样设备无法满足需求。

而通信信号等这些人造的信号具有一定的结构和特征，即在频域内是具有稀疏特性的。针对稀疏的信号处理，压缩感知技术是近年来的一项革命性的技术，是一种有效的信号获取方式，能够在远低于奈奎斯特采样频率的情况下对模拟信号进行采样，获取非常少量的观测值，并且能够根据这些少量的观测值对原始信号进行恢复。

观测矩阵是压缩感知理论中实现数据采集的核心，观测矩阵的设计是一个重要的研究方向。为了保证信号在观测测量的过程中不丢失信息，观测矩阵的设计需要满足一定的性质。为了能够精确重建原始信号，约束等距性质(Restricted Isometry Property，RIP)是观测矩阵需要满足的充分条件。然而，证明观测矩阵满足RIP性质是一个组合问题，目前没有有效的方法能够在多项式时间内验证观测矩阵是否满足RIP性质。一种可行替代的方法是评价观测矩阵的互相关(Mutual coherence)。理论证明，观测矩阵的互相关值越小，观测矩阵更容易满足RIP性质。其中广泛使用的满足RIP性质的观测矩阵有高斯随机矩阵和伯努利观测矩阵，虽然其维数可以任意产生，性能良好，但是由于元素的随机性，需要占用较大的存储空间，在硬件实现的复杂度很高，限制了压缩感知在实际的应用；以部分傅里叶矩阵、部分哈达玛矩阵等为代表的部分正交矩阵和以Toeplitz矩阵为代表的结构化矩阵，硬件实现复杂度较随机观测矩阵大大降低，但是其观测矩阵维数固定，限制了在信号采集、欠采样等实际工程的应用范围，而且其对信号恢复的概率有待于进一步提升。如何设计一种易于硬件实现，矩阵维度可变且性能优异的观测矩阵，成为压缩感知在调制宽带转换器等欠采样设备以及工程应用的关键，同时具有一定的理论意义；

综上，现有技术存在的问题是：

1、现有随机观测矩阵随机性导致的存储空间大、硬件要求高的问题；

2、现有部分正交矩阵和结构化矩阵维数固定不可改变，限制了其应用范围；

3、使用随机观测矩阵、部分正交矩阵和结构化矩阵进行信号恢复，其信号重构性能和重构概率有待于进一步提升；

4、现有欠采样设备所需要的观测矩阵维数大，需要在硬件中有更多的通道数导致的硬件复杂度高的问题。

发明内容

本发明的目的是为了解决现有压缩感知观测矩阵硬件存储空间大、应用范围受限、其信号重构性能和重构概率低以及硬件复杂度高的缺点，而提出一种基于m序列的可变维度的压缩感知观测矩阵构造方法。

一种基于m序列的可变维度的压缩感知观测矩阵构造方法具体步骤为：

步骤一、根据信号处理和通信领域中实际的待采样信号长度N选择子矩阵维数P×P，满足N由待采样信号的实际长度来确定；

步骤二、根据步骤一得到的子矩阵维数P×P得到m序列优选对，P满足P＝2^r-1，其中r为移位寄存器的级数；

所述，m序列是最长线性移位寄存器序列的简称，是一种伪随机序列；

步骤三、利用步骤二得到的m序列优选对分别产生第一组m序列a＝{a₁ a₂ a₃ …a_P-1}和第二组m序列b＝{b₁ b₂ b₃ … b_P-1}；

其中，a为m序列优选对产生的第一组m序列；a₁、a₂、a₃、a_P-1为周期为P的m序列a包含的元素；b为m序列优选对产生的第二组m序列；b₁、b₂、b₃、b_P-1为周期为P的m序列b包含的元素；

步骤四、对步骤三得到的第一组m序列和第二组m序列循环位移，分别得到Toeblitz矩阵A和B；

步骤五、将步骤四得到的Toeblitz矩阵A和Toeblitz矩阵B合并为观测矩阵Φ₁，大小为P×2P；

步骤六、选取步骤五得到的观测矩阵Φ₁中前N列构成观测矩阵Φ₂；

步骤七、在步骤六得到的观测矩阵Φ₂的基础上，根据实际所需要的观测量M，随机生成S个子集满足|Γ_i|＝M，S取值范围为100≤S≤10000；M＜N；

步骤八、根据步骤七得到的每一个子集Γ_i，计算互相关μ_i；

步骤九、根据步骤八得到的互相关μ_i，得到最优子集Γ_opt；

步骤十、根据步骤九得到的最优子集Γ_opt中元素，选取步骤六得到的观测矩阵Φ₂对应的行序号，构成本发明的观测矩阵Θ。

本发明的有益效果为：

1、本发明方法构造出来的观测矩阵具有伪随机特性。由于观测矩阵中的每一个元素是伪随机码序列以及延时移位产生的，即能够保证观测矩阵中的各个元素的出现具有随机性，从而为精确恢复原始信号提供保证；

2、本发明方法中观测矩阵的元素为±1序列，可以在硬件电路中直接实现，有利于节省硬件的存储空间；

3、本发明方法中观测矩阵的是由一对m序列优选对构成，可以通过移位寄存器来实现，而且可以拓展到其他低互相关的码序列，例如RS码序列等；

4、本发明方法中观测矩阵的每一行是由m序列延时得到的，能够保证观测矩阵元素的正确性，易于硬件实现；

5、本发明方法中观测矩阵的维数可以自由选择，克服了部分正交矩阵以及Toeplitz矩阵维数固定不能改变的缺点，拓展了观测矩阵的应用范围；

6、本发明方法中观测矩阵具有一定的确定性，仅需要较少的输入参数即可得到所需的观测矩阵；

7、本发明方法中观测矩阵对模拟信号进行采集，对信号的重构概率高于随机观测矩阵，保证了实际应用时的可靠度和精度；

8、本发明方法中观测矩阵可以拓展到部分正交矩阵以及Toeplitz矩阵等以实现维数可变，从而提供了一种可变维度的观测矩阵的构造框架。

如图4a、图4b所示使用本发明的观测矩阵用于调制宽带转换器进行信号采集，恢复后的信号时域和频域波形示意图。得出，采用本发明构造的观测矩阵作为混合函数，在采样频率仅有奈奎斯特频率的四分之一的时候，可以几乎没有失真的恢复原始模拟信号，保证了对原始模拟信号恢复的精度。

图5是给出了本发明的观测矩阵和伯努利观测矩阵用于调制宽带转换器进行信号采集，模拟信号在不同的稀疏度下恢复成功概率的对比示意图。可以看到，在模拟信号的频带数量为6的时候，使用本发明的观测矩阵用于调制宽带转换器进行信号采集，依然能够近乎100％的恢复信号。当原始模拟信号的频带数量继续增大时，使用本发明的观测矩阵用于调制宽带转换器进行信号采集，模拟信号恢复的成功率依然高于使用伯努利观测矩阵的调制宽带转换器。本发明设计的观测矩阵与随机伯努利观测矩阵相比，占用的存储资源和硬件资源少，同时对于模拟信号的恢复，表现出了更好的性能。本发明方法中观测矩阵对模拟信号采集，保证了实际应用时的高恢复概率和可靠度。

图6和图7是本发明设计的观测矩阵和随机高斯观测矩阵用于信号观测，分别对高斯信号和0-1信号，在不同的稀疏度下恢复成功概率的对比示意图。观测矩阵的维数为50×195。对于高斯稀疏信号，在信号的稀疏度为40的时候，随机高斯观测矩阵已经不能保证对于信号的精确重构，使用本发明涉及的观测矩阵依然能够保证对于原始信号100％的重构成功率。对于0-1信号，信号的重构性能较高斯稀疏信号有所下降。然而，较随机高斯观测矩阵，使用本发明设计的观测矩阵，能够获得最大45％的信号重构概率增益。可以看到，本发明设计的观测矩阵和非常适合作为压缩感知的观测矩阵使用。

综合考虑上述仿真实验,本发明设计的观测矩阵性能优于随机高斯观测矩阵和随机伯努利矩阵，同时克服了部分正交矩阵以及Toeplitz矩阵等维数固定不能改变的缺点，拓展了观测矩阵的应用范围。对高斯和0-1稀疏信号进行采集，信号的恢复概率高于随机高斯矩阵；在实际的调制宽带转换器欠奈奎斯特采样系统中，模拟信号恢复的成功率高于使用伯努利观测矩阵的调制宽带转换器，而且节省了硬件存储资源。本发明适用于用于压缩感知观测矩阵的构造以及在调制宽带转换器等欠采样设备以提高信号恢复质量。本发明方法中观测矩阵可以拓展到部分正交矩阵以及Toeplitz矩阵等以实现维数可变，从而提供了一种可变维度的观测矩阵的构造框架。

附图说明

图1是本发明中的观测矩阵的构造流程图；

图2是调制宽带转换器的示意图，p₁(t)为第1个通道的混频函数，p_i(t)为第i个通道的混频函数，P_I(t)为第I个通道的混频函数，h(t)为低通滤波器，t＝nT_s为采样过程，t为采样时刻，n为离散的采样序号，T_s为采样时间间隔，y₁(t)为第1个通道的采样值，y_i(t)为第i个通道的采样值，y_I(t)为第I个通道的采样值，I为通道的数量，10≤I≤100；

图3a是原始模拟信号的时域波形示意图；

图3b是原始模拟信号的频域波形示意图；

图4a是使用本发明的观测矩阵用于调制宽带转换器进行信号采集，恢复后的信号时域波形示意图；

图4b是使用本发明的观测矩阵用于调制宽带转换器进行信号采集，恢复后的信号频域波形示意图；

图5是伯努利观测矩阵和本发明的观测矩阵用于调制宽带转换器进行信号采集，模拟信号在不同的稀疏度下恢复成功概率的对比示意图；

图6是随机观测矩阵和本发明的观测矩阵用于信号观测，高斯信号在不同的稀疏度下恢复成功概率的对比示意图；

图7是随机观测矩阵和本发明的观测矩阵用于信号观测，0-1信号在不同的稀疏度下恢复成功概率的对比示意图。

具体实施方式

具体实施方式一、结合图1至图6说明本发明具体实施方式，一种基于m序列的可变维度的压缩感知观测矩阵构造方法具体步骤为：

步骤一、根据信号处理和通信领域中实际的待采样信号长度N选择合适的子矩阵维数P×P，满足N是一个确定的值，由待采样信号的实际长度来确定(比如待采样信号的实际长度和N相等)；

步骤二、根据步骤一得到的子矩阵维数P×P得到m序列优选对，P满足P＝2^r-1，其中r为移位寄存器的级数；

所述，m序列是最长线性移位寄存器序列的简称，是一种广泛应用的伪随机序列；

步骤三、利用步骤二得到的m序列优选对分别产生第一组m序列a＝{a₁ a₂ a₃ …a_P-1}和第二组m序列b＝{b₁ b₂ b₃ … b_P-1}；

其中，a为m序列优选对产生的第一组m序列；a₁、a₂、a₃、a_P-1为周期为P的m序列a包含的元素；b为m序列优选对产生的第二组m序列；b₁、b₂、b₃、b_P-1为周期为P的m序列b包含的元素；

步骤四、对步骤三得到的第一组m序列和第二组m序列循环位移，分别得到Toeblitz矩阵A和B；

步骤五、将步骤四得到的Toeblitz矩阵A和Toeblitz矩阵B合并为观测矩阵Φ₁，大小为P×2P；

步骤六、选取步骤五得到的观测矩阵Φ₁中前N列构成观测矩阵Φ₂；

步骤七、在步骤六得到的观测矩阵Φ₂的基础上，根据实际所需要的观测量M，随机生成S个子集满足|Γ_i|＝M，考虑到准确性和计算复杂度S取值范围为100≤S≤10000；M＜N；

步骤八、根据步骤七得到的每一个子集Γ_i，计算互相关μ_i；

步骤九、根据步骤八得到的互相关μ_i，得到最优子集Γ_opt(找到最小的μ_i对应的i即为最优值opt)；

步骤十、根据步骤九得到的最优子集Γ_opt中元素，选取步骤六得到的观测矩阵Φ₂对应的行序号，构成本发明的观测矩阵Θ。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤二中根据步骤一得到的子矩阵维数P×P得到m序列优选对，P满足P＝2^r-1，其中r为移位寄存器的级数；具体过程为：

m序列优选对的寻找方法为：若a是2^r阶有限域GF(2)的一个本原元，f₁(x)与f_t(x)是2^r阶有限域GF(2)上的r次本原多项式，a是f₁(x)的首根，取

t为最大互相关函数绝对值；

使a^t为r次本原多项式f_t(x)的一个根，则以r次本原多项式f₁(x)与f_t(x)为特征多项式所产生的m序列构成m序列优选对。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤四中对步骤三得到的第一组m序列和第二组m序列循环位移，分别得到Toeblitz矩阵A和B；具体过程为：

Toeblitz矩阵A和B为：

式中，a₀、a₁、a₂、a₃、a₄、…、a_P-2、a_P-1为周期为P的m序列a包含的元素；

b₀、b₁、b₂、b₃、…、b_P-2、b_P-1为周期为P的m序列b包含的元素。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤五中将步骤四得到的Toeblitz矩阵A和Toeblitz矩阵B合并为观测矩阵Φ₁，大小为P×2P；具体过程为：

式中，为观测矩阵Φ₁中列向量，即

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤六中选取步骤五得到的观测矩阵Φ₁中前N列构成观测矩阵Φ₂；具体过程为：

式中，为观测矩阵Φ₂中列向量，即

式中，a_N-1、a_N、a_N+1、a_N-2为周期为P的m序列a包含的元素；b_N-P-1、b_N-P、b_N-P+1、b_(N-2)modP为周期为P的m序列b包含的元素。

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：所述步骤七中在步骤六得到的观测矩阵Φ₂的基础上，根据实际所需要的观测量M，随机生成S个子集满足|Γ_i|＝M，考虑到准确性和计算复杂度S取值范围为100≤S≤10000；M＜N；具体过程为：

随机地在集合{1,2,…,N}中选择S个子集合Γ_i，S个子集合Γ_i中元素的个数为M个。

其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：所述步骤八中根据步骤七得到的每一个子集Γ_i，计算互相关μ_i；具体过程为：

式中，i＝1,2,…,S，1≤k≤N，1≤l≤N且k≠l；为选取列向量中序号在集合Γ_i中的元素组成的列向量，为选取列向量中序号在集合Γ_i中的元素组成的列向量。

其它步骤及参数与具体实施方式一至六之一相同。

具体实施方式八：本实施方式与具体实施方式一至七之一不同的是：所述步骤九中根据步骤八得到的互相关μ_i，得到最优子集Γ_opt；具体过程为：

根据步骤八得到的互相关μ_i，通过计算opt＝argmin{μ_i},i＝1,2,…,S，得到最优子集Γ_opt(找到最小的μ_i对应的i即为最优值opt)。

其它步骤及参数与具体实施方式一至七之一相同。

采用以下实施例验证本发明的有益效果：

实施例一：

本实施例一种基于m序列的可变维度的压缩感知观测矩阵构造方法具体是按照以下步骤制备的：

以欠奈奎斯特采样系统宽带调制转换器为例，其框图如图2所示。待采样信号经过I个并行通道，信号在每个通道分别和调制序列相乘，经过低通滤波器，最后以低速率进行采样。而本发明构成的观测矩阵的每一行对应于每一个并行通道，从而实现信号的采集过程。图3a和图3b所示的原始模拟信号的函数表达式为：

其中，能量系数E_n和时延τ_n随机设定，信号带宽B_n＝50MHz，载波频率f_n随机分布在[0,5]GHz范围内，则信号的奈奎斯特采样频率至少为f_N＝10GHz。设置调制宽带转换器的并行通道数量为50，原始模拟信号分为195个等效的频谱块，从而设置每个通道内低通滤波器的截止带宽和低速模数转换器的采样频率f_p＝f_s＝f_NYQ/195≈51.3MHz，f_NYQ为信号的奈奎斯特频率。所以调制宽带转换器所需要的总采样率为50×51.3≈2.565GHz。此时调制宽带转换器所需要的观测矩阵维数为50×195。此时选取m序列周期为P＝127＝2⁷-1，即需要7级移位寄存器。选取7阶m序列优选对，对应的生成多项式为x⁷+x³+x²+x+1和x⁷+x³+1。根据本发明构造出观测矩阵来进行信号的欠奈奎斯特采样。根据本发明的构造步骤，本发明设计的观测矩阵则只需要两对循环移位寄存器来生成m序列，其他通道的混合序列根据已生成的m序列循环移位即可得到。与传统的随机伯努利观测矩阵相比，本发明构造的观测矩阵可以大大减少所需要的存储空间，易于硬件实现。

图4a、图4b给出了使用本发明的观测矩阵用于调制宽带转换器进行信号采集，恢复后的信号时域和频域波形示意图。可以看到，采用本发明构造的观测矩阵作为混合函数，在采样频率仅有奈奎斯特频率的四分之一的时候，可以几乎没有失真的恢复原始模拟信号，保证了对原始模拟信号恢复的精度。

图5是给出了本发明的观测矩阵和伯努利观测矩阵用于调制宽带转换器进行信号采集，模拟信号在不同的稀疏度下恢复成功概率的对比示意图。可以看到，在模拟信号的频带数量为6的时候，使用本发明的观测矩阵用于调制宽带转换器进行信号采集，依然能够近乎100％的恢复信号。当原始模拟信号的频带数量继续增大时，使用本发明的观测矩阵用于调制宽带转换器进行信号采集，模拟信号恢复的成功率依然高于使用伯努利观测矩阵的调制宽带转换器。本发明设计的观测矩阵与随机伯努利观测矩阵相比，占用的存储资源和硬件资源少，同时对于模拟信号的恢复，表现出了更好的性能。本发明方法中观测矩阵对模拟信号采集，保证了实际应用时的高恢复概率和可靠度。

图6和图7是本发明设计的观测矩阵和随机高斯观测矩阵用于信号观测，分别对高斯信号和0-1信号，在不同的稀疏度下恢复成功概率的对比示意图。观测矩阵的维数为50×195。对于高斯稀疏信号，在信号的稀疏度为40的时候，随机高斯观测矩阵已经不能保证对于信号的精确重构，使用本发明涉及的观测矩阵依然能够保证对于原始信号100％的重构成功率。对于0-1信号，信号的重构性能较高斯稀疏信号有所下降。然而，较随机高斯观测矩阵，使用本发明设计的观测矩阵，能够获得最大45％的信号重构概率增益。可以看到，本发明设计的观测矩阵和非常适合作为压缩感知的观测矩阵使用。

综合考虑上述仿真实验,本发明设计的观测矩阵性能优于随机高斯观测矩阵和随机伯努利矩阵，同时克服了部分正交矩阵以及Toeplitz矩阵等维数固定不能改变的缺点，拓展了观测矩阵的应用范围。对高斯和0-1稀疏信号进行采集，信号的恢复概率高于随机高斯矩阵；在实际的调制宽带转换器欠奈奎斯特采样系统中，模拟信号恢复的成功率高于使用伯努利观测矩阵的调制宽带转换器，而且节省了硬件存储资源。本发明适用于用于压缩感知观测矩阵的构造以及在调制宽带转换器等欠采样设备以提高信号恢复质量。本发明方法中观测矩阵可以拓展到部分正交矩阵以及Toeplitz矩阵等以实现维数可变，从而提供了一种可变维度的观测矩阵的构造框架。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (1)

1.一种基于m序列的可变维度的压缩感知观测矩阵构造方法，其特征在于：一种基于m序列的可变维度的压缩感知观测矩阵构造方法具体步骤为：

步骤一、根据信号处理和通信领域中实际的待采样信号长度N选择子矩阵维数P×P，满足

步骤二、根据步骤一得到的子矩阵维数P×P得到m序列优选对，P满足P＝2^r-1，其中r为移位寄存器的级数；

所述，m序列是最长线性移位寄存器序列的简称，是一种伪随机序列；

步骤三、利用步骤二得到的m序列优选对分别产生第一组m序列a＝{a₁ a₂ a₃ … a_P-1}和第二组m序列b＝{b₁ b₂ b₃ … b_P-1}；

其中，a为m序列优选对产生的第一组m序列；a₁、a₂、a₃、a_P-1为周期为P的m序列a包含的元素；b为m序列优选对产生的第二组m序列；b₁、b₂、b₃、b_P-1为周期为P的m序列b包含的元素；

步骤四、对步骤三得到的第一组m序列和第二组m序列循环位移，分别得到Toeblitz矩阵A和B；

步骤五、将步骤四得到的Toeblitz矩阵A和Toeblitz矩阵B合并为观测矩阵Φ₁，大小为P×2P；

步骤六、选取步骤五得到的观测矩阵Φ₁中前N列构成观测矩阵Φ₂；

步骤七、在步骤六得到的观测矩阵Φ₂的基础上，根据实际所需要的观测量M，随机生成S个子集满足|Γ_i|＝M，S取值范围为100≤S≤10000；M＜N；

步骤八、根据步骤七得到的S个子集中每一个子集Γ_i，计算互相关μ_i；

步骤九、根据步骤八得到的互相关μ_i，得到最优子集Γ_opt；

步骤十、根据步骤九得到的最优子集Γ_opt中元素，选取步骤六得到的观测矩阵Φ₂对应的行序号，构成观测矩阵Θ；

所述步骤二中根据步骤一得到的子矩阵维数P×P得到m序列优选对，P满足P＝2^r-1，其中r为移位寄存器的级数；具体过程为：

m序列优选对的寻找方法为：若a是2^r阶有限域GF(2)的一个本原元，f₁(x)与f_t(x)是2^r阶有限域GF(2)上的r次本原多项式，a是f₁(x)的首根，取

t为最大互相关函数绝对值；

使a^t为r次本原多项式f_t(x)的一个根，则以r次本原多项式f₁(x)与f_t(x)为特征多项式所产生的m序列构成m序列优选对；

所述步骤四中对步骤三得到的第一组m序列和第二组m序列循环位移，分别得到Toeblitz矩阵A和B；具体过程为：

Toeblitz矩阵A和B为：

式中，a₀、a₁、a₂、a₃、a₄、…、a_P-2、a_P-1为周期为P的m序列a包含的元素；

b₀、b₁、b₂、b₃、…、b_P-2、b_P-1为周期为P的m序列b包含的元素；

所述步骤五中将步骤四得到的Toeblitz矩阵A和Toeblitz矩阵B合并为观测矩阵Φ₁，大小为P×2P；具体过程为：

式中，为观测矩阵Φ₁中列向量，即

所述步骤六中选取步骤五得到的观测矩阵Φ₁中前N列构成观测矩阵Φ₂；具体过程为：

式中，为观测矩阵Φ₂中列向量，即1＜N≤P；P＜N≤2P；

式中，a_N-1、a_N、a_N+1、a_N-2为周期为P的m序列a包含的元素；b_N-P-1、b_N-P、b_N-P+1、b_(N-2)modP为周期为P的m序列b包含的元素；

所述步骤七中在步骤六得到的观测矩阵Φ₂的基础上，根据实际所需要的观测量M，随机生成S个子集满足|Γ_i|＝M，S取值范围为100≤S≤10000；M＜N；具体过程为：

随机地在集合{1,2,…,N}中选择S个子集合Γ_i，S个子集合Γ_i中元素的个数为M个；

所述步骤八中根据步骤七得到的每一个子集Γ_i，计算互相关μ_i；具体过程为：

式中，i＝1,2,…,S，1≤k≤N，1≤l≤N且k≠l；为选取列向量中序号在集合Γ_i中的元素组成的列向量，为选取列向量中序号在集合Γ_i中的元素组成的列向量；

所述步骤九中根据步骤八得到的互相关μ_i，得到最优子集Γ_opt；具体过程为：

根据步骤八得到的互相关μ_i，通过计算opt＝argmin{μ_i},i＝1,2,…,S，得到最优子集Γ_opt。