Blind Source Separation with Python

This post is an overview of the article Blind source separation by multiresolution analysis using AMUSE algorithm, but here the goal is a implementation in Python.

Algorithms for blind source separation have been extensively studied in the last years. This paper proposes the use of multiresolution analysis in three decomposition levels of the wavelet transform, such as a preprocessing step, and the AMUSE algorithm to separate the source signals in distinct levels of resolution.

AMUSE algorithm explores the temporal structure of sources projecting the observation vector in an orthogonal space. In this space, the n largest singular values of the covariance matrix of the observation vector are distinct.

Let $s(t)$ a stochastic signal where $s(t) = [s_{1}(t), ..., s_{n}(t)]$ and each $s_{i}(t) (i = 1, ..., n)$ is named as "sources".

Exists a process which mixture this sources building a new signal $x(t)$. Whether consider a additive random noise $n(t)$ then $x(t)$ can be written: $x(t) = A s(t) + n(t)$.

How to get the sources $s(t)$ knowing only $x(t)$ , without information about the mixture process, $A$, and the noise $n(t)$?

The AMUSE algorithm solve the problem!

1. Estimate the covariance matrix of the $x(t)$: $R_x = E[x(t) x(t)^T]$.

2. Decompose the matrix $R_x$ by singular value decomposition (SVD).

3. Perform the orthogonal transformation of the data, by product: $y(t) = B x(t)$ where $B = (\table ψ_1, ... , 0; 0, ⋱ , 0; 0, ..., ψ_n)$ and $v_i$ is $i$-th singular value.

4. For a time delay $k$ compute the covariance matrix of the $R_y (k) = [R_y (k) + R_y (k)^T]/2$. The eigenvectors are stored in the columns of the $V$ matrix.

5. Estimate the sources: $ŝ = V^T B x(t)$.

The script Python of the AMUSE algorithm is show below.

#!/usr/bin/python

# -*- coding: UTF-8 -*-

import numpy as np

class AMUSE:

    def __init__(self, x, n_sources, tau):

        self.x = x

        self.n_sources = n_sources

        self.tau = tau

        self.__calc()

        

    def __calc(self):

       #BSS using eigenvalue value decomposition

       #Program written by A. Cichocki and R. Szupiluk at MATLAB
      

        R, N = self.x.shape

        Rxx = np.cov(self.x)

        U, S, U = np.linalg.svd(Rxx)

        

        if R > self.n_sources:

            noise_var = np.sum(self.x[self.n_sources+1:R+1])/(R - (self.n_sources + 1) + 1)

        else:

            noise_var = 0

        

        h = U[:,0:self.n_sources]

        T = np.zeros((R, self.n_sources))

        

        for m in range(0, self.n_sources):

            T[:, m] = np.dot((S[m] - noise_var)**(-0.5) ,  h[:,m])

        

        T = T.T

        y = np.dot(T, self.x)

        R1, N1 = y.shape

        Ryy = np.dot(y ,  np.hstack((np.zeros((R1, self.tau)), y[:,0:N1 - self.tau])).T) / N1

        Ryy = (Ryy + Ryy.T)/2

        D, B  = np.linalg.eig(Ryy)

        

        self.W = np.dot(B.T, T)

        self.sources = np.dot(self.W, self.x)

Example of use of the class AMUSE.

#!/usr/bin/env python

# -*- coding: utf-8 -*-

import sys

sys.path.append('/home/bruno/PESQUISAS/Python/') #directory of the AMUSE class

import numpy as np

import pylab as pl

from AMUSE import AMUSE

from scipy.io import wavfile

    

#read sources

fs, s1 = wavfile.read('/home/bruno/MESTRADO/Pesquisa/ArquivosAudio/drums.wav')

fs, s2 = wavfile.read('/home/bruno/MESTRADO/Pesquisa/ArquivosAudio/bass.wav')

#5 sec of sources

s1 = s1[0:220500]

s2 = s2[0:220500]

s = np.c_[s1, s2].T

A = np.array([[10.2,12.5],[5.7,2.4]])  #Matrix mixture. Merges the sources

x = np.dot(A, s) #observed signal

amuse = AMUSE(x, 2, 1)

s_hat = amuse.sources #estimate sources. 1/8 size of the original signal

W = amuse.W #separation matrix

The article Blind source separation by multiresolution analysis using AMUSE algorithm addresses the multiresolution analysis of wavelet transform.

The DWT (Discrete Wavelet Transform) decomposes a signal $x(t)$ in sub-bands with distinct frequencies. Every new level of resolution less noise appear in the approximation signal. In this way, a signal decomposed by DWT and then apply the AMUSE algorithm should improve the estimation.

Caution! The decomposition level must be chosen such that it does not change the waveform of the signal, otherwise the estimation fails.

Repeating the last example more now with applying the DWT.

Example of use of the class AMUSE with DWT.

#!/usr/bin/env python

# -*- coding: utf-8 -*-

import sys

sys.path.append('/home/bruno/PESQUISAS/Python/') #directory of the AMUSE class

import numpy as np

import pylab as pl

from AMUSE import AMUSE

import pywt

from scipy.io import wavfile

    

#read sources

fs, s1 = wavfile.read('/home/bruno/MESTRADO/Pesquisa/ArquivosAudio/drums.wav')

fs, s2 = wavfile.read('/home/bruno/MESTRADO/Pesquisa/ArquivosAudio/bass.wav')

#5 sec of sources

s1 = s1[0:220500]

s2 = s2[0:220500]

s = np.c_[s1, s2].T

A = np.array([[10.2,12.5],[5.7,2.4]])  #Matrix mixture. Merges the sources

x = np.dot(A, s) #observed signal

Wavelet = pywt.Wavelet('db2') #set wavelet

