CN105206277A - 基于单比特压缩感知的语音压缩方法 - Google Patents

Abstract

本发明属于语言信号压缩领域，具体涉及一种基于单比特压缩感知的语音压缩方法。本发明的有益效果如下：1、语音信号的采样速率要求并不高，用于语音采样的AD并不是语音信号采样的瓶颈问题。其他需要利用压缩感知进行采样以降低AD采样速率的场景为高速采样场景，语音信号采样不属于高速采样，故可以不改变系统原有的采样器件与采样结构，只需要进行软件的修改，大大节约了系统改造的成本；2、目前流行的mp3技术，仅能压缩10余倍，而利用本发明的方法，可轻易地将语音压缩32倍以上，其音质也能达到近似于mp3标准音质。譬如电信领域的录音系统，需要存储海量的语音数据，如果利用本发明的压缩方法将会十分有利。

Description

基于单比特压缩感知的语音压缩方法

技术领域

本发明属于语言信号压缩领域，具体涉及一种基于单比特压缩感知的语音压缩方法。

背景技术

为了更方便的存储或者传输，必须对语音信号进行压缩。语音信号能够进行压缩的基本依据是语音信号中存在的冗余和人类的听觉感知系统的独特性。语音信号存在多种多样的冗余，可分别从时间域和频率域描述。从时间域分析，幅度的非均匀分布，即语音中的小幅度样本出现的概率高，信息主要集中在低功率上；采样数据间的相关，相邻的语音信号间有很强的相关性。从频率域分析，非均匀的长时功率谱密度，从相当长的时间内统计平均，语音信号的功率谱呈现强烈的非平坦性，这说明语音信号对给定的频段利用的不充分，存在固定的冗余度。语音压缩的第二个依据是人的听觉生理和心理特性对语音的感知，人的听觉对低频端比较敏感，而对高频端不太敏感，这主要是因为浊音的周期和共振峰在这里。此外，人耳对语音信号的相位变化不敏感，人耳的听觉不像人的视觉对感知的信号相位那样敏感，一定程度的相位失真对听觉来说是可以容忍的。人耳听不到或感知极不灵敏的语音信号都可以认为是冗余，可以利用这些特性进行语音压缩。

传统的信号采集以奈奎斯特采样定理作为理论依据，即信号的采样频率必须高于信号最高频率的两倍。近年来，E.J.Candes，J.Romberg，T.Tao和D.L.Donoho等学者提出了压缩感知理论(Compressedsensing,CS)，该理论指出，对于稀疏信号，采样时可以以远低于奈奎斯特采样频率进行信号采样，实现信号的压缩。由于语音信号在频域是可稀疏的，故可以通过压缩感知理论对其进行压缩和重构。

现有技术的技术方案

①对一段连续的语音信号采样N个点，得到长度为N的离散语音信号；

②对离散的语音信号进行DCT变换，将时域信号变换到频域，得到了长度为N的稀疏的语音信号；

③利用高斯分布生成M行N列的矩阵，2K≤M<<N；

④观测矩阵与频域信号相乘，由矩阵乘法易知，最后得到了M个语音信号，由于M<<N，故实现了语音信号的压缩。

现有技术的缺点

1、现有的基于压缩感知的语音压缩常见的为“非1bit压缩”方式，单比特的语音压缩还很少见；2、没有巧妙的应用基于压缩感知方式的心理声学模型，保存的数据存在冗余；3、重构精度难以保证；4、观测矩阵没有针对性。

缩略语和关键术语定义

稀疏度K：一帧频域语音信号中幅度不为0的个数；

AD：模数转换器；

DCT：DiscreteCosineTransform，离散余弦变换，简称DCT变换，由N.Ahmed等人在1974年提出的正交变换方法，它常被认为是对语音和图像信号进行变换的最佳方法；

IDCT变换：DCT变换的逆变换；

N：一帧离散语音信号的长度，因为语音是离散的，所以长度也即采样点的个数；

M：压缩后语音信号的长度，一般2K≤M<<N。

发明内容

本发明针对现有技术语音信号压缩时为“非1bit”且重构精度低而提出了一种基于单比特压缩感知的语音压缩方法。

为解决以上问题，本发明的技术方法是提供一种基于单比特压缩感知的语音压缩方法，包括以下步骤：

