季策, 孙梦雪, 张君
东北大学 计算机科学与工程学院, 辽宁 沈阳 110169
收稿日期：2017-04-07
基金项目：国家自然科学基金资助项目(61370152, 61671141, 61501038)；沈阳市科技计划项目(F16-205-1-01)。
作者简介：季策(1969-)，女，辽宁沈阳人, 东北大学副教授。

摘要：为改善拟合效果, 针对经验模态分解(empirical mode decomposition, EMD)算法存在的端点效应, 提出一种改进的EMD算法——端点极值延拓方法.利用改进的EMD算法对观测信号进行分解, 将分解分量连同之前的观测信号构成新的观测信号, 从而将欠定情况转化为超定情况, 最后利用独立成分分析(independent component analysis, ICA)算法得到源信号的估计.通过仿真实验对比, 证明了本文算法的有效性.
关键词：欠定盲源分离经验模态分解端点效应极值延拓独立成分分析
Underdetermined Blind Separation Based on Improved EMD Algorithm
JI Ce, SUN Meng-xue, ZHANG Jun
School of Computer Science & Engineering, Northeastern University, Shenyang 110169, China
Corresponding author: SUN Meng-xue, E-mail: smx_1564@163.com
Abstract: In order to improve the effect of data fitting, a method of extremum extension on endpoints is proposed, which is aimed at the endpoint effect of empirical mode decomposition (EMD) algorithm. The improved EMD algorithm is used to decompose the observed signals, and then the decomposed components together with the prior observed signals are regarded as new observed signals. Thus the underdetermined situation is changed into an overdetermined case. Finally, we use independent component analysis (ICA) algorithm to obtain the estimation of source signals. Simulation result shows that the proposed algorithm is effective.
Key words: underdetermined blind source separationEMD(empirical mode decomposition)endpoint effectextremum extensionICA(independent component analysis)
盲源分离(blind source separation, BSS)是信号处理领域的一个研究热点, 在生物医学工程、地震检测、信号增强、雷达和声纳等领域都有极其广泛的应用, 而欠定盲分离因其观测信号数目小于源信号数目, 是盲分离中的难点, 也更符合实际问题.
目前, 欠定盲分离算法大体分为两类, 一类为基于信号的稀疏性的方法, 另一类为基于源信号的统计特性的方法, 这两类方法都是以信号稀疏性为前提的^[1].为解决欠定盲分离问题, 研究人员提出了一种利用经验模态分解(empirical mode decomposition, EMD)方法将欠定混合信号升维的思想.EMD方法是1998年Norden E.Huang等人提出的一种用于非线性、非平稳信号分解的方法, 目前大多用于机械故障诊断、管道泄漏检测等领域^[2-3].由于EMD算法具有的端点效应会使分解结果失真, 严重影响算法的精度, 因此一些新的信号延拓技术被相继提出.文献[4]提出一种镜像延拓法, 将信号对称地延拓为一个闭合环形信号来避免端点效应, 但是由于延拓后的信号趋势与原信号不一致, 导致EMD分解效果不好.还有较常用的神经网络延拓方法^[5], 该方法分解结果十分理想, 但是速度太慢且编程较复杂.
针对上述问题, 本文在多项式拟合算法的基础上提出了一种新的端点极值延拓方法, 根据信号序列极值点的规律求出端点处的近似值作为极值点估计, 与序列中的最大值(或最小值)作比较, 避免了延拓的极值不符合插值要求的情况, 最后利用独立成分分析(independent component analysis, ICA)算法得到源信号的估计.通过与镜像延拓EMD-ICA算法的仿真对比证明了本文提出的改进算法应用于欠定混合信号盲分离是可行、有效的.
1 欠定盲源分离及EMD端点效应分析1.1 欠定盲源分离盲源分离是指在源信号和传输信道参数未知的情况下, 仅由观测信号恢复源信号的过程, 其数学模型为

(1)

式中:s(t)=[s₁(t), s₂(t), …, s_n(t)]为n维独立源信号; x(t)=[x₁(t), x₂(t), …, x_m(t)]为m维观测信号; A为m×n维混合矩阵.在欠定条件下, 混合系统是不可逆的, 不能通过传统的盲分离算法来估计混合矩阵, 即使估计出了混合矩阵, 源信号的估计值也不是唯一的^[6]; 因此, 传统算法不能处理这种“病态”盲分离情况.本文从处理欠定盲分离问题入手, 利用改进EMD算法将欠定问题转化为超定问题, 再结合ICA算法实现混合信号的分离.
1.2 EMD算法及端点效应依据信号的时间尺度特征, EMD通过把复杂的信号分解为有限个不同时间尺度的本征模态函数分量(intrinsic mode function, IMF)来实现对信号不同模态的分离.IMF分量须满足如下两个性质^[7-8]:
① 整个数据长度中, 极值点个数和过零点个数必须相等或最多相差一个;
② 任何数据点处, 局部极大值的上包络和局部极小值的下包络的均值必须为零.
图 1给出了EMD算法的流程图.
图 1(Fig. 1)