level_max = 3 #set level decomposition

#compute the coefficients

#the aproximation signal is given by coeffs1[0] and coeffs2[0]

coeffs1 = pywt.wavedec(x[0], Wavelet, level=level_max)

coeffs2 = pywt.wavedec(x[1], Wavelet, level=level_max)

x_ = np.c_[coeffs1[0], coeffs2[0]].T

amuse = AMUSE(x_, 2, 1)

s_hat_dwt = amuse.sources #estimate sources. 1/8 size of the original signal

W = amuse.W #separation matrix

s_hat = np.dot(W, x) #estimate sources. Same size of the original signal

Note: coeffs1[0] and coeffs2[0] has less samples which the original signal (1/8 of the samples), nevertheless the estimation will succeed because the waveform has not changed.

To show the results for both examples:

t = np.arange(0, 5, 1./fs)

fig0 = pl.figure()

pl.subplot(211)

pl.title('Original Sources')

pl.plot(t, s1)

pl.subplot(212)

pl.plot(t, s2)

pl.show()

fig1 = pl.figure()

pl.subplot(211)

pl.title('Observed Signal')

pl.plot(t, x[0])

pl.subplot(212)

pl.plot(t, x[1])

pl.show()

fig2 = pl.figure()

pl.subplot(211)

pl.title('Estimated Sources')

pl.plot(t, s_hat[0])

pl.subplot(212)

pl.plot(t, s_hat[1])

for fig in [fig0, fig1, fig2]:   

    for ax in fig.get_axes():

        ax.titlesize = 10

        ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))

        ax.set_xlabel('Time (sec)')

        ax.set_ylabel('Amplitude')

        ax.grid()

            

pl.show()

It should be noted that the amplitude phase estimation and order are ambiguous. This is normal. Is a characteristic common of the algorithms of blind source separation.

References

CHICHOCKI, A.; AMARI, S. Adaptive Blind Signal and Image Processing: Learning Algorithms and Applications. West Sussex: Wiley, 2003.

COMON, P.; JUTTEN, C. Handbook of Blind Source Separation: Independent Component Analysis and Applications. Oxford: Elsevier, 2010.

DAUBECHIES, I. Ten Lectures on Wavelets. CBMS-NSF regional conference series in applied mathematics, Philadelphia, 1992.

HYVÄRINEN, A.; KARHUNEN, J.; OJA, E. Independent Component Analysis. New York: Wiley, 2001.

HUANG, R.; CHEUNG, Y.; ZHU, S. A Parallel Architecture Using Discrete Wavelet Transform for Fast ICA Implementation. IEEE Int. Conf. Neural Networks & Signal Processing on, v. 14-17, p. 1358-1361, 2003.

KAWAGUCHI, A.; TRUONG, Y. K.; HUANG, X. Application of Polynomial Spline Independent Component Analysis to fMRI Data. In: NAIK, G. R. Independent Component Analysis for Audio and Biosignal Applications. Rijeka: InTech, 2012. p. 209-220.

LEO, M.; D'ORAZIO, T.; DISTANTE, A. Feature extraction for automatic ball recognition: comparison between wavelet and ICA preprocessing. Image and Signal Processing and Analysis, ISPA 2003. Proceedings of the 3rd International Symposium on, vol. 2, p. 587-592, 2003.

LÓ, P. M. G.; LOZANO, H. M.; SÁNCHEZ, F. L. P.; MORENO, L. N. O. Blind Source Separation of audio signals using independent component analysis and wavelets. Electrical Communications and Computers CONIELECOMP 21st International Conference on, p.152-157, 2011.

MISSAOUI, I.; LACHIRI, Z. Blind speech separation based on undecimated wavelet packet-perceptual filterbanks and independent component analysis. IJCSI International Journal of Computer Science Issues on, v. 8, 2011.

MIJOVIC, B.; VOS, M. D.; GLIGORIJEVIC, I.; TAELMAN, J.; HUFFEL, S. V. Using Wavelet Transformation in Blind Sources Separation of the Fetal Electrocardiogram. Majlesi Journal of Electrical Engineering on, v. 5, p. 2188- 2196, 2011.

TALBI, M.; AICHA, A. B.; SALHI, L.; CHERIF, A. Bionic Wavelet Based Denoising Using Source Separation. Int J Comput Commun on, v. 7, p. 574-585, 2012.

RADU, M.; HULLE, M. M. V. A comparative survey on adaptive neural network algorithms for independent component analysis. Romanian Reports in Physics on, v. 55.1, p. 49-74, 2003.

SHAYESTEH, M.; FALLAHIAN, J. Source Separation From Single-Channel Recordings by Combining Empirical-Mode Decomposition and Independent Component Analysis. IEEE Transactions on Biomedical Enginnering on, v. 57, p. 33-37, 2010.

STRANG, G.; NGUYEN, T. Wavelets and Filter Banks. Cambridge: Wellesley, 1997. TONG, L.; LIU, R.; SOON, V. C.; HUANG, Y. Indeterminacy and identifiability of blind identification. Circuits and Systems, IEEE Transactions on, v. 38, p. 499-509, 1991.

YEREDOR, A. Second-order methods based on color. In: COMON, P.; JUTTEN, C. Handbook of Blind Source Separation: Independent Component Analysis and Applications. Oxford: Elsevier, 2010. p. 227-278. WEEKS , M.

Digital Signal Processing using MATLAB and Wavelets. Massachusetts: Jones and Bartlett Publishers, 2007.


Links to paper:

1. Blind source separation by multiresolution analysis using AMUSE algorithm
2. Blind source separation by multiresolution analysis using AMUSE algorithm
3. Blind source separation by multiresolution analysis using AMUSE algorithm