S1.将已采样的离散的N个语音信号从时域信号变换到频域信号，得到长度为N的稀疏的语音信号；

S2.将低于静音门限的信号幅度置0；

S3.查找幅度不为0的信号，记录其个数和位置，总个数即为稀疏度K，位置信息的集合记为I，并根据构造更优的观测矩阵；

S4.将S3生成的更优的观测矩阵与频域信号相乘，再对乘积结果取符号即可得到压缩语音信号，此时的信号长度由N变成了M，且每个信号占用1bit存储；

S5.将稀疏度K、位置信息的集合I与压缩的语音信号中的“1”登记在一起，组帧为数据流；

S6.将S5得到的数据流进行霍夫曼编码，进一步压缩。

作为优选，S1所述的将语音信号从时域信号变换到频域信号采用的方法为DCT变换、小波变换或者傅里叶变换。

作为优选，S3所述的记录语音信号位置信息的方法如下：以一帧长为N＝2048语音信号为例，如果要直接记录位置，那么记录一个位置则使用10bit，假设最后一个非0信号的位置为P，通过0,1标识则需要P个bit，记录位置则需要10K个bit，若10K>P，则采用方式0，反之采用方式1。

作为优选，S3所述的利用信号的位置信息，构造更优的观测矩阵方法如下：

S31.输入：稀疏基D，随机矩阵的行M，随机矩阵φ，位置集合I；

S32.初始化：感知矩阵F＝φD，迭代次数为n，i＝0，i是用来记录后续步骤中迭代次数的，大于n就结束迭代，否则继续；F^T是F的共轭转置矩阵，实数范围内即为F的转置矩阵；

S33.对F'每一列归一化，并把F'的值赋给F，F＝F'；

S34.计算矩阵，G＝F^TF；

S35.在位置集合I中任取其中的两个元素I_i和I_j，令G(i,j)＝0,要求集合I中的所有元素取遍；

S36.约束G中对角线外的元素大小:

如果G(i,j)的绝对值小于则保持不变，否则， $G(i,j)=sign\left(G2\right(i,j\left)\right)\sqrt{\left(\right(N-M)/(M(N-1)\left)\right)};$ sign表示取符号，如果G(i,j)大于0，则sign(G2(i,j))＝1，如果G(i,j)小于0，则sign(G2(i,j))＝-1，如果G(i,j)等于0，则sign(G2(i,j))＝0；

S37.SVD分解，G＝UQV；

S38.更新G，将G的秩降为M：G＝UQU^T；

S39.计算S＝Q^1/2U^T， $S^{\&prime;}=\sqrt{N/M}{\left(S{S}^T\right)}^{-\frac{1}{2}}S;$

S310.F'＝S'，更新表示矩阵D的伪逆；

S311.i＝i+1,i从0变为1，当i>n，则进入下一步，否则返回S33；

S312.迭代结束，输出优选之后的观测矩阵φ。

为解决以上问题，本发明还提出了一种基于单比特压缩感知的语音压缩方法的解压方法，包括以下步骤：

S7.对压缩的数据进行霍夫曼解码，得到0和1的数据流；

S8.根据压缩时组帧数据流的方法从数据流中取出稀疏度K和位置信息，构造观测矩阵；

S9.将信号的位置信息加入重构算法BITH，重构出频域的语音信号；

S10.将频域的语音信号从频域信号变换到时域信号。

作为优选，S9的具体步骤如下：

S91.输入一个M行1列的元素只有0,1的观测矩阵Y，即可得出更新步长α值，位置集合I和最大迭代次数t；

S92.初始化稀疏信号X，即X为一个N行1列的零矩阵；

S93.迭代更新X＝X+αφ'(Y-sign(φX))，φ'为φ的转置矩阵；

S94.X中位置集合I位置以外的元素都置零；

S95.若n≤t，则重复以上步骤，每迭代一次，n就加1，即n＝n+1；

S96.当n>t时，根据单位能量约束条件X＝X/||X||₂，输出X的最优估计值。

本发明的有益效果如下：