图 1 EMD算法流程图Fig.1 Flow chart of EMD algorithm

在EMD分解步骤中, 通过对原始数据中极大值点和极小值点分别进行样条插值拟合，再取平均值来求包络平均，是最重要的一步.因为要分析的信号长度有限, 且不能确定信号的两个端点是否为极值点, 因此会在样条插值时产生数据的拟合误差^[9].由于除第一个IMF分量外的每一个分量都是以原始数据减去前一个IMF分量得到的残余项作为原始数据进行分解得到的, 因此只要有一个分量产生误差, 随着分解的不断进行, 误差就会从端点处向内依次传播, 使分解得到的每个信号分量都存在一定误差, 接下来利用ICA进行信号分离时就会影响源信号的分离精度^[10].
2 基于端点极值延拓的欠定盲源分离基于1.2节的误差分析, 本文提出一种新的端点极值延拓算法来抑制EMD的端点效应, 提高信号分离精度.算法的求解步骤如下.
1) 将端点处函数值与靠近端点的第一个极值作对比, 前者大, 则端点处的函数值可能为极大值点, 故进行极大值序列延拓; 反之, 进行极小值序列延拓.
2) 根据上一步的判断结果, 在极值点序列(若进行极大值序列延拓, 则取极大值序列; 进行极小值序列延拓, 取极小值序列)中取最靠近端点处的3个极值点(若所要取的极值点序列中极值点个数小于3个, 则取序列中所有元素), 进行多项式拟合, 求出多项式在端点处的函数值, 把此函数值看作该端点处近似极值点, 即实现极值延拓.其中求解拟合多项式的步骤如下.
① 确定拟合多项式的次数n.
② 计算出S_r和t_r:

(2)

式中:0, m分别表示离散数据起、终点的下标; x_i, y_i分别为离散数据点的横、纵坐标；r=0, 1, …, n.
③ 写出正规方程组:

(3)

④ 解正规方程组, 求出a₀, a₁, …, a_n.

(4)

3) 找出极大值序列中的最大值P_max, 极小值序列中的最小值P_min.若延拓的极大值小于P_min, 将该延拓极大值定义为(端点值+P_max)/2;若延拓的极小值大于P_max, 将该延拓极大值定义为(端点值+P_min)/2, 这样可避免延拓极值不符合插值要求的情况.
4) 利用三次样条函数插值新的极值点序列获得上下包络线.这样通过三次样条函数插值时在端点处就做到了有值可依, 避免上下包络线摆动的情况.
利用改进EMD算法对观测信号进行分解, 所得各信号分量和原观测信号构成新的观测信号, 再对新的观测信号进行ICA处理, 得到源信号的估计.
3 仿真实验与算法分析3.1 性能指标1) 相似系数:描述估计信号与源信号之间相似性的一个参数^[11].设源信号s中的第i个元素为s_i, s_i经过分离后的估计信号为s_i′, 则s_i与s_i′之间的相似系数为

(5)

式中cov(·)表示方差.
由每对估计信号与源信号相似系数组成相似系数矩阵, 若每行每列中的元素有且仅有一个接近于1, 其他都接近于0, 则可认为能实现分离效果.
2) 输出信噪比:描述估计信号s_i′与某个源信号s_i之间的对应关系, 计算公式为

(6)

估计信号与源信号之间的输出信噪比越大, 说明该信号分离效果越好^[11].
3.2 仿真实验3.2.1 正弦信号、方波信号与随机噪声实验首先随机选取三路源信号：s₁=2sin(20πt), s₂=2square(100 t, 50)，一路随机噪声信号(不要求0均值信号).其中采样频率f_s=1 000 Hz, 采样周期T=2.图 2为源信号的波形图.将三路源信号随机生成的2×3矩阵A混合,
图 2(Fig. 2)

图 2 三路源信号Fig.2 Three-channel signals

得到两路观测信号, 如图 3所示.分别用文献[4]给出的镜像延拓EMD算法和本文的EMD算法对观测信号进行分解, 进一步利用ICA算法对重新构成的观测信号进行分离，结果如图 4所示.
图 3(Fig. 3)

图 3 混合观测信号Fig.3 Mixed observation signals

图 4(Fig. 4)

图 4 两种算法的分离结果Fig.4 Separation results of two algorithms (a)—镜像延拓EMD-ICA；(b)—端点极值延拓EMD-ICA.

两种算法的相似系数如表 1所示, 其中s₁, s₂, s₃表示源信号, z₁, z₂, z₃表示分离信号.两种算法的输出信噪比见表 2.
表 1(Table 1)

