Sunday, September 11, 2016

QRS detection in 12-lead ECG based Multiresolution Singular Value Decomposition

This is a summary of work "Detecção de Complexos QRS do ECG pela Decomposição em Valores Singulares em Multirresolução".

Multiresolution Singular Value Decomposition (MRSVD) 

The Singular Value Decomposition (SVD) is a method to decompose a matrix in three other matrices: two orthogonal and one diagonal. Let $X$ $m×n$ be a real-valued matrix, with $m ≤ n$, then the SVD of $X$ may be written [4]:

$X = U Γ V^{T}$ (1)

where $U$ $m×m$ and $V$ $n×m$ are orthogonal matrix and $Γ$ $m×m$ a diagonal matrix.

The columns of $U$, called left singular vectors, are eigenvectors of $XX^{T}$ and the columns of $V$, called right singular vectors, are eigenvectors of $X^{T}X$. The diagonal entries of $Γ$ are the singular values, square root of the eigenvalues of $XX^{T}$ or $X^{T}X$ and are ordered so that $γ 1 ≥ γ 2 ≥ · · ·≥0 $. The orthogonality of $U$ also allows to write equation (1) as $X̂ = ΓV T$ , where

$X̂ = U^{T}X$ (2)

The Principal Component Analysis (PCA) is a method related to SVD which decomposes a data matrix in uncorrelated principal components. The columns of the decomposition matrix are formed by the eigenvec-tors of the covariance matrix. Let $X$ be a data matrix and $X̄$ is the same matrix with the mean removed. The eigendecomposition of the scatter matrix, estimation of the covariance matrix, may be written $X̄X̄^{T} = UΛU^{T}$ [4], where the columns of $U$ are eigenvectors associated with eigenvalues of $Λ$. The principal components are given by equation $Y = U^{T}X̄$, which is the same as (2), except for the mean removed. This implies that, when the mean of $X$ is removed before the SVD, the principal components of $X̂$ are uncorrelated. This way, the first principal component $φ = X̂(1, ·)$ is the approximation signal, and the other $ψ = X̂(2, ·)$, the detail component, where $X̂(n, ·)$ is n-th row of matrix $X̂$. Relating to the structure of a two-channel filter bank, $U$ represent the filtering process and $X̂$ a sub-band decomposition of $X$ [4].

In order to implement a multiresolution SVD, for the one-dimensional case, [4] propose the construction of a signal $X = [x(1) · · · x(n)]$, being n divisible by $2^L$ for any integer $L ≥ 1$, as a data matrix with the first row containing the odd-numbered entries, and the other ones containing the even-numbered entries. Therefore, the first resolution level the data matrix becomes

$X_{1} = (\table x(1), ... , x(n-1); x(2), ... , x(n))$ (3)

The next step is the removing of the mean, multiplying $X_1$ by $H_1$ , be $H_l = I_l − ( 1/l )e_l e_{l}^{T}$ for the $l$th level, where $e_l$ is the $l × 1$ vector containing all unit entries. For the decorrelation filter $U$ the eigen-decomposition of the scatter matrix is computed. The components in subband are obtained by multiplying of this filter by the data matrix with the mean removed (equation 7).

The following equations define the MRSVD for $l = 1, · · · , L − 1$ levels, $n_l = n/2^l$ and $φ_0 = X$:

$X_{l} = (\table φ_{l-1} (1), ... , φ_{l-1} (n-1); φ_{l-1} (2), ... , φ_{l-1} (n))$ (4)
$ X̄ _{l}= X _{l} H_{n_l} $ (5)
$ X̄ _{l} X̄ _{l}^{T} = U_l Γ_{l}^{2} U_{l}^{T}$ (6)
$X̂_{l} = U_{l}^{T} X̄ _{l}$ (7)
$ φ_{l} = X̂_{l}(1, cdot)$ (8)
$ ψ_{l} = X̂_{l}(2, cdot)$ (9)

