-
与小波变换相比,小波包变换为信号提供了更精细的分析方法。对仍含有丰富信息特征的高频信号也进行分解,并能根据信号的具体特征自适应地选择频带来提高高频信号的分辨率。以信号的3层小波包分解为例来分析其分解过程,如图 1所示。
图 1中,用L和H分别表示小波分解后的低频和高频部分,则原始信号S0经第1层小波包分解后得到低频信号S1, L和高频信号S1, H;对S1, L再进行一层小波包分解后得到下一尺度的低频信号S2, LL和高频信号S2, HL,对S1, H进行一层小波包分解后得到下一尺度的低频信号S2, LH和高频信号S2, HH;对低频部分S2, LL进行第3层分解后的低频S3, LL和高频S3, HLL部分……。依次类推,可以得到N层小波包分解的低频信号和高频信号。假设正交尺度函数、正交小波函数分别为ϕ(t), φ(t),则小波包变换的双尺度方程[10]可表示为:
$ \left\{ {\begin{array}{*{20}{l}} {\phi (t) = \sqrt 2 \sum\limits_k h (k)\phi (2t - k)}\\ {\varphi (t) = \sqrt 2 \sum\limits_k g (k)\varphi (2t - k)} \end{array}} \right. $
(1) 式中,t为时间参量;k为平移参量;h(k), g(k)分别为多尺度分析中低、高通滤波器系数,且g(k)=(-1)k×h(1-k),即两系数具有正交性。
为进一步推广,由(1)式得到下列递推关系:
$ \left\{ {\begin{array}{*{20}{l}} {{\omega _{2n}}(t) = \sqrt 2 \sum\limits_k h (k){\omega _n}(2t - k)}\\ {{\omega _{2n + 1}}(t) = \sqrt 2 \sum\limits_k g (k){\omega _n}(2t - k)} \end{array}} \right. $
(2) 式中,当n=0时,递归函数ω0(t)退化为尺度函数ϕ(t),ω1(t)退化为小波函数φ(t),即ω0(t)=ϕ(t),ω1(t)=φ(t)。所以,(2)式所表示的函数集合{ωn(t)}(n=0, 1, …, n)可以看成是基于ω0(t)=ϕ(t)所确定的小波包。或者说,小波包就是一个包括ϕ(t), φ(t)在内的,满足双尺度递推关系的函数集合。其中,小波包系数的递推关系为:
$ \left\{ {\begin{array}{*{20}{l}} {{d_{j + 1, 2n}} = \sum\limits_k h (k - 2t){d_{j, n}}(k)}\\ {{d_{j + 1, 2n + 1}} = \sum\limits_k g (k - 2t){d_{j, n}}(k)} \end{array}} \right. $
(3) 式中,dj, n(k)为小波包经分解后第(j, n)节点所对应的第k个系数,节点(j, n)表示第j层分解的第n个频带。
由上述分析可知,小波包分解过程中对低频部分和高频部分进行了同步分解,对高频部分的分解提供了更为精细的信号分析方法[11]。
-
对同一信号来说,其小波包变换后在频域的稀疏度越高,其压缩率越高,且重构精度越高。而频域的稀疏度与小波包变换过程中所选取的小波基类型关系密切。目前,选择什么样的小波基能使信号在频域具有较高的稀疏度的方法研究还很少,还没有统一的选取标准。本文中采用实验的方法来选择适合对光纤振动信号进行稀疏分解的小波基。首先,从一套光纤周界报警系统中采样(采样频率20kHz)得到包括4种常见扰动(分别为正常信号、小扰动信号、大扰动信号、蓄意入侵信号,数据长度为2018×256个采样点)类型的光纤扰动信号作为测试样本。从测试样本中任意截取一段长度为256个采样点的振动信号,分别采用几种常用的小波基进行4尺度小波包分解,分别计算4尺度下的数学期望Ep(p=1, 2, 3, 4)作为阈值Tp对小波包系数进行置零处理,并分别计算其稀疏度kp(即系数中为零系数占所有系数的百分比),实验结果如表 1所示。比较表 1中数据,本文中选择sym7小波基对光纤振动信号进行小波包变换。
Table 1. Sparse degree of different wavelets
wavelet basis DB4 DB6 sym3 sym5 sym7 sparse degree/% 64.65 64.73 63.64 65.72 67.21 -
利用小波包对信号进行分解的层数与小波包系数的稀疏度密切相关。在实际应用中,应根据信号的特征选择合适的分解层数,以便获得最大信号稀疏度。本文中用实验的方法来分析小波包的分解层数与小波包系数稀疏度之间的关系。在第1.2节中的测试样本任意截取50段长度为256个采样点的光纤振动信号,分别采用sym7小波基对其进行2~6层小波包分解,分别计算各尺度下高频系数的数据期望作为阈值,对该尺度下的小波包系数进行置零处理,然后计算各尺度下小波包系数的平均稀疏度, 如表 2所示。
Table 2. Sparse degree of different decomposition layers/%
decomposition layers 2 3 4 5 6 sparse degree 57.17 62.81 66.32 64.47 63.26 由表 2中数据可以看出,随着分解层数的增加,小波包系数的稀疏度也增加;但到达一定层数时,稀疏度开始下降。所以,在实际应用中,必须对小波包变换的分解层数进行合理的选择,才能使信号在频域的稀疏度达到最大。
-
设有限长、非稀疏信号x∈RN(RN表示由N个元素所表示的实数集合)经某种变换后可由一组正交基ψ={ψi}i=1N线性表示,即:
$ \mathit{\boldsymbol{x}} = \sum\limits_{i = 1}^N {{s_i}} {\psi _i} $
(4) 式中,s=〈x, ψi〉=ψiTx是该信号在ψ域的稀疏系数,T表示向量的转置,包含有k个非零项,且k < < N;x为在正交基矩阵ψ上的k稀疏表示。
构建一个和正交基矩阵ψ不相关的观测矩阵ΦM×N={ϕ1T, ϕ1T, …, ϕMT},M为ΦM×N矩阵的行向量数,且M < < N,将信号x由高维空间投射到低维空间[12]:
$ y = \mathit{\boldsymbol{ \boldsymbol{\varPhi} x}} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{\boldsymbol{\psi }}_s} $
(5) 式中,y∈RM×1称为观测向量;ΦΨ称为传感矩阵,s的含义同(4)式。
根据压缩感知理论,信号的重构是一个欠采样条件下病态的求逆问题,即信号x是下式l0最小化问题的解,即:
$ \mathit{\boldsymbol{\widehat x}} = \arg \min ||\mathit{\boldsymbol{x}}|{|_0}\quad {\rm{ s}}{\rm{.t}}{\rm{. }}\mathit{\boldsymbol{y = \boldsymbol{\varPhi} x}} $
(6) 为了更容易求解,将(6)式的NP-hard问题转化为求解最小化l1的问题[13-14]。目前常用的方法有正交匹配追踪(orthogonal matching pursuit,OMP)算法[15]、基追踪(basis pursuit,BP)算法[16]等。
-
为了充分体现小波包系数中大系数的分布情况,需重新定义数学期望的计算公式[17]。设数据序列x={xi}(i=1, 2, …, n),定义其数学期望为:
$ \left\{ {\begin{array}{*{20}{l}} {E = \sum\limits_i^n {\left| {{x_i}} \right|} {p_i}}\\ {p = 1/n} \end{array}} \right. $
(7) 对光纤振动信号进行多尺度小波包分解后,利用(7)式计算小波包系数高频部分的数学期望作为阈值,对小波包系数进行置零处理。例如,任意截取第1.2节中所述测试样本中的一段长度为256个采样点的光纤振动信号进行2尺度小波包分解,2尺度下小波包系数各系数块的数学期望值为{613.54, 19.44, 3.72, 4.36},其中613.54为低频部分的数学期望值,高频各部分的数学期望中最大值为19.44,取该值作为2尺度下小波包系数的阈值T,对该尺度下的小波包系数进行置零处理,即:
$ x = \left\{ {\begin{array}{*{20}{l}} {x, (|x| \ge T)}\\ {0, (|x| < T)} \end{array}} \right. $
(8) 关于阈值选取对信号重构精度及压缩率的影响,在后面的实验分析中还将进一步讨论。
-
光纤振动信号经小波包分解后,小波包系数的稀疏度可用数学期望表示,当小波包系数的数学期望越大时,说明大系数越多,其稀疏度越低[10, 18]。根据第1.3节中的分析,采用小波包对光纤振动信号进行多层分解后,各尺度下的小波包系数的稀疏度是不同的,为了使信号在频域具有较高的稀疏度,采用相同的分解层数显然是不合理的。因此,本文中基于小波包系数的数学期望提出自适应小波包变换。方法如下:首先对光纤振动信号进行3层~5层小波包分解,然后计算各尺度下小波包系数高频部分的数学期望值,按照(5)式对各尺度下的小波包系数进行阈值处理。比较各尺度下的稀疏度,选择稀疏度最大的尺度作为信号的最终分解层数,从而完成自适应小波变换。
-
为了进一步提高信号的压缩率和重构精度,在对小波包系数进行阈值处理后,重新计算该尺度下各小波包系数块高频部分的数学期望E和信息熵H[19],并根据E和H对小波包系数块进行分类。
(1) 低频系数块。是指小波包系数块中的低频部分,包含着信号的最主要信息,包含了大量的大系数,属于不稀疏信号,如果仍采用压缩感知方法进行处理,不但会增加较大的计算量,也会造成信号重构精度的降低。所以这部分系数块采用线性直接传输方式,一方面保证了信号的重构精度,另一方面也降低了算法的运算时间。
(2) 无价值系数块。当系数块的数学期望E=0时,说明该系数块全为0,这类系数块在信号的数据压缩、传输和重构过程中都没有利用价值,所以没有考虑的必要,称为无价值系数块,这类系数块不参与数据的传输与存储。
(3) 特殊系数块。如果系数块的信息熵为零或数学期望小于某一特定值C时,归类为特殊系数块。
根据信息熵的定义,当系数块的信息熵为零时,说明该系数块要么全为0,要么只有一个不为零的大系数;又因为其数学期望不为0,则易知该系数块中只有一个不为零的大系数。
当系数块的数学期望小于某一特定值C时,C值大小应根据实际需要进行确定。对有N个系数的系数块,设阈值为T,其中包含有k个大系数;为了确定C值,假定该系数块中的k个大系数均为阈值T,此时该系数块的数学期望为kT/N,那么该系数块的C=kT/N;显然,当某系数块的数学期望小于C值时,则其大系数的个数必然小于k。
由以上分析可知,特殊系数块就是指只包含有少数几个大系数的小波包系数块。在传输时,特殊系数块只需要传输大系数值及其相应的位置,即采用编码方式进行网络传输。
(4) 压缩感知系数块。除无价值系数块和特殊系数块外的其系数块都归类为压缩感知系数块。这类系数块的数学期望和信息熵均不为零,说明其中包含的大系数较多,采用压缩感知方法进行处理比较合适。
-
信号的重构工作包括小波包系数重构和小波包重构。小波包系数重构根据系数块的类型选择相应的重构方法:(1)低频系数块直接恢复,并摆放在小波包系数的低频部分;(2)特殊系数块按系数值及位置编码恢复小波包系数及其所在位置,其余位置用0填充;(3)压缩感知系数块采用压缩感知重构方法(例如OMP)恢复数据,并摆放在小波包系数块的对应位置;(4)小波包系数块的剩余位置用0填充。
完成小波包系数恢复(重构)后,再对小波包系数进行小波包逆变换,完成信号重构。
-
为了比较光纤振动信号的压缩感知和重构性能,引入信号压缩比(compression ratio,CR)RCR评价信号的压缩效果:
$ {R_{{\rm{CR}}}} = \frac{{{n_1} - {n_2}}}{{{n_1}}} $
(9) 式中,n1, n2分别为信号压缩前、后的数据量。
引入PSNR(RPSNR)作为信号重构精度的评价指标:
$ \left\{ {\begin{array}{*{20}{l}} {{R_{{\rm{PSNR}}}} = 10\lg \left( {\frac{{{y_{{{\max }^2}}}}}{{{\delta ^2}}}} \right)}\\ {{\delta ^2} = \frac{1}{m}\sum\limits_{i = 1}^m {{{\left( {{{\hat y}_i} - {y_i}} \right)}^2}} } \end{array}} \right. $
(10) 式中,$y,\hat y$分别为信号重构前、后的采样数据;m为信号的采样点数。
通过实验得到,对周界报警信号来说,采用PSNR作为信号重构评价指标时,当重构后信号的RPSNR>40dB时,信号的特征信息能够得到完整的保留。
-
从第1.2节中所述的测试样本中任意截取5段长度为256个采样点的光纤振动信号,分别对其进行4层小波分解,计算4尺度下小波系数高频部分的数学期望E,分别选用1.2E,E,0.8E作为阈值T,对小波包系数高频部分进行置零处理。然后采用本文中的小波系数分类和信号重构方法得到重构后的光纤振动信号。分别计算3种阈值T下的信号压缩率/重构精度, 结果如表 3所示。
Table 3. Effect of threshold size on signal compression/reconstruction accuracy
different signal segment the zero threshold of wavelet packet coefficient 1.2E E 0.8E signal segment 1/(%/dB) 73.8/40.4 68.5/44.7 65.2/44.6 signal segment 2/(%/dB) 73.4/41.2 70.1/43.9 66.1/43.7 signal segment 3/(%/dB) 67.5/42.7 63.8/45.2 59.9/44.3 signal segment 4/(%/dB) 72.2/40.5 69.2/42.0 63.6/42.2 signal segment 5/(%/dB) 71.6/39.2 67.3/43.3 62.4/41.9 由表 3中数据可以看出,阈值较小时,小波包系数稀疏度较低,压缩率也低;阈值较大时,小波包系数稀疏度提高,压缩率提高。特别是,当阈值为1.2E时,信号的重构精度较其它两种方法低,这是因为阈值过大,基于该阈值处理小波包系数会使稀疏过度,丢失原始信号中一些有用的特征信息,造成信号重构精度下降;而当阈值为0.8E时,信号重构精度较阈值为E时稍低,这是因为前者虽然保留了更多信号的特征信息,但由于稀疏度下降,所以重构精度也有所下降。综上所述,选择小波包系数的数学期望E为阈值时,其压缩率和精度都比较高,所以该阈值大小选取最为合理。
-
仍采用第1.2节中的测试样本信号,任意截取256个采样点,分别采用小波变换和小波包变换进行固定4尺度分解,并分别选取4尺度下小波系数、小波包系数据的数学期望作为阈值对小波系数、小波包系数进行置零处理。再将该段信号采用本文中所述算法进行处理。在保证3种算法具有相同压缩率的前提下,对信号进行重构的结果如图 3所示。
Figure 3. Signal reconstruction effect of different algorithms a—original signal b—the algorithm in this paper c—wavelet transform d—wavelet packet transform
从图 3可以看出,3种算法在相同压缩率前提下,小波变换方法的信号重构效果最差;而小波包算法和本文算法的重构效果相差不大,但将波形放大后进行比较,还是可以看出本文中所述算法的重构效果优于小波包算法。
-
从第1.2节中的测试样本中分别截取光纤传感器常见的4种输出信号类型:正常信号、小扰动信号、大扰动信号和入侵信号各50段,每段信号长度均为256个采样点。分别采用本文中算法对信号处理,并计算4种信号的平均压缩率, 如表 4所示。
Table 4. Compressibility of fiber vibration signal
signal types normal signal small perturbation large disturbance intrusion signal compressibility 76.67% 76.65% 75.82% 75.37% 从表 4中的数据可以看出,4种光纤振动信号的压缩率基本一致,可见本文中算法是适合对光纤周界报警系统的各种输出信号进行处理的。
-
从第1.2节中的测试样本中任意截取1500段长度为256个采样点的光纤振动信号,首先采用sym7小波基分别对1500段信号进行4尺度小波、4尺度小波包稀疏变换,并计算4尺度下的小波系数、小波包系数的数学期望作为阈值,并对小波系数、小波包系数进行置零处理,并采用OMP方法完成信号重构;然后采用sym7小波基利用本文中算法对1500段信号进行处理。分别计算3种算法的平均压缩率、信号的平均重构精度以及各算法处理1500段信号总消耗时间,计算结果如表 5所示。
Table 5. Comprehensive performance of various algorithms
algorithm types compressibility/% reconstruction accuracy/dB elapsed time/ms wavelet 53.36 38.28 1342 wavelet packet 55.46 40.14 1514 this paper 56.25 41.46 891 比较表 5中的数据明显可以看出,小波包算法的压缩率和重构精度优于小波算法,这是因为小波包算法对信号的高、低频部分进行同频分解,信号的稀疏度提高,且保留了更多信号的特征信息。但小波包算法由于增加了信号的处理过程,因而在计算耗时上比小波算法略有增加;本文中算法采用数学期望作为阈值对小波包系数进行置零处理,并自适应地选取了最佳分解尺度,进一步提高了信号的稀疏度,因而算法的压缩率和重构精度比小波包算法更好。同时,本文中算法利用数学期望及信息熵对小波包系数进行分类处理,优化了信号处理方法,提高了信号处理速度,因而运算速度得到提高。
由于时域下的光纤报警信号表现复杂,各分块信号的信息特征差异较大,采用固定的稀疏字典对各分块信号进行稀疏表示仍不可能使信号的稀疏度达到最优。因而可以考虑根据信号的具体特征构建过完备字典,使信号在变换域的稀疏度进一步提高,从而可以进一步提高信号的压缩率和重构精度。这是下一步工作的研究重点。
光纤周界报警信号自适应压缩感知
Adaptive compression sensing of optical fiber perimeter alarm signal
-
摘要: 在光纤周界报警系统中,对光纤振动信号进行分析与辨识时,针对高频、大规模信号在采样、存储、传输与信号处理过程中存在网络宽带、存储容量、计算速度等一系列限制问题,提出基于小波包的光纤周界报警信号自适应压缩感知方法。对光纤振动信号进行多尺度小波包分解,通过计算各尺度下小波包系数高频部分的数学期望作为阈值,对小波包系数进行置零处理,自适应地选择小波包分解尺度,使信号在频域得到更高的稀疏度;根据小波包系数块的数学期望和信息熵对小波包系数块进行分类,并针对不同系数块的类型设计对应的处理方法,提高信号的传输与处理速度。结果表明,该方法能够有效减少光纤振动信号的观测数据,并在相同采样率前提下,能够提高信号的重构精度和重构速度。Abstract: In an optical fiber perimeter alarm system, when analyzing and identifying the fiber vibration signal, there are a series of limitations such as network broadband, storage capacity and computing speed in the process of sampling, storage, transmission and signal processing of high-frequency large-scale signal. In order to solve this problem, an adaptive compression sensing method of optical fiber perimeter alarm signal based on wavelet packet was proposed. Firstly, multi-scale wavelet packet decomposition was used to decompose the optical fiber vibration signal. By calculating the mathematical expectation of the high frequency part of the wavelet packet coefficients at different scales as the threshold value, wavelet packet coefficients were set to zero. The wavelet packet decomposition scale was adaptively selected so that the signal sparse in frequency domain. Then, wavelet packet coefficients were classified according to the mathematical expectation and information entropy of the wavelet packet coefficients. According to different types of coefficient blocks, the corresponding processing methods were designed to improve the speed of signal transmission and processing. The results show that this method can effectively reduce the observation data of optical fiber vibration signal. At the same sampling rate, it can improve the accuracy and speed of signal reconstruction.
-
Table 1. Sparse degree of different wavelets
wavelet basis DB4 DB6 sym3 sym5 sym7 sparse degree/% 64.65 64.73 63.64 65.72 67.21 Table 2. Sparse degree of different decomposition layers/%
decomposition layers 2 3 4 5 6 sparse degree 57.17 62.81 66.32 64.47 63.26 Table 3. Effect of threshold size on signal compression/reconstruction accuracy
different signal segment the zero threshold of wavelet packet coefficient 1.2E E 0.8E signal segment 1/(%/dB) 73.8/40.4 68.5/44.7 65.2/44.6 signal segment 2/(%/dB) 73.4/41.2 70.1/43.9 66.1/43.7 signal segment 3/(%/dB) 67.5/42.7 63.8/45.2 59.9/44.3 signal segment 4/(%/dB) 72.2/40.5 69.2/42.0 63.6/42.2 signal segment 5/(%/dB) 71.6/39.2 67.3/43.3 62.4/41.9 Table 4. Compressibility of fiber vibration signal
signal types normal signal small perturbation large disturbance intrusion signal compressibility 76.67% 76.65% 75.82% 75.37% Table 5. Comprehensive performance of various algorithms
algorithm types compressibility/% reconstruction accuracy/dB elapsed time/ms wavelet 53.36 38.28 1342 wavelet packet 55.46 40.14 1514 this paper 56.25 41.46 891 -
[1] WANG W, ZHANG T Y, CHEN H B, et al. High-speed demodulation system of fiber Bragg grating temperature sensor based on FPGA[J]. Laser Technology, 2018, 42(6): 849-853(in Chinese). [2] CHEN Z, SAMARABANDU J, RODRIGO R. Recent advances in simultaneous localization and map-building using computer vision[J]. Advanced Robotics, 2007, 21(3/4): 233-265. [3] DARDARI D, CLOSAS P, DJURIC P M. Indoor tracking theory, methods, and technologies[J]. IEEE Transactions on Vehicular Technology, 2015, 64(4): 1263-1278. doi: 10.1109/TVT.2015.2403868 [4] YANG X F, CHENG Y Y, WANG G. Super-resolution algorithm for remote sensing images based on compressive sensing in wavelet domain[J]. Journal of Computer Applications, 2017, 37(5): 1430-1433(in Chinese). [5] WANG Q, ZHANG P L, WANG H G, et al. Adaptive blocked compressed sensing algorithm for the machinery vibration signal[J]. Ch-inese Journal of Scientific Instrument, 2017, 38(2):313-319(in Chinese). [6] YAN G, SUN H. Identification of impact force for composite structure using Bayesian compressing sensing[J]. Journal of Vibration Engineering, 2018, 31(3): 484-489(in Chinese). [7] WU J N, XU H D, WANG J. A new joint reconstruction algorithm of compressed sensing for multichannel EEG signals based on over complete dictionary approach[J]. Journal of Electronics & Information Technology, 2016, 38(7): 1666-1673(in Chinese). [8] GAO Y, ZANG M X, GUO F Y. Research on speech compression based on wavelet transform and compressed sensing Application[J]. Research of Computers, 2017, 34(12): 3672-3680(in Chinese). [9] LIU Y, TANG W, LIU B Q. Data sparse analysis and improved reconstruction algorithm of power quality disturbance based on compressed sensing[J]. Transactions of China Electro-technical Society, 2018, 15(8): 3461-3470(in Chinese). [10] YUAN L F, SUN Y Sh, HE Y G, et al. Fault feature extraction method for analog circuit based on preferred wavelet packet[J]. Transactions of China Electro-technical Society, 2018, 33(1): 158-168(in Chinese). [11] XIONG N, CAO Y L. Detection and analysis of power system harmonics based on wavelet packet decomposition[J]. Journal of Shanghai University of Electric Power, 2017, 33(1): 60-64(in Chinese). [12] CHOI H, YANG K W, KIM E. Simultaneous global localization and mapping[J]. IEEE/ASME Transactions on Mechatronics, 2014, 19(4): 1160-1170. doi: 10.1109/TMECH.2013.2274822 [13] MAECHLER P, STUDER C, BELLASO D, et al. VLSI design of approximate message passing for signal restoration and compressive sensing[J]. IEEE Journal Selected Topic on Signal Processing, 2012, 2(3): 570-590. [14] GRIBONVAL R, CEVHER V, DAVIES M E. Compressible distributions for high-dimensional statistics[J]. IEEE Transactions on Information Theory, 2012, 58(8): 5016-5034. doi: 10.1109/TIT.2012.2197174 [15] LU Y M, DO M N. Sampling signals from a union of subspaces[J]. IEEE Signal Processing Magazine, Special Issue on Compressive Sampling, 2008, 25(2): 41-47. doi: 10.1109/MSP.2007.914999 [16] DU Zh H, CHEN X F, ZHANG H, et al. Feature identification with compressive measurements for machine fault diagnosis[J]. IEEE Transactions on Instrumentation and Measurement, 2016, 65(5): 977-987. doi: 10.1109/TIM.2016.2521223 [17] HU K Y, WANG J, WANG Y. Image encryption based on block compression sensing and the improved magic square transformation[J]. Laser Technology, 2019, 43(4): 532-538 (in Chinese). [18] ZHA Sh M, ZHU Zh J, CHI H. Research of key issues of photonic-assisted compressive sensing technology[J]. Laser Technology, 2016, 40(4): 565-570(in Chinese). [19] GUO H, YANG K M, ZHANG W W, et al. Classification of a hyperspectral image based on wavelet packet entropy feature vector angle[J]. Journal of Image and Graphics, 2017, 22(2): 205-211(in Chinese). [20] CANDES E J, STROHMER T, VORONINSKI V. Phaselift: Exact and stable signal recovery from magnitude measurement via convex programming[J]. Communications on Pure and Applied Mathema-tics, 2013, 66(8): 1241-1274. doi: 10.1002/cpa.21432