1、语音信号的采样速率要求并不高，用于语音采样的AD并不是语音信号采样的瓶颈问题。其他需要利用压缩感知进行采样以降低AD采样速率的场景为高速采样场景，语音信号采样不属于高速采样，故可以不改变系统原有的采样器件与采样结构，只需要进行软件的修改，大大节约了系统改造的成本；

2、目前流行的mp3技术，仅能压缩10余倍，而利用本发明的方法，可轻易地将语音压缩32倍以上，其音质也能达到近似于mp3标准音质。譬如电信领域的录音系统，需要存储海量的语音数据，如果利用本发明的压缩方法将会十分有利。

附图说明

图1为现有技术流程示意图；

图2为压缩流程示意图；

图3为重构流程示意图；

图4为记录语音信号位置的流程示意图；

图5为改进的BITH算法流程示意图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下参照附图并举实施例，对本发明做进一步详细说明。

压缩部分：如图2所示

①压缩部分的处理对象是已采样的离散语音信号，对离散的N个语音信号进行DCT变换，将时域信号变换到频域，得到长度为N的稀疏的语音信号；

②人的耳朵对声音的感知能力是有限的，有许多声音是人耳无法感知的。根据心理声学模型中的静音门槛曲线可知，低于静音门限的声音人耳是感知不到的，所以将低于静音门限的信号幅度置零，这样既不会影响音质，也实现了信号的压缩；

③查找幅度不为0的信号，记录其个数和位置，总个数即为稀疏度K，位置信息的集合I则是生成更优的观测矩阵的关键信息，也是高精度重构的重要保障。

④将生成的更优的观测矩阵与频域信号相乘，再对乘积结果取符号即可得到压缩语音信号，此时的信号不仅长度由N变成了M，而且每个信号仅需要占用1bit进行存储，实现了大幅度地压缩；

⑤将稀疏度K、位置集合I与压缩的语音信号中的“1”登记在一起，组帧为数据流。

⑥将⑤得到的数据流进行霍夫曼编码，进一步压缩。

重构部分(解压部分)：如图3所示

①对压缩的数据进行霍夫曼解码，得到0,1数据流；

②根据压缩部分组帧数据流的方法从数据流中取出K和位置信息，构造观测矩阵；

③利用改进的BITH算法重构出频域的语音信号；

④对频域的语音信号进行IDCT变换(DCT变换的逆变换)，得到时域的语音信号。

实施例1：记录语音信号的位置流程如图4所示：

本实施例以一帧长为N＝2048语音信号为例，如果要直接记录位置，那么记录一个位置则使用10bit，假设最后一个非0信号的位置为P，通过0,1标识则需要P个bit，记录位置则需要10K个bit，若10K>P，则采用方式0，反之采用方式1。例如，K＝3，非0位置是1,3,7，方式0需要7bit表示，即1010001，方式1需要3*10＝30bit表示，故采用方式0；同样K＝3，非0位置是1,3,100，方式0则需要100个bit表示，方式1还是只需要30个bit表示，故此时采用方式1。将方式号(0或者1)以及位置信息与压缩语音信号中的“1”一起组帧成数据流进行存储，如上述第一个例子，假设压缩语音信号中有4个“1”，组帧数据流为010100011111。

实施例2.利用信号的位置信息，构造更优的观测矩阵：

①输入：稀疏基D，随机矩阵的行M，随机矩阵φ，位置集合I。例如，D取为N行N列的单位矩阵(单位矩阵：是个方阵，从左上角到右下角的对角线(称为主对角线)上的元素均为1以外全都为0)， $D=\&lsqb;$

$\left.\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right],$

M＝2， $\mathrm{\&phi;}=\&lsqb;$

$\left.\begin{array}{cccc}0.5377& -2.2588& 0.3188& -0.4336\\ 1.8339& 0.8622& -1.3077& 0.3426\end{array}\right],$ I＝[1,2]

②初始化：感知矩阵F＝φD，迭代次数为n，i＝0。 $F=\&lsqb;$

$\left.\begin{array}{cccc}0.5377& -2.2588& 0.3188& -0.4336\\ 1.8339& 0.8622& -1.3077& 0.3426\end{array}\right],$ n＝100,i是用来记录后续步骤中迭代次数的，大于n就结束迭代，否则继续。

F^T是F的共轭转置矩阵，实数范围内即为F的转置矩阵， $F^T=\&lsqb;$