表 1 两种算法的相似系数Table 1 Similarity coefficients of two algorithms 算法 相似系数 s1 s2 s3 镜像延拓EMD-ICA z1 0.994 3 0.011 2 0.031 1 z2 0.002 6 0.038 5 0.796 0 z3 0.014 2 0.997 2 0.000 7 端点极值延拓EMD-ICA z1 0.995 4 0.004 9 0.013 4 z2 0.010 6 0.999 9 0.040 7 z3 0.004 3 0.014 7 0.842 0

表 1 两种算法的相似系数 Table 1 Similarity coefficients of two algorithms

表 2(Table 2)

表 2 两种算法的输出信噪比Table 2 Output SNR of the two algorithms 算法 输入信号 输出信噪比 镜像延拓EMD-ICA s1 52.849 2 s2 69.657 2 s3 51.402 9 端点极值延拓EMD-ICA s1 54.376 0 s2 71.043 8 s3 53.293 8

表 2 两种算法的输出信噪比 Table 2 Output SNR of the two algorithms

由表 1可以看出, 两个相似系数矩阵中每行每列元素有且仅有一个接近1, 其他都接近0, 说明都实现了分离, 证明了算法的可行性.通过对比两个表中的数据可以看出, 端点极值延拓EMD-ICA算法分离信号与源信号的相似系数比镜像延拓EMD-ICA算法更接近于1, 说明本文提出的端点极值延拓EMD-ICA算法提高了估计信号与源信号的相似性, 分离效果明显改善.比如表 1中, z₂是源信号s₃的估计信号, 相似系数为0.796 0;对应表 2中z₃为源信号s₃的估计信号, 相似系数为0.842 0.由表 2可看出, 端点极值延拓EMD-ICA算法分离三路信号的输出信噪比大于镜像延拓EMD-ICA算法, 说明本文算法提高了分离精度.该实验共随机选取参数50次, 每次实验中由本文提出的端点极值延拓EMD-ICA算法分离的三路信号对应的相似系数及输出信噪比均比镜像延拓EMD-ICA算法要高, 进一步证明了算法的可行性、有效性.
3.2.2 语音信号实验选用三路语音信号, 音频采样大小为16位.采样率为44 kHz, 读入点数为30 000, 三路语音信号如图 5所示.随机选取一个2×3矩阵
图 5(Fig. 5)

图 5 三路原始声音信号时域波形及频谱图Fig.5 Time domain waveforms and spectrograms of three primitive sound signals (a)—时域波形；(b)—频谱图.

与语音信号构成的矩阵进行混合, 得到两路混合观测信号如图 6所示.利用镜像延拓EMD-ICA分离混合语音信号得到的解混声音信号如图 7所示, 本文改进的EMD-ICA算法得到的解混声音信号如图 8所示.
图 6(Fig. 6)

图 6 混合声音信号及其频谱图Fig.6 Mixed sound signals and its spectrograms (a)—时域波形;(b)—频谱图

图 7(Fig. 7)

图 7 镜像延拓EMD-ICA算法分离的声音信号及频谱图Fig.7 Separated sound signals and spectrograms of EMD-ICA algorithm for image extension (a)—时域波形；(b)—频谱图.

图 8(Fig. 8)

图 8 端点极值延拓EMD-ICA算法分离声音信号及频谱图Fig.8 Separated sound signals and spectrograms of EMD-ICA algorithm for extremum extension on endpoints (a)—时域波形;(b)—频谱图

两种算法的相似系数如表 3所示, 其中s₁, s₂, s₃表示原三路语音信号；y₁, y₂, y₃表示分离的语音信号.
表 3(Table 3)

表 3 两种算法的相似系数Table 3 Similarity coefficients of two algorithms 算法 相似系数 s1 s2 s3 镜像延拓EMD-ICA y1 0.587 2 0.014 9 0.110 4 y2 0.000 7 0.968 3 0.027 3 y3 0.014 4 0.005 5 0.969 8 端点极值延拓EMD-ICA y1 0.002 1 0.011 9 0.999 6 y2 0.038 3 0.989 6 0.014 0 y3 0.999 3 0.002 5 0.001 8

表 3 两种算法的相似系数 Table 3 Similarity coefficients of two algorithms

由表 3可以看出, 端点极值延拓EMD-ICA算法分离的三路语音信号与源信号相似系数比镜像延拓EMD-ICA算法分离语音信号与源信号相似系数都增大了, 更接近于1, 说明新算法提高了分离精度, 能更好地实现欠定语音信号的分离.由表 4可看出, 端点极值延拓EMD-ICA算法分离语音信号后的输出信噪比大于镜像延拓EMD-ICA算法, 说明分离效果有一定改善.
表 4(Table 4)

