一、ICA原理

独立成分分析（Independent Component Analysis, ICA），来源于著名的“鸡尾酒会”事件。关于它的背景不再过多赘述，只是在一场所，有N个人同时说话，并且有N个话筒同时录音，如何从这N段录音中分离出原始N个人单独的语音呢？

对于我们所得到的观测信号，每一路都是由这几种独立源 $s\left ( t \right )$ 成分组成的，只是不同通道采出的观测信号中各项信源的权重系数不同。假设只有两个独立源 $s_{1}s_{2}$ ,分别代表母体心电和胎儿心电，由于电极片位置的不同，靠近母体心脏的电极采出的混合信号 $x_{1}$ 与靠近腹部的电极采出的混合信号 $x_{2}$ 中的独立源 $s_{1}s_{2}$ 的权重系数肯定是不同的，腹部胎儿心电所占权重要大。

而实际上，我们只知道观测信号 $X\left ( t \right )$ ,对于混合矩阵 $A$ 和信源 $S\left ( t \right )$ 的成分（个数）是未知的。但是我们已知要分离的信号中信源是相互独立的，是独立源。因此，我们可以利用源信号的独立性和非高斯性，得到一个混合矩阵（这里的混合矩阵不是 $A$ ，或者称之为解混矩阵 $B$ ），使之与源信号对原观测信号有着最佳逼近。求解独立源问题便转化为最优化问题，我们所要求的就是解混矩阵B。

独立成分分析的过程如下图所示：

原ICA算法经过不断改进，从最大熵、最小互信息、极大似然估计以及负熵最大等角度提出了一系列优化算法。其中最常用的就是FastICA优化算法，利用负熵最大判据，该算法具有简单、收敛快的优点。

二、FastICA算法步骤

1、步骤

FastICA的学习规则是找出一个方向使 $W.T * X$ 有最大的非高斯性。这一思想主要体现步骤6~9步。

1）对观测信号 $X$ 取均值，中心化处理；

2）白化矩阵估计，对中心化后的观测信号进行白化得到 $Z$ ；

3）估计要分离的独立源个数R；分离母体和胎儿心电，通常我们要用四通道来分离。

4）根据分量个数来确定初始化权重矩阵 $W$ 、最大迭代次数 $Maxcount$ 、收敛误差 $Critical$ ；

5）对 $W$ 按列取向量 $w$ ,重组为R行1列的矩阵；