$\left.\begin{array}{cc}0.5377& 1.8339\\ -2.2588& 0.8622\\ 0.3188& -1.3077\\ -0.4336& 0.3426\end{array}\right],$ 所以 $F^{\&prime;}=\sqrt{N/M}{\left({\mathrm{FF}}^T\right)}^{-\frac{1}{2}}F=\sqrt{4/2}{\left({\mathrm{FF}}^T\right)}^{-\frac{1}{2}}F=\&lsqb;$

$\left.\begin{array}{cccc}0.4758& -1.3076& 0.0888& -0.2367\\ 1.1369& 0.3318& -0.7541& 0.1693\end{array}\right]$

③对F'每一列归一化，并把F'的值赋给F，F＝F'。比如F'的第一列，0.4758,1.1369，归一化过程：即为0.3860，0.9225，同理可得其他列。最后， $F=\&lsqb;$

$\left.\begin{array}{cccc}0.3860& -0.9693& 0.1170& -0.8133\\ 0.9225& 0.2460& -0.9931& 0.5818\end{array}\right]$

④计算矩阵 $G=F^{T}F=\&lsqb;$

$\begin{array}{cccc}1.0000& -0.1472& -0.8710& 0.2228\end{array}$

$\left.\begin{array}{cccc}-0.1472& 1.0001& -0.3577& 0.9315\\ -0.8710& -0.3577& 0.9999& -0.6729\\ 0.2228& 0.9315& -0.6729& 0.9999\end{array}\right]$

⑤在位置集合I中取遍两个元素I_i，I_j，令G(i,j)＝0,取的时候不重不漏。I＝[1,2]那么G(1,2),G(2,1)都置为0。 $G=\&lsqb;$

$\left.\begin{array}{cccc}1.0000& 0& -0.8710& 0.2228\\ 0& 1.0001& -0.3577& 0.9315\\ -0.8710& -0.3577& 0.9999& -0.6729\\ 0.2228& 0.9315& -0.6729& 0.9999\end{array}\right]$

⑥约束G中对角线外的元素大小:如果G(i,j)的绝对值小则保持不变，否则， $G(i,j)=sign\left(G2\right(i,j\left)\right)\sqrt{\left(\right(N-M)/(M(N-1)\left)\right)}.$ sign表示取符号，如果G(i,j)大于0，则sign(G2(i,j))＝1，如果G(i,j)小于0，则sign(G2(i,j))＝-1，如果G(i,j)等于0，则sign(G2(i,j))＝0。 $G=\&lsqb;$

$\left.\begin{array}{cccc}1.0000& 0& -0.5774& 0.2228\\ 0& 1.0001& -0.3577& 0.5774\\ -0.5774& -0.3577& 0.9999& -0.5774\\ 0.2228& 0.5774& -0.5774& 0.9999\end{array}\right]$

⑦SVD分解， $G=UQV,U=\&lsqb;$

$\left.\begin{array}{cccc}-0.3835& 0.7073& -0.3668& -0.4671\\ -0.4425& -0.5925& -0.6732& -0.0052\\ 0.5840& -0.2593& -0.1498& -0.7545\\ -0.5622& -0.2855& 0.6244& -0.4610\end{array}\right],$

$Q=\&lsqb;$

$\left.\begin{array}{cccc}2.2059& 0& 0& 0\\ 0& 1.1217& 0& 0\\ 0& 0& 0.3850& 0\\ 0& 0& 0& 0.2873\end{array}\right],$

$\begin{array}{c}V=\&lsqb;\\ \left.\begin{array}{cccc}-0.3835& 0.7073& -0.3668& -0.4671\\ -0.4425& -0.5925& -0.6732& -0.0052\\ 0.5840& -0.2593& -0.1498& -0.7545\\ -0.5622& -0.2855& 0.6244& -0.4610\end{array}\right]\end{array}$

⑧更新G，将G的秩降为M： $G={\mathrm{UQU}}^T=\&lsqb;$

$\left.\begin{array}{cccc}1.0000& -0.0000& -0.5774& 0.2228\\ -0.0000& 1.0001& -0.3577& 0.5774\\ -0.5774& -0.3577& 0.9999& -0.5774\\ 0.2228& 0.5574& -0.5774& 0.9999\end{array}\right]$