表 4 两种算法的输出信噪比Table 4 Output SNR of the two algorithms 算法 输入信号 输出信噪比 镜像延拓EMD-ICA s1 52.590 2 s2 51.402 9 s3 69.657 2 端点极值延拓EMD-ICA s1 54.376 0 s2 53.293 8 s3 71.043 8

表 4 两种算法的输出信噪比 Table 4 Output SNR of the two algorithms

4 结语本文通过分析EMD算法存在的端点效应, 提出一种端点极值延拓改进方案, 利用改进的EMD算法对混合观测信号进行分解, 构造新的观测信号, 将欠定问题转化为超定问题, 进一步利用ICA算法进行信号分离, 得到源信号的估计.将本文提出的改进算法和镜像延拓改进算法通过仿真实验作了对比分析, 可以看出本文算法改善了EMD中数据拟合的效果, 从而进一步提高了信号的分离精度.
参考文献

[1]	Tichavsky P, Koldovsky Z. Weight adjusted tensor method for blind separation of underdetermined mixtures of nonstationary sources[J].IEEE Transactions on Signal Processing, 2011, 59(3): 1037–1047.DOI:10.1109/TSP.2010.2096221
[2]	李志农, 吕亚平, 范涛, 等. 基于经验模态分解的机械故障欠定盲源分离方法[J].航空动力学报, 2009, 24(8): 1886–1892. ( Li Zhi-nong, Lyu Ya-ping, Fan Tao, et al. Underdetermined blind source separation method for mechanical fault based on empirical mode decomposition[J].Journal of Aerospace Power, 2009, 24(8): 1886–1892.)
[3]	孙洁娣, 郝雅立, 温江涛, 等. 基于EMD的高压燃气管道泄漏信号欠定盲分离方法[J].振动与冲击, 2013, 32(18): 81–86. ( Sun Jie-di, Hao Ya-li, Wen Jiang-tao, et al. Underdetermined blind source separation method of pipeline leakage signals based on empirical mode decomposition[J].Journal of Vibration & Shock, 2013, 32(18): 81–86.DOI:10.3969/j.issn.1000-3835.2013.18.015)
[4]	Zhu X L, Liu J P, Zhang X N. Analysis of motor imagery EEG based on Hilbert-Huang Transform[J].Advanced Materials Research, 2014, 998/999: 833–837.DOI:10.4028/www.scientific.net/AMR.998-999
[5]	孟宗, 顾海燕, 李姗姗. 基于神经网络集成的B样条经验模态分解端点效应抑制方法[J].机械工程学报, 2013, 49(4): 106–112. ( Meng Zong, Gu Hai-yan, Li Shan-shan. An empirical mode decomposition method based on neural network ensemble for endpoint effect suppression of B spline[J].Journal of Mechanical Engineering, 2013, 49(4): 106–112.)
[6]	Zhang Y, Wu K, Tan G, et al.An online adaptive algorithm for underdetermined blind source separation[C]// 2014 12th International Conference on Signal Processing.Hangzhou, 2014: 467-472.http://www.researchgate.net/publication/282374329_An_online_adaptive_algorithm_for_underdetermined_blind_source_separation
[7]	Wang J, Liu Y, Chao Z, et al.A modified single-channel blind separation method using EMD and ICA[C]// International Conference on Trustworthy Computing and Services.Berlin: Springer, 2014: 78-85.http://rd.springer.com/chapter/10.1007/978-3-662-43908-1_10
[8]	Gu Q W, Jin W D, Yu Z B. Blind source separation of single-channel train signal based on EEMD and ICA[J].Application Research of Computers, 2014, 31(5): 1551–1553.
[9]	齐扬阳, 于淼. 基于EMD的单通道盲源分离跳频通信抗干扰方法[J].计算机科学, 2016, 43(1): 149–153. ( Qi Yang-yang, Yu Miao. Anti-jamming method for frequency hopping communication based on single channel BSS and EMD[J].Computer Science, 2016, 43(1): 149–153.)
[10]	王传菲, 安钢, 王凯, 等. 基于镜像延拓和神经网络的EMD端点效应改进方法[J].装甲兵工程学院学报, 2010, 24(2): 66–69. ( Wang Chuan-fei, An Gang, Wang Kai, et al. An improved method of EMD endpoint effect based on mirror extension and neural network[J].Journal of the Academy of Armored Forces Engineering, 2010, 24(2): 66–69.DOI:10.3969/j.issn.1672-1497.2010.02.015)
[11]	杨杰, 俞文文, 田昊, 等. 基于独立分量分析的欠定盲源分离方法[J].振动与冲击, 2013, 32(7): 30–33. ( Yang Jie, Yu Wen-wen, Tian Hao, et al. Underdetermined blind source separation method based on independent component analysis[J].Journal of Vibration & Shock, 2013, 32(7): 30–33.DOI:10.3969/j.issn.1000-3835.2013.07.006)