对于鸡尾酒会问题,一种简单的情况如下:有n个人在同时说话,同时又m个声音接收器捕捉到了信号之间的线性组合,于是我们可以得到m组声音数据。那么,如何利用这m组接收到的声音信号恢复成原来的n组独立信号呢?
在上个世纪末,鸡尾酒会问题催生了各种盲源分离问题,也有不少机器学习算法被应用于此问题。其中,独立成分分析是简便且有效的一种。
##问题分析 设原始信号为S1、S2、...、Sn,一定存在一个mn的权重矩阵A,使得每个声音接收器对应A中一行。如m=n=2时,令A=[[0.6, 0.4], [0.45, 0.55]],则第一个麦克风接收到的信号为x1=0.6S1+0.4S2,第二个麦克风接收到的信号为x2=0.45S1+0.55S2。即$$AS=x$$ 将该公式左右同乘以A的逆矩阵$$W=A^{-1}$$,得$$A^{-1}AS=A^{-1}x$$,即$$S=Wx$$,这样问题就解决了(当m≠n时,A不是方阵,此时可以乘A的伪逆)。于是问题来了:如何求A的逆矩阵W?
##独立成分分析(Independent Component Analysis) ICA算法始于Bell和Sejnowski,下面用极大似然估计简单推导一下ICA的求解过程:
###引理:已知x的概率分布$$f_{x}(x)$$,若$$Y=Ax$$,则对其累积分布函数有$$F_{Y}(y)=P{Y\leq y}=P{Ax\leq y}=P{x\leq y/A}=F_{X}(y/A)$$。求导得$$f_{Y}(y)=frac{f_{X}(x/A)}{|A|}$$。
每个信号源$$S_{i}$$都有其概率密度$$p_{i}(s)$$,在第j时刻,向量s的联合概率密度为$p(s_{j})=\prod_{i=1}^{n} p_{i}(s_{i,j})$。通过引理,我们可以将该式转化为$$L(x_{j})=|W|p_{i}(Wx_{j})=|W|\prod_{i=1}^{n} p_{i}(w_{i}^{T}x_{j})$$。 于是可以得到似然函数$$L(x)=\prod_{j=1}^{m} p(x_{j})=\prod_{j=1}^{m} [|W|\prod_{i=1}^{n} p_{i}(w_{i}^{T}x_{j})]$$ 现在似然函数还是个连乘的形式,不便计算。我们将其取对数,得到对数似然函数$$l(x)=\sum_{j=1}^{m} [\sum_{i=1}^{n} ln p_{i}(w_{i}^{T}x_{j}) + ln|W|]$$ 这种形式的似然函数就可以当做我们的成本函数cost function了,即$$J(W)=\sum_{j=1}^{m} [\sum_{i=1}^{n} ln p_{i}(w_{i}^{T}x_{j}) + ln|W|]$$ 但是,式中的$$p_{i}$$,即先验概率$$p(s)$$是未知的,但是我们可以假定它服从某种已知的分布。一般的,我们可以假设$$p(s)$$的累积分布函数$$F(s)$$是sigmoid函数。当然了,假设它们服从正态分布(erf函数),甚至可以使用tanh、arctan当做s的累积分布函数,这些先验假设都是可行的。 以sigmoid函数为例,令$$g(s)=frac{1}{1+e^{-s}}$$,易知$$g'(s)=g(s)(1-g(s))$$、$$g''(s)=g'(s)(1-2g(s))$$。带入对数似然函数得$$J(W)=\sum_{i=1}^{m}[\sum_{j=1}^{n} ln g'(wx) + ln|W|]$$。 将对数似然函数对$$W$$求偏导,得$$\frac{\partial J(W)}{\partial W}=[1-2g(Wx)]x + (W^{-1})^{T}$$。 到了这一步,相信大家已经知道下面怎么做了:梯度上升到最大值,求出此时的$$W$$即可。也就是$$W += \alpha * (1-2g(W*x) + (W^{-1})^{T})$$。这时候一般梯度上升、随机梯度上升和mini-batch算法就可以派上用场了。 用python实现出来的代码如下:
def sigmoid(x):
return 1.0 / (1 + np.exp(-x))
def ica(x):
# x = np.array(x)
n = len(x[0])
m = len(x)
w = np.eye(n)
iw = np.zeros([n, n])
w1 = np.zeros([n, n])
alpha = 0.001
for time in range(200):
for i in range(m):
for j in range(n):
# 使用sigmoid函数作为CDF
t = 1 - 2 * sigmoid(np.dot(w[j], x[i]))
w1[j] = t * x[i]
# 使用arctan函数作为CDF
# t = -2 * np.dot(w[j], x[i]) / (1 + np.dot(w[j] ** 2, x[i] ** 2))
# w1[j] = t * x[i]
# 使用tanh函数作为CDF
# t = -2 * x[i] * np.tanh(np.dot(w[j], x[i]))
# w1[j] = t
# 使用erf函数(正态分布)的CDF
# t = -2 * np.dot(w[j], x[i])
# w1[j] = t * x[i]
iw = np.linalg.inv(w.transpose())
w1 += iw * alpha
w1 = w1.transpose() * alpha
w += w1
# print time, ":\t"
# print w
return w
##仿真结果 首先,使用python模拟,生成一个扫描波、一个正弦波如下;使用权重矩阵$$A=\begin{pmatrix} 0.6 & 0.4 \ 0.45 & 0.55 \end{pmatrix}\$$对原始信号进行加权,得到的混合信号如下:
原始信号:
混合信号:
import numpy as np
from matplotlib import pyplot as plt
def show(s1, s2, label1, label2):
x = [i for i in range(1000)]
plt.plot(x, s1, 'r', label = label1)
plt.hold(True)
plt.plot(x, s2, 'g', label = label2)
plt.hold()
plt.legend()
plt.show()
if __name__ == "__main__":
s1 = np.array([np.sin(float(x) / 20) for x in range(1000)])
s2 = np.array([float(x) / 50 for x in range(50)] * 20)
show(s1, s2, 'Signal 1', 'Signal 2')
A = np.array([[0.6, 0.4], [0.45, 0.55]])
x = np.dot(A, np.array([s1, s2])).transpose()
show(x[:,0], x[:,1], 'Mic 1', 'Mic 2')
# x -= x.mean()
w = ica(x)
[ps1, ps2] = np.dot(x, w).transpose()
show(ps1, ps2, 'Component 1', 'Component 2')
将混合信号x输入给ICA算法,并使用得到的W进行还原,得到的结果如下:
虽然有些局部失真,但得到这种结果已经很不错了。下面是使用其他CDF函数得到的结果:
陈政/arc001 原创作品转载请注明出处
我的博客即将搬运同步至腾讯云+社区,邀请大家一同入驻:https://cloud.tencent.com/developer/support-plan?invite_code=6klbojckkbuw