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

本文链接：https://blog.csdn.net/weixin_47750895/article/details/115680198

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

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

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

二、FastICA算法步骤

1、步骤

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

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

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

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

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

5）对按列取向量,重组为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）设一初始零列向量，重组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、白化矩阵估计

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

plt.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.5
    w = 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 = w
        gx = 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 += 1
        if (count == Maxcount):
            print("reach Maxcount，exit loop", 1-LA.norm(absCos, 1))
            break
    print("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峰提取计算心率以及相关特征值提取。

四、小结

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