6）令 $Wp=E{{ Z*g(Z.T * w) }} - E{ g'(Z.T * w) }*w$ ,其中非线性函数 $g(ax)$ ，可取 $tanh(x)$ 或 $x*exp(-x^{2}/2)$ 或 $x^{3}$ ,a通常取1；

7）设一初始零列向量 $wn$ ，重组R行1列， $w = Wp - W * W.T * w$ ；

8）找到一方向 $Cos = w.T * wn$ ，计算其1范记为norm；

9）判断 $1-norm<Critical$ ;

10）若收敛，则步骤7得到的 $w= W[:,n]$ ;若不收敛，则返回第六步迭代执行至收敛。

2、白化矩阵估计

一般情况下，我们所获得的观测信号都具有相关性，所以首先要对数据进行去相关处理，即白化/球化处理。白化后的观测信号去除了各通道信号间的相关性，极大简化了后续独立分量的提取过程，并且增强了算法的收敛能力；而在白化之前，对数据进行中心化处理，即去均值，其主要目的是是白化后的信号 $Z$ 的方差为1、协方差矩阵为单位矩阵。

如何估计白化矩阵？

利用主成分分析（PCA），使一变换 $W = \Lambda ^{1/2}U^{T}X$ 。其中

$\Lambda$ 为X的协方差矩阵的特征值构成的对角矩阵

U为特征值对应的特征向量构成的正交矩阵

使得 $W = WX = \Lambda ^{1/2}U^{T}X$ 。

如何证明Z白化？

三、Python实现FastICA胎心信号分离

'''
Maximum criterion of negative entropy
样本源：mat
fs=250Hz
可选channels：34
'''import matplotlib.pyplot as plt
import numpy as np
from numpy import linalg as LA
from pyemd import EMD, Visualisation
from scipy import signal
import wfdb
import scipy.io
'''
channels = [2:6]
samprate = 250HZ
sampcount = 2500
time = 10s
'''
record1 = scipy.io.loadmat('FetalDatabase/sub01_snr03dB_l2_c1_mecgm.mat')
record2 = scipy.io.loadmat('FetalDatabase/sub01_snr03dB_l2_c1_fecg1m.mat')
record3 = scipy.io.loadmat('FetalDatabase/sub01_snr03dB_l2_c1_noise1m.mat')
record4 = scipy.io.loadmat('FetalDatabase/sub01_snr03dB_l2_c1_noise2m.mat')
demo1 = record1['val'].astype(float)
demo2 = record2['val'].astype(float)
demo3 = record3['val'].astype(float)
demo4 = record4['val'].astype(float)
sampcount = 2500
samprate = 250
data = demo1+demo2+demo3+demo4
data = np.array([data[i][0:sampcount] for i in range(2, 6)])
t = np.linspace(0, samprate/samprate, sampcount)'''
归一化处理
Output：四路归一化数据mix
'''
for i in range(4):for index, val in enumerate(data[i]):x = float((val - np.min(data[i])) / (np.max(data[i]) - np.min(data[i])))data[i][index] = x
mix = np.array([data[0], data[1], data[2], data[3]])
R, C = mix.shapeplt.figure(1)
plt.subplot(411)
plt.plot(mix[0])
plt.title('Observed signal')
plt.subplot(412)
plt.plot(mix[1])
plt.subplot(413)
plt.plot(mix[2])
plt.subplot(414)
plt.plot(mix[3])'''
中心化/去均值化
Output：在同一基线上的观察信号mix
'''
average = np.mean(mix, axis=1)  # 计算行均值
for i in range(R):mix[i, :] -= average[i]'''
估计白化矩阵
Output：白化矩阵White、白化后的观测信号（正交）Z
'''
Cx = np.cov(mix)
value, eigvector = np.linalg.eig(Cx)  # 计算协方差阵的特征值
val = value ** (-1 / 2) * np.eye(R, dtype=float)
White = np.dot(val, eigvector.T)  # 白化矩阵
Z = np.dot(White, mix)'''
优化后的基于负熵最大判据的FastICA算法
1、初始化权重矩阵
2、非线性函数gx=tanh(ax)  系数a=1
3、最大迭代次数Maxcount = 1000
4、收敛判据Critical = 0.0001
Output：ICA分离的四路源信号
'''
W = np.zeros((R, R))
a = 1
Maxcount = 1000
Critical = 0.0001
for n in range(R):count = 0# 加速收敛w = W[:, n].reshape(R, 1) - 0.5w = w/LA.norm(w)wOld = np.zeros(w.shape)absCos = np.abs(np.dot(w.T, wOld))while (1-LA.norm(absCos, 1)) > Critical:wOld = wgx = np.tanh(np.dot(a*Z.T, w))e = np.mean(1-gx ** 2)w = np.dot(Z, gx)/C - np.dot(a*e, w)w = w - np.dot(W.dot(W.T), w)w = w / LA.norm(w)absCos = np.abs(np.dot(w.T, wOld))count += 1if (count == Maxcount):print("reach Maxcount，exit loop", 1-LA.norm(absCos, 1))breakprint("loop count:", count)W[:, n] = w.reshape(R,)
ICA = np.dot(W.T, Z)plt.figure(2)
plt.subplot(411)
plt.plot(ICA[0])
plt.title('Signal after separation by ICA')
plt.subplot(412)
plt.plot(ICA[1])
plt.subplot(413)
plt.plot(ICA[2])
plt.subplot(414)
plt.plot(ICA[3])
plt.show()

原始四路观测信号

经ICA分离后的独立源信号

单独拿出1、2两路信号

可见第一路为胎心信号，母体信号基本被剔除；第二路为母体信号。后续经过滤波处理后，可进一步增强胎心信号、抑制母体信号，进行R峰提取计算心率以及相关特征值提取。

四、小结

理想情况下，所有样本数据都应该参与计算，但实际上采出的数据并没有理想中那么完美，且大多数情况都是在某一段时间内的信号质量相对较好。现实情况下，我们要选取一部分样本来平均估计，并且样本点的个数对最后分离的结果有很大影响，理论上收敛不理想可以增加样本点数量；收敛误差并不是固定不变的，收敛误差的不同对同一个数据的结果影响也很大，并不是误差越小越好；不同样本、不同样本点的收敛误差设置也是不同的，可以在算法中对收敛误差进行寻优，通过样本熵对分离信号质量进行判断，在一定范围内找到最优的收敛误差。

Python实现基于负熵最大判据的FastICA胎心信号分离相关推荐

Python实现基于朴素贝叶斯的垃圾邮件分类标签： python朴素贝叶斯垃圾邮件分类 2016-04-20 15:09 2750人阅读评论(1) 收藏举报分类：机器学习（19）听说
Python实现基于朴素贝叶斯的垃圾邮件分类标签: python朴素贝叶斯垃圾邮件分类 2016-04-20 15:09 2750人阅读评论(1) 收藏举报分类: 机器学习(19) 听说朴 ...
Python构建基于elkan优化算法的K-Means聚类模型
Python构建基于elkan优化算法的K-Means聚类模型目录 Python构建基于elkan优化算法的K-Means聚类模型 #elkan优化算法
如何用 Python 进行基于深度学习的计算机视觉项目开发？
令人惊喜的"智能"年代深度学习有着广阔的前景我们正处在一个"智能"的年代,比如智能手机中的语音助手.机器翻译和人脸识别:战胜过日本将棋冠军.国际象棋冠军, ...
手把手教你python实现量价形态选股知乎_【手把手教你】Python实现基于事件驱动的量化回测...
01引言使用矢量化方法(pandas)建立的基于研究的量化回测框架,不考虑交易的委托成交行为,与真实市场情况差距比较大.今天为大家介绍的是基于事件驱动的回测框架,这是一种十分复杂的回测系统,力图模拟 ...
ICA独立成分分析—FastICA基于负熵最大
1. 概念官方解释:利用统计原理进行计算的方法,是一种线性变换. ICA分为基于信息论准则的迭代算法和基于统计学的代数方法两大类,如FastICA算法,Infomax算法,最大似然估计算法等. 这里 ...
postgresql 远程用户_构建Python pandas基于SSH远程MySQL和PostgreSQL的数据分析
背景知识视频教程 Python中使用Pandas教程 - 国外课栈viadean.com Pandas数据分析与探索 - 国外课栈viadean.com 如果您无法从外部环境直接访问数据库,则可能 ...
python数据框去重_【Python】基于某些列删除数据框中的重复值
Python按照某些列去重,可用drop_duplicates函数轻松处理.本文致力用简洁的语言介绍该函数. 一.drop_duplicates函数介绍 drop_duplicates函数可以按某列去 ...
python绘制星空图_【Python】基于某些列删除数据框中的重复值
阿黎逸阳精选Python.SQL.R.MATLAB等相关知识,让你的学习和工作更出彩(可提供风控建模干货经验). Python按照某些列去重 ,可用 drop_duplicates函数轻松处理 . ...
python实现蒙特卡洛算法_用Python实现基于蒙特卡洛算法小实验
用Python实现基于蒙特卡洛算法小实验蒙特卡洛算法思想蒙特卡洛(Monte Carlo)法是一类随机算法的统称,提出者是大名鼎鼎的数学家冯· 诺伊曼 ,他在20世纪40年代中期用驰名世界的赌城- ...
Python实现基于HDFS的云盘系统
Python实现基于HDFS的云盘系统一.云盘系统二.功能需求 2.1.用户管理 2.2.文件管理 2.3.界面设计三.用户代码 3.1 用户登录 3.2 用户注册 3.3 用户退出四.文件代 ...

Python实现基于负熵最大判据的FastICA胎心信号分离