Proposed Method 

The sum this detail components provides a signal whose QRS complexes are in highlighted, however due the oversampling operation the location of the start point and end in that sum is not same of the QRS complex in the original signal, therefore an oversampling operation is need. Besides it is appropriate to calculate the absolute value of the oversampled detail components, to avoid the calculation of two thresholds.

The ECG signal after the pre-processing stage is:

$ψ_{sum} = |(↑2^{2})ψ_2| + |(↑2^{3})ψ_3|$ (11)

In the decision stage, the threshold is calculated by equation:

$thr = 0.3 max(ψ_{sum})$ (12)

and checked if the amplitude of $ψ_sum$, for each sample, is higher than $thr$. In the affirmative case the beginning of the QRS complex is found. After detecting the starting point, if the amplitude of $ψ_sum$ is less than $thr$, then is set the end of the QRS complex. However due the oversampling operation (equation (11)) have zeros between samples and this causes the algorithm finds the end of one QRS complex when in reality it is not the end.

The solution is using the physiological information of the heart to avoid the premature detection of end of the QRS complex, based on the fact which depolarization of the ventricles lasts from 0.06 to 0.10 (or 0.12) seconds [11]. Thus it cannot be found the end of the QRS complex before $0.06 Fs$ samples, where $Fs$ is the sampling rate. Let $qrs_ s$ and $qrs_e$ be the location of the starting samples and end, respectively, of each QRS complex, so the difference between them must satisfy the relation:

 $| qrs_s - qrs_e| >= 0.12 Fs$ (13)

It is clear that the limitation in equation (13) can returns an interval greater than the QRS interval, perhaps until ST segment. Nevertheless the maximum point it will still be a peak R, because in $ψ_sum$ the waves of low frequencies, waves T and P, have little energy, soon their amplitudes are attenuated. After found the intervals of the QRS complexes, a search is made, in the original signal, by the maximum point corresponding to the R peak.

The decision stage is executed to check if the R peak detected is authentic or not. The initial verification is relating the distance of the first R peak, which must occur after the depo-larization of the atria, i.e, after the P wave. Let $r(t)$ be the vector which stores the R peaks found, where $t = (t_1 , · · · , t_n )$ is the vector with the localization of all R peaks. For first R peak, the following equation must be satisfied:

$t_1 > 0.2 Fs$ (14)

When the algorithm has detected at least one R peak, it executes the next verifying based in the distance between two R peaks that is related to rate heart $RR = 60/HR$ [6], where $HR$ is a rate heart, measured in bpm (beats per minute) and $RR$ seconds. The maximum rate heart was obtained by experiments of [12] and is expressed by equation $HR_{max} =208−0.7a$, where $a$ is the patient age. Taking ${HR_max} = 208$, since there is no information about the age of the patients in selected database, the minimum samples amount between two R peaks is

$RR = 60/208 Fs$ (15)

Therefore, if RR, than the minimum, it should be chosen the one with the largest amplitude, removing the other which does not satisfy the criterium. Otherwise, both peaks are stored as correct. After detected all R peaks, its analyzed, in the correction stage, the need to search other R peaks. When exists a R peak 30% larger than the most of the peaks the detection fails because the threshold, equation (12), will be above the outhers. To circumvent this inconvenience a second threshold should be computed, removing the highest R peak and computing a new threshold for the resulting signal.

To decide when this calculation should be performed, it is considered the minimum amount of R peaks for any ECG signal, which is expressed by equation

$R_min = τ HR_min / 60$ (16)
where $τ$ is the signal duration, in seconds, and $HR_min = 35$ [5] is a minimum rate heart. The manner as this is done in practice, is by copying the original signal and replacing the of QRS interval by a null vector, and compute a new threshold.

Python scripts

 file MSVD.py

#!/usr/bin/python
# -*- coding: UTF-8 -*-

import numpy as np