⑨计算 $S=Q^{1/2}{U}^T=\&lsqb;$

$\left.\begin{array}{cccc}-0.5696& -0.6571& 0.8674& -0.8351\\ 0.7491& -0.6275& -0.2746& -0.3024\\ -0.2276& -0.4177& -0.0929& 0.3874\\ -0.2503& -0.0028& -0.4044& -0.2471\end{array}\right],$

$\begin{array}{c}{S}^{\&prime;}=\sqrt{N/M}{\left(S{S}^T\right)}^{-\frac{1}{2}}S=\&lsqb;\\ \begin{array}{cccc}-0.5423& -0.6257& 0.8259& -0.7951\\ 1.0002& -0.8379& -0.3667& -0.4037\\ -0.5187& -0.9520& -0.2118& 0.8830\end{array}\end{array}$

$\left.\begin{array}{cccc}-0.6605& -0.0074& -1.0670& -0.6520\end{array}\right]$

⑩F'＝S'，更新

表示矩阵D的伪逆

i＝i+1,i从0变为1，小于n＝100，返回第③步，直到i>n。

迭代结束，输出优选之后的观测矩阵φ。

实施例3.巧妙地将信号的位置信息加入重构算法BITH，保证了重构信号的高精度，改进的BITH算法流程如图5所示：

1)输入一个4行1列(M行1列)的元素只有0,1的矩阵Y

数据流 $Y=\left[\begin{array}{c}1\\ 1\\ 0\\ 1\end{array}\right];$

观测矩阵 $\mathrm{\&Phi;}=\left[\begin{array}{cccccccc}0.6892& 0.2290& 0.5383& 0.1067& 0.8173& 0.2599& 0.1818& 0.8693\\ 0.7482& 0.9133& 0.9961& 0.9619& 0.8687& 0.8001& 0.2638& 0.5797\\ 0.4505& 01524& 0.0782& 0.0046& 0.0844& 0.4314& 0.1455& 0.5499\\ 0.0838& 0.8258& 0.4427& 0.7749& 0.3998& 0.9106& 0.1361& 0.1450\end{array}\right];$

更新步长α＝1；

位置集合I＝[13]；

最大迭代次数t＝300。

2)初始化稀疏信号X，即X为一个8行1列(N行1列)的零矩阵

$X=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right];$

更新迭代次数n＝1。

3)迭代更新X＝X+αφ'(Y-sign(φX))；

φ'为φ的转置，即 ${\mathrm{\&phi;}}^{\&prime;}=\left[\begin{array}{cccc}0.6892& 0.7482& 0.4505& 0.0838\\ 0.2290& 0.9133& 0.1524& 0.8258\\ 0.5383& 0.9961& 0.1782& 0.4427\\ 0.1067& 0.9619& 0.0046& 0.7749\\ 0.8173& 0.8687& 0.0844& 0.3998\\ 0.2599& 0.8001& 0.4314& 0.9106\\ 0.1818& 0.2638& 0.1455& 0.1361\\ 0.8693& 0.5797& 0.5499& 0.1450\end{array}\right];$

第一次迭代更新的 $X=\left[\begin{array}{c}1.5212\\ 1.9681\\ 1.9772\\ 1.8435\\ 2.0858\\ 1.9706\\ 0.5817\\ 1.5940\end{array}\right].$

4)X中位置集合I＝[137]位置以外的元素都置零

$X=\left[\begin{array}{c}1.5212\\ 0\\ 1.9772\\ 0\\ 0\\ 0\\ 0.5817\\ 0\end{array}\right].$

5)若n≤t，则重复以上步骤，每迭代一次，n就加1(n＝n+1)。

6)当n>t时，根据单位能量约束条件X＝X/||X||₂，输出X的最优估计值。比如对于1.5212，单位能量约束即为：

如果没有获取先验信息，随机生成观测矩阵，也能将压缩率提高，但是重构音质会大打折扣。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的实施方法，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (6)

1.一种基于单比特压缩感知的语音压缩方法，其特征在于，包括以下步骤：

S1.将已采样的离散的N个语音信号从时域信号变换到频域信号，得到长度为N的稀疏的语音信号；

