我只是好奇

AuxIVA Algorithm

发表于 分类于
本文字数： 3.7k 阅读时长 ≈ 3 分钟

这篇文章简单说一下AuxIVA算法，前段时间接触并实现了一下，有些结论和体会，在这里做一下小结。

BSS（Blind Source/Speech/Signal Separation）问题中根据信号混合的方式，可以分为instantaneous/convolutive两种，前者仅仅是简单的线性叠加，和convolutive不同，对信号源到麦克风之间的冲击响应不加以建模；如果根据信号源数目和观测（麦克风）数目划分，有determined和underdetermined两种，前者表示观测信道足够（大于等于信号源数目），后者对应数量不足的情况。因此，BSS问题可以组合成四种case，比较显而易见，underdetermined/convolutive的情况最为复杂。

BSS的解决方法，比较主流的是ICA[2]及其衍生算法，比如IVA[3]，TRINICON[4]等等。相比Frequency-ICA，IVA的优点是不需要解决frequency上的置换问题，显得十分简洁。

Problem setup

在频域，convolutive mixture情况下（$N$路信号和观测），观测信号$\mathbf{x}_{f,t}$可以写成源信号$\mathbf{s}_{f,t}$的线性叠加：

其中$ \boldsymbol{A}_f \in \mathbf{C}^{N \times N}$表示子带$f$上的混合矩阵。算法最终需要估计一个对应的解混矩阵$\boldsymbol{A}_f $，使得：

中的各路信号尽可能相互独立，对应的代价函数用KL散度定义：

其中$\mathbf{y}_t^n = [\mathbf{y}_{0,t}^n, \cdots, \mathbf{y}_{N_f,t}^n]^T$，$p$表示以频率为随机变量的分布密度函数，$q$为联合分布密度函数。$(1)$式最终可以简化成：

blind情况下，对于$p$一般是事先指定的。

Estimate Demixing Matrix

根据式$(1)$，只要估计出了$\boldsymbol{W} = {\boldsymbol{W}_f}_{f=0}^{N_f}$，做iSTFT就可以获得分离的时域信号：

对于AuxIVA而言，解混矩阵的估计通过构造辅助函数进行，原始版本是基于自然梯度的梯度下降法，学习率/步长的设置对于最终结果的影响较大，相比之，AuxIVA的收敛更快，结果更加稳定。

令$G(\mathbf{y}_t^n) =- \log p(\mathbf{y}_t^n)$，AuxIVA的更新法则为：

  1. update auxiliary variables

    其中$r_{n,t} = \Vert \mathbf{y}_t^n \Vert_2$，$G’(r_{n,t})$一般取1。

  2. update demixing matrix

    $\mathbf{w}_{n,f}$表示$\boldsymbol{W}_f$的第$n$行。

参考pyroomacoustics的实现，上述过程可以写成：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
def auxiva(X, epochs=20):
    """
    Arguments:
        X: shape in N x T x F
    Return
        Y: same shape as X
    """
    N, T, F = X.shape
    # X: F x T x N
    X = X.transpose([2, 1, 0])
    # F x N x N
    W = np.array([np.eye(N, dtype=np.complex) for f in range(F)])
    I = np.eye(N)
    # Y: F x T x N
    Y = np.einsum("...tn,...nx->...tx", X, np.conj(W))

    for _ in range(epochs):
        # T x N
        R = np.sqrt(np.sum(np.abs(Y)**2, axis=0))
        # N x T
        Gr = 1 / (R.T + EPSILON)
        for f in range(F):
            for n in range(N):
                # compute V
                V = (np.dot(np.expand_dims(Gr[n], 0) * X[f].T, np.conj(
                    X[f]))) / T
                # update W
                w = np.linalg.solve(np.conj(W[f].T) @ V, I[n])
                W[f, :, n] = w / np.inner(np.conj(w), V @ w)

        Y = np.einsum("...tn,...nx->...tx", X, np.conj(W))
    # F x T x N => N x T x F
    Y = np.transpose(Y, [2, 1, 0])
    return Y

实际表现，IVA的结果给我的感觉不是非常理想（但是相比原始信号，还是可以体会到显著的分离效果），毕竟是完全blind的方法。另外需要注意的是，直接对$(1)$式做iSTFT合成的时域信号在能量上和原始信号不匹配，因此，很多地方都提到了project back这个概念，计算一个scale矩阵，对分离信号的每个bin上的能量进行缩放。pyroomacoustics中的算法为：

ICA系列比较重要的问题就是如何衡量信号的独立性，依此设计cost function，希望深纠的同学可以详细看一下参考文献[2]中的一些基本概念，对于了解这系列算法挺有帮助。

Reference

[1]. Ono N. Stable and fast update rules for independent vector analysis based on auxiliary function technique[C]//Applications of Signal Processing to Audio and Acoustics (WASPAA), 2011 IEEE Workshop on. IEEE, 2011: 189-192.
[2]. Hyvärinen A, Oja E. Independent component analysis: algorithms and applications[J]. Neural networks, 2000, 13(4-5): 411-430.
[3]. Kim T, Eltoft T, Lee T W. Independent vector analysis: An extension of ICA to multivariate components[C]//International Conference on Independent Component Analysis and Signal Separation. Springer, Berlin, Heidelberg, 2006: 165-172.
[4]. Buchner H, Aichner R, Kellermann W. TRINICON: A versatile framework for multichannel blind signal processing[C]//Acoustics, Speech, and Signal Processing, 2004. Proceedings.(ICASSP’04). IEEE International Conference on. IEEE, 2004, 3: iii-889.