class MSVD:
    
    def __init__(self, signal):
        
        self.signal = signal
    
    '''
    Decomposes the signal at a certain level
    If level_dec = 0, decompose in max level
    '''    
    def dec(self, level_dec, mean_removal = False):
    
        phi = signal
        
        if level_dec == 0:
            level_dec = self.levelmax(self.signal)
    
        for l in range(1, level_dec + 1):
            N = phi.size
            even = phi[0:N:2]
            odd = phi[1:N:2]
            X = np.vstack((even , odd))
            
            if (mean_removal):
                Xbar = X - X.mean()
            else:
                Xbar = X
            
            U, S, Vt = np.linalg.svd(Xbar)
            Xhat = np.dot(U.T , Xbar)
            phi = Xhat[0]
            psi = Xhat[1]
    
        return phi, psi, U
        
    '''
    Returns max level decomposition
    '''
    def levelmax(self, signal):
        return int(np.log(signal.size) / np.log(2)) - 1

------------------------------------------------------------------------------------------------------------------

file QRSdetect.py

#!/usr/bin/python
# -*- coding: UTF-8 -*-

import numpy as np
from MSVD import MSVD
from multirate import multirate 


class QRSdetect:
    
    def __init__(self, ecg, Fs, HR_min = 30):
        
        #taxa de amostragem
        self.Fs = Fs
        #sinal a ser analisadao
        self.ecg = ecg
        #batimentos minimo
        self.HR_min = HR_min
            

    def detect(self, level_dec = 2):


        msvd = MSVD(self.ecg)
        phi, psi, U = msvd.dec(level_dec)
        mt = multirate(psi)
        signal0 = mt.upsample(2**level_dec, -2**level_dec)
        
        msvd = MSVD(self.ecg)
        phi, psi, U = msvd.dec(level_dec + 1)
        mt = multirate(psi)
        signal1 = mt.upsample(2**(level_dec + 1), -2**(level_dec + 1) -2**level_dec)
        
        self.signal = np.abs(signal0) + np.abs(signal1)
                        
        signal_copy = 0
        R_min = (self.HR_min/60.0) * (self.signal.size/self.F       
        self.thr = 0.3*np.max(self.signal)
        
        while True:
            
            start = 0
            end = 0
            Rpeak = []
            
            for x in range(0, self.signal.size):
                
                if self.signal[x] >= self.thr:
                    if start == 0:
                        start = x
                else:
                    if start != 0:
                  
                        if abs(start - x) >= self.Fs * 0.12:
                            end = x - 1
                
                if (start != 0) and (end != 0):
            
                    
                    v_max = max(self.ecg[start:end])
                    v = v_max    
                    
                    t = np.where(self.ecg == v)
                    
                    
                    if t[0].size > 1:
                        idx = t[0]
                        for ii in range(0, t[0].size):
                            if (idx[ii] >= start) and (idx[ii] <= end):
                                t = idx[ii]
                    else:
                        t = t[0]
                    
                    
                    if type(t).__name__ == 'ndarray':
                        t = np.copy(t[0])                    
                    
                    
                    l = len(Rpeak)
                    attr = 0
                    if l > 1:
                        tp, vp = Rpeak[l - 1]
                        if (t - tp) < self.Fs * 0.29:   
                            
                            if v > vp:
                                Rpeak.pop()
                                attr = 1
                        else:
                            attr = 1
                    else:
                        if t > self.Fs * 0.20:    
                            attr = 1
                    
                    if attr == 1:
                        Rpeak.append((t, v))
                        
                    start = 0
                    end = 0
            
            if (len(Rpeak) > R_min):
                break;
            else:
                if (signal_copy == 0):
                    signal_copy = np.copy(self.signal)
                
                lw = int(self.Fs * 0.12)
                wz = np.zeros((1, lw))
                
                l = np.where(self.signal == np.max(self.signal))
                
                signal_copy[abs(int(lw/2) - l[0]):l[0] + int(lw/2)] = wz
                
                self.thr = 0.3*np.max(signal_copy)
           
        self.Rpeak = np.array(Rpeak)
        
        return self.Rpeak
       
    def serie(self):
        self.rr = []
        for i, peak in enumerate(self.Rpeak[:,0]):
            if (i + 1) < self.Rpeak[:,0].size:
                d = self.Rpeak[i + 1, 0] - peak
                self.rr.append(d)
        
        return self.rr
       
  
file multirate.py (link)

Example of detection

#!/usr/bin/python
# -*- coding: UTF-8 -*-
import sys
sys.path.append('/home/.../Scripts/Python')

import numpy as np
import matplotlib
import matplotlib.pylab as pl
from QRSdetect import QRSdetect
import scipy.io
import time

#Configurações do gráficos
matplotlib.rc('xtick', labelsize=8)
matplotlib.rc('ytick', labelsize=8)
   
''' 
ECG
'''
Fs = 500.0
mat = scipy.io.loadmat('/home/.../ECG/C2011/set-a/1002867.mat', mat_dtype=True)
dados = mat['val']

t = np.arange(0, 8.192, 1/Fs)

V1 = dados[6]
V2 = dados[7]
V3 = dados[8]
V4 = dados[9]
V5 = dados[10]
V6 = dados[11]

X_Grupo2 = (np.c_[V1[0:4096], V2[0:4096], V3[0:4096], V4[0:4096], V5[0:4096], V6[0:4096]].T)/100.0

I = dados[0]
II = dados[1]
III = dados[2]
aVR = (-1)*dados[3]
aVF = dados[4]
aVL = dados[5]

X_Grupo1 = (np.c_[I[0:4096], II[0:4096], III[0:4096], aVR[0:4096], aVF[0:4096], aVL[0:4096]].T)/100.0


derivacoes = ['I','II','III','aVR','aVF','aVL','V1','V2','V3','V4','V5','V6']
X = np.c_[X_Grupo1.T, X_Grupo2.T].T

fig1 = pl.figure(1)
fig1.suptitle('Sinais ECG originais', fontsize=14)
fig1.text(0.46,0.03, 'Tempo (s)')
fig1.text(0.07,0.55, 'Amplitude', rotation='vertical')
for i in range(0, 12):
    pl.subplot(4,3,i+1)    
    pl.plot(t, X[i]/np.max(np.abs(X[i])))
    pl.title(derivacoes[i])
    pl.xlim((0, 8.192))
fig0 = pl.figure(2)
fig0.suptitle('Espectrograms', fontsize=14)
fig0.text(0.46,0.03, 'Time (s)')
fig0.text(0.07,0.60, u'Frequency (Hz)', rotation='vertical')
for i in range(0, 12):
    pl.subplot(4,3,i+1)    
    pl.specgram(X[i], Fs = Fs)
    pl.title(derivacoes[i])
    pl.xlim((0, 8))      
       
fig2 = pl.figure(3)
fig2.suptitle(u'R peaks', fontsize=14)
fig2.text(0.46,0.03, 'Sample')
fig2.text(0.07,0.55, 'Amplitude', rotation='vertical')
sigma = 0
media = 0
tempo_medio = 0
for i in range(0, 12):
    pl.subplot(4,3,i+1)

    inicio = time.time()
    
    qrs = QRSdetect(X[i]/np.max(np.abs(X[i])), Fs)
    Rpeak = qrs.detect()   
    fim = time.time()
    serie = qrs.serie()
    print Rpeak[:,0]
    
    tempo_medio += fim - inicio
        
    pl.stem(Rpeak[:,0], Rpeak[:,1])
    pl.title(derivacoes[i])
    pl.xlim((0, 4096))
    pl.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
    med = np.mean(serie)
    media = med + media
    sig = np.std(Rpeak[:,0])/med
    sigma = sig + sigma
    print derivacoes[i]
    print u'Time execution: ', fim - inicio
    print 'Beats: ', Rpeak[:,1].size

figs = [fig1, fig2]

for fig in figs:   
    for ax in fig.get_axes():
        ax.titlesize = 10
        for line in ax.get_lines():
            line.set_color('black')
            
pl.show()



References 

1. Bashir M. E. A., Lee D. G., Aksha M., Yi G. M., Cha E., Bae J., Cho M. C., Ryu K. H, Highlighting the current is- sues with pride suggestions for improving the performance of real time cardiac health monitoring, Information Tech- nology in Bio-and Medical Informatics ITBAM 2010, 6266, 226-233 (2010).
2. Goldberger A. L., Amaral L. A. N., Glass L., Hausdorff J. M., Ivanov P. Ch., Mark R. G., Mietus J. E., Moody G. B., Peng C. K., Stanley H. E., PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals.Circulation101(23): e215-e220 [Circulation Electronic Pages; http://circ.ahajournals.org/cgi/content/full/101/23/e215], (2000).
3. Gussak, I., Antzelevitch, C. and Wilde, A. A. M. and Friedman, P. A. and Ackerman, M. J. and Shen, W.- K., Electrical Diseases of the Heart: Genetics, Mecha- nisms, Treatment, Prevention, 600. Springer-Verlag, Lon- don (2008).
4. Kakarala R., Ogunbona P. O., Signal Analysis Using a Multiresolution Form of the Singular Value Decomposition. IEEE Transaction on Image Processing, 10, 724-735 (2001).
5. Kokkinos, P., Physical Activity and Cardiovascular Dis- ease Prevention, 418. Jones and Bartlett Publishers, On- tario (2010).
6. Lee V. S., Cardiovascular MRI: Physical Principles to Practical Protocols, Lippincott Willians & Wilkins, Philadelphia (1966).
7. Lepage R., Boucher J. M., Blanc J. J., Cornilly, J. C., ECG Segmentation and P-wave Feature Extraction: Application to Patients Prone To Atrial Fibrillation, 23rd Annual Inter- national Conference of the IEEE Engineering in Medicine and Biology Society, 1, 298-301 (2001).
8. Mehta S. S., Lingayat N. S., Comparative Study of QRS detection in Single Lead and 12-Lead ECG Based on En- tropy and Combined Entropy Criteria Using Support Vec- tor Machine, Journal of Theoretical and Applied Informa- tion Technology, 8-18 (2007).
9. Rastogi, N. Mehra, R., Analysis of Savitzky-Golay Filter for Baseline Wander Cancellation in ECG using Wavelets, International Journal of Engineering Sciences & Emerging Technologies, 6, 15-23 (2013).
10. Pan J., Tompkins W. J., A Real-time QRS detection al- gorithm, IEEE Transaction on Biomedical Engineering, 32, 230-236 (1985).
11. Randall D. C., El-Wazir Y. M., ECG Interpretation, Hayes Barton Press, Raleigh (2004).
12. Tanaka H., Monahan K. D., Seals D. R., Age-predicted maximal heart rate revisited. Journal of the American Col- lege of Cardiology, 37, 153-156 (2001).
13. Zidelman Z., Amirou A., Adnane M., Belochrani A., QRS detection base wavelet coefficients. Computer Methods and Programs and Biomedicine, 07, 490-496 (2012).

Links to paper: 

Detecção de Complexos QRS do ECG pela Decomposição em Valores Singulares em Multirresolução - UNESP 

Detecção de Complexos QRS do ECG pela Decomposição em Valores Singulares em Multirresolução - IX ENAMA 

Detecção de Complexos QRS do ECG pela Decomposição em Valores Singulares em Multirresolução - Academia.edu 

QRS detection in 12-lead ECG based Multiresolution Singular Value Decomposition - Academia.edu (english version) 

Detecção de Complexos QRS do ECG pela Decomposição em Valores Singulares em Multirresolução - GoogleScholar

Thursday, December 24, 2015

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.