S2.将低于静音门限的信号幅度置0；

S3.查找幅度不为0的信号，记录其个数和位置，总个数即为稀疏度K，位置信息的集合记为I，并根据构造更优的观测矩阵；

S4.将S3生成的更优的观测矩阵与频域信号相乘，再对乘积结果取符号即可得到压缩语音信号，此时的信号长度由N变成了M，且每个信号占用1bit存储；

S5.将稀疏度K、位置信息的集合I与压缩的语音信号中的“1”登记在一起，组帧为数据流；

S6.将S5得到的数据流进行霍夫曼编码，进一步压缩。

2.根据权利要求1所述的一种基于单比特压缩感知的语音压缩方法，其特征在于，S1所述的将语音信号从时域信号变换到频域信号采用的方法为DCT变换、小波变换或者傅里叶变换。

3.根据权利要求1或2所述的一种基于单比特压缩感知的语音压缩方法，其特征在于，S3所述的记录语音信号位置信息的方法如下：以一帧长为N＝2048语音信号为例，如果要直接记录位置，那么记录一个位置则使用10bit，假设最后一个非0信号的位置为P，通过0,1标识则需要P个bit，记录位置则需要10K个bit，若10K>P，则采用方式0，反之采用方式1。

4.根据权利要求3所述的一种基于单比特压缩感知的语音压缩方法，其特征在于，S3所述的利用信号的位置信息，构造更优的观测矩阵方法如下：

S31.输入：稀疏基D，随机矩阵的行M，随机矩阵Ф，位置集合I；

S32.初始化：感知矩阵F＝φD，迭代次数为n，i＝0，i是用来记录后续步骤中迭代次数的，大于n就结束迭代，否则继续；F^T是F的共轭转置矩阵，实数范围内即为F的转置矩阵；

S33.对F'每一列归一化，并把F'的值赋给F，F＝F'；

S34.计算矩阵，G＝F^TF；

S35.在位置集合I中任取其中的两个元素I_i和I_j，令G(i,j)＝0,要求集合I中的所有元素取遍；

S36.约束G中对角线外的元素大小:

如果G(i,j)的绝对值小于则保持不变，否则，

$G(i,j)=sign\left(G2\right(i,j\left)\right)\sqrt{\left(\right(N-M\left)\right)/\left(M\right(N-1\left)\right)};$ sign表示取符号，如果G(i,j)大于0，则sign(G2(i,j))＝1，如果G(i,j)小于0，则sign(G2(i,j))＝-1，如果G(i,j)等于0，则sign(G2(i,j))＝0；

S37.SVD分解，G＝UQV；

S38.更新G，将G的秩降为M：G＝UQU^T；

S39.计算S＝Q^1/2U^T， $S^{\&prime;}=\sqrt{N/M}{\left(S{S}^T\right)}^{-\frac{1}{2}}S;$

S310.F'＝S'，更新表示矩阵D的伪逆；

S311.i＝i+1,i从0变为1，当i>n，则进入下一步，否则返回S33；

S312.迭代结束，输出优选之后的观测矩阵φ。

5.一种根据权利要求1所述的基于单比特语音压缩方法的解压方法，其特征在于，包括以下步骤：

S7.对压缩的数据进行霍夫曼解码，得到0和1的数据流；

S8.根据压缩时组帧数据流的方法从数据流中取出稀疏度K和位置信息，构造观测矩阵；

S9.将信号的位置信息加入重构算法BITH，重构出频域的语音信号；

S10.将频域的语音信号从频域信号变换到时域信号。

6.根据权利要求5所述的一种解压方法，其特征在于，S9的具体步骤如下：

S91.输入一个M行1列的元素只有0,1的观测矩阵Y，即可得出更新步长α值，位置集合I和最大迭代次数t；

S92.初始化稀疏信号X，即X为一个N行1列的零矩阵；

S93.迭代更新X＝X+αφ'(Y-sign(φX))，φ'为φ的转置矩阵；

S94.X中位置集合I位置以外的元素都置零；

S95.若n<t，则重复以上步骤，每迭代一次，n就加1，即n＝n+1；

S96.当n>t时，根据单位能量约束条件X＝X/||X||₂，输出X的最优估计值。