This article was published in an Elsevier journal. The attached copy is furnished to the author for non-commercial research and education use, including for instruction at the author’s institution, sharing with colleagues and providing to institution administration. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited. In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier’s archiving and manuscript policies are encouraged to visit: http://www.elsevier.com/copyright
Transcript
This article was published in an Elsevier journal. The attached copyis furnished to the author for non-commercial research and
education use, including for instruction at the author’s institution,sharing with colleagues and providing to institution administration.
Other uses, including reproduction and distribution, or selling orlicensing copies, or posting to personal, institutional or third party
websites are prohibited.
In most cases authors are permitted to post their version of thearticle (e.g. in Word or Tex form) to their personal website orinstitutional repository. Authors requiring further information
regarding Elsevier’s archiving and manuscript policies areencouraged to visit:
Wavelet packets approach to blind separation of statisticallydependent sources
Ivica Koprivaa,�, Damir Sersicb
aRudjer Boskovic Institute, Bijenicka cesta 54, P.O. Box 180, 10002 Zagreb, CroatiabDepartment of Electronic Systems and Information Processing, Faculty of Electrical Engineering and Computing, University of Zagreb,
Unska 3, 10000 Zagreb, Croatia
Received 25 September 2006; received in revised form 2 April 2007; accepted 13 April 2007
Communicated by E.W. Lang
Available online 22 April 2007
Abstract
Sub-band decomposition independent component analysis (SDICA) assumes that wide-band source signals can be dependent but some
of their sub-components are independent. Thus, it extends applicability of standard independent component analysis (ICA) through the
relaxation of the independence assumption. In this paper, firstly, we introduce novel wavelet packets (WPs) based approach to SDICA
obtaining adaptive sub-band decomposition of the wideband signals. Secondly, we introduce small cumulant based approximation of the
mutual information (MI) as a criterion for the selection of the sub-band with the least-dependent components. Although MI is estimated
for measured signals only, we have provided a proof that shows that index of the sub-band with least dependent components of the
measured signals will correspond with the index of the sub-band with least dependent components of the sources. Unlike in the case of
the competing methods, we demonstrate consistent performance in terms of accuracy and robustness as well as computational efficiency
of WP SDICA algorithm.
r 2007 Elsevier B.V. All rights reserved.
Keywords: Sub-band decomposition; Independent component analysis; Wavelet packets; Mutual information
1. Introduction
Independent component analysis (ICA) is a statisticaltechnique to extract non-Gaussian and statistically inde-pendent source signals given only the observed ormeasured data [12,23]. The problem is known as blindsource separation (BSS) and is formally described as
xðtÞ ¼ AsðtÞ, (1)
where x 2 RN represents vector of measured signals, A 2RN�M represents an unknown mixing matrix and s 2 RM
represents unknown vector of the source signals. In thesubsequent derivation, we shall assume M ¼ N. For M4N
the BSS problem is reduced on the squared problem by thedimensionality reduction technique, which is realized bythe principal component analysis or singular value decom-position based transforms. The case MoN leads tounderdetermined BSS problem, which is not solvable undergeneral ICA assumptions (non-Gaussianity and statisticalindependence of the source signals). Some additional a
priori information about source signals, such as sparseness[20,29,30,44], must be known in order to solve theunderdetermined BSS problem. Underdetermined case willnot be treated in this paper.Statistical independence assumption of the source signals
is satisfied in many situations which leads to the successfulapplication of the ICA in various fields functional magneticresonance imaging signal processing [33], processing ofradio signals from multiantenna based stations in wirelesscommunication systems [36], multispectral astrono-mical and remotely sensed images [19,34], multiframe
ARTICLE IN PRESS
www.elsevier.com/locate/neucom
0925-2312/$ - see front matter r 2007 Elsevier B.V. All rights reserved.
blind image deconvolution with non-stationary blurringprocess [10,27], speech separation and enhancement[2,28], etc.
However, in practice the independence property of thesource signals does not always hold. The example is brainsource signals observed by EEG [15]. One approach tosolve this problem and relax statistical independenceassumption is to assume that the wideband source signalsare dependent, but there exist some sub-bands where theyare independent. This assumption leads to the SDICA[13–15,38,40,41]. In [14,38], measured signals had to bepassed through the filter bank in order to achieve sub-banddecomposition. The problem remains how to select the sub-band with the least dependent components. Some a priori
knowledge in a form of statistical measure such as kurtosiswas used in [14] while in [38] it has been additionallyassumed the existence of at least two sub-bands where sub-components of the source signals are independent. Theseassumptions may not necessarily hold in practical situa-tions. In addition to that, in the later approach the problemstill remains how to select the best sub-band.
A simplified version of the filter bank approach is to usethe high pass (HP) filter only [14,15] in preprocessing theobserved signals and apply standard ICA algorithm onfiltered data in order to learn the mixing matrix. This ismotivated by the fact that high pass filtered version isusually more independent than original signals. At thesame time, this approach is computationally very efficient,which makes it attractive for BSS problems with statisti-cally dependent sources. Nevertheless, we shall illustratethat due to its simplicity the method is not robust undercertain interference conditions. The similar commentsapply for innovation-based blind separation of statisticallydependent sources [21,24]. The arguments for usinginnovations are that they are usually more independentfrom each other and more non-Gaussian than originalprocesses [21]. In relation to [21] where innovationsrepresentation was proposed as a preprocessing methodto increase accuracy of the linear instantaneous ICAalgorithms, in [24] this was additionally extended to thepost-nonlinear instantaneous ICA problems. In [24] it wasassumed that source signals are statistically independentand temporally correlated. Usage of the linear filter totemporally decorrelate sources will make them more non-Gaussian and consequently increase accuracy of the ICAalgorithms. We comment here that in the context of theBSS problem with statistically dependent sources, we useboth properties of innovations: to be more independentand more non-Gaussain than original processes. Anotherreason to use innovations representation for discussedproblem is its computational efficiency, because it isimplemented by simple linear time invariant filter. How-ever, as in the case of HP filter, we shall illustrate thatinnovations approach is also not robust under certaininterference conditions.
Another approach has been formulated in [40,41] wheremeasured signals were preprocessed by an adaptive filter
with adaptation based on the minimum of mutualinformation (MI) among the filter outputs. In the firstversion of the algorithm [40], prefilter adaptation and de-mixing matrix adaptation of instantaneous ICA stage areimplemented simultaneously. In the second version of thealgorithm [41], the prefilter is adapted alone on the subsetof source signals, which are assumed to be available. ICA isthen applied on filtered observed data. We find theassumption about the availability of the subset ofthe source signals unfair and somewhat unrealistic in thecontext of the blind processing scenario. Additionalcritique of this approach is problems associated with theadaptation of the prefilter. Standard score functions areneeded in learning. It is known that a priori knowledge ofat least statistical class to which source signals belong isrequired for their estimation. In order to obtain optimalscore functions, some form of the kernel-based estimationof the unknown density functions is required [5,35], whichis used in the implementation of the adaptive filteralgorithm [40]. This is one part that adds to thecomputational complexity of this algorithm. Moreover,estimation of the joint score function is also requiredduring filter adaptation phase. The joint score functiondepends on joint density function, which is known to becomputationally very demanding to estimate, especiallywhen the BSS problem is highly dimensional. The tradeoffhas been found in [40] by minimizing pairwise MI.Anyway, this is another part that adds to the computa-tional complexity of this algorithm. We shall demonstratecomputational inefficiency of this algorithm as well asdifficulties associated with the stability of the adaptivelearning process in Section 4.Alternative approach for separation of statistically
dependent sources has been proposed in [7] wheremaximization of the non-Gaussianity measure equivalentto the minimization of the Shannon entropy is used.Performance of this approach strongly depends on thecorrelation level between the sources. In addition to that,there is also relatively high computational cost involvedwith the entropy-based approximation of MI.We propose an approach to SDICA, which is based on
multiresolution decomposition using wavelet packets(WPs) based iterative filter banks [31,39]. In order toenable the filter bank to adaptively select the sub-band withthe least dependent sub-components of the source signals,we have introduced a criterion based on small cumulantapproximation of MI. As it has been demonstrated in theAppendix A the cumulant based approximation is con-sistent estimator of the MI and computationally moreefficient than entropy-based estimator. Although MI isestimated for measured signals only, we provide a proof atthe end of Appendix A that shows that index of the sub-band with least dependent components of the measuredsignals will correspond with the index of the sub-band withleast dependent components of the sources. Thus, ourapproach is adaptive as in [40,41] and computationallyefficient as in [14,15,21,38]. Moreover, it is directly
ARTICLE IN PRESSI. Kopriva, D. Sersic / Neurocomputing 71 (2008) 1642–1655 1643
Author's personal copy
extendable to time-frequency representation. Additionally,as demonstrated in Section 4, it exhibits consistentperformance in terms of robustness and accuracy forvarious characters of interfering signals that act asdependent components.
The rest of the paper is organized as follows. Weintroduce the sub-band decomposition independent com-ponent analysis (SDICA) model in Section 2. The WP-based approach to SDICA is introduced in Section 3.Simulation results with comparative performance evalua-tion are presented in Section 4. Conclusion is given inSection 5.
2. SDICA model
In the BSS problem (1) to be solved by the ICAalgorithms, it is assumed that sources sm(t), mA{1,y,M},are non-Gaussian and statistically independent. A powerfulextension and generalization of this basic ICA model isSDICA. It assumes that wide-band source signals can bedependent, but some of their sub-components are inde-pendent. Thus, source signals can be represented as[14,15,38,40,41]:
smðtÞ ¼ sm;1ðtÞ þ sm;2ðtÞ þ . . .þ sm;LðtÞ, (2)
where sm,k(t), k ¼ 1,y,L, are narrow-band sub-compo-nents. The ICA problem is to find a separation or demixingmatrix W9A�1, which recovers independent componentsor source signals
yðtÞ ¼WxðtÞ, (3)
where y 2 RM and y9s. Like [40,41] we shall assume thatfor certain set of k, sub-components in (2) are leastdependent or possibly independent. This significantlyrelaxes assumptions on the sub-band model (2) in relationto assumptions made in [14], where some a priori
information about the sources is assumed. That contradictsto blindness assumption used to solve the BSS problem (1),or in [38], where existence of at least two sets of k sub-components is assumed. Under presented assumption thestandard linear ICA algorithms can be applied to theselected set of k sub-components in order to learn de-mixing matrix W:
ykðtÞ ¼WxkðtÞ. (4)
The problems to be resolved are which preprocessingtransform should be used to obtain sub-band representa-tion of the original wideband BSS problem (1) and whichcriteria should be used to select the set k with leastdependent sub-components.
3. Multiresolution SDICA
In order to obtain sub-band representation of theoriginal wideband BSS problem (1), we can use any linearoperator Tk which will extract a set k of sub-components
skðtÞ ¼ Tk sðtÞ½ �, (5)
where Tk can, for example, represent a linear time-invariant bandpass filter. Using (5) and sub-band repre-sentation of the sources (2), application of the operator Tk
on the wideband BSS model (1) yields
xkðtÞ ¼ Tk AsðtÞ½ � ¼ ATk sðtÞ½ � ¼ AskðtÞ. (6)
For this purpose a fixed filter bank has been used in[14,38]. A high pass filtering can be seen as a special case ofthe filter bank approach. An adaptive pre-filter has beenapplied on xn, nA{1,y,N} in [40,41] and trained byminimizing MI information among the filter outputs, thuseliminating dependent sub-components sk. In this paper wepropose WP transform to be used for Tk in order to obtainsub-band representation of the wideband BSS problem (1).The main reason is existence of the WP transform in a formof iterative filter bank and multiresolution property of theWP transform which allows isolation of the fine detailswithin each decomposition level and enable adaptive sub-band decomposition [31,39]. In the context of SDICA, itmeans that an independent sub-band that is arbitrarilynarrow can be isolated, by progressing to the higherdecomposition levels. Also, computationally efficient im-plementations of the WP transform exist for both two-dimension (2D) and 1D signals. This will be illustrated inthe next section by applying WP-based SDICA approachto blind separation of human faces, as well as to 1D signalswith a sine wave as interferer, i.e. dependent component.We further use properties of the WP transform in order toconfirm existence of (6). That property was extensivelyexploited in various versions of the sparse ICA, where ithas been found that either WP or short-time Fouriertransform are very useful for obtaining new representationof data which is sparser than original formulation. As ithas been shown, executing ICA in sparse domain producedmore accurate solutions in solving linear instantaneousBSS problem and enabled to solve underdetermined (moresources than sensors) BSS problem [20,25,29,30,44]. Inparticular case of WP, we express each source signal(image) in terms of its decomposition coefficients:
sjkmðtÞ ¼
Xl
cjkmljjlðtÞ, (7)
where j represents scale level, k represents sub-band index,m represents source index and l represents shift index. jj(t)is the chosen wavelet also called atom or element of therepresentation space and c
jkml are decomposition coeffi-
cients. In our implementation of the described SDICAalgorithm, we have used shift-invariant 2D WP decom-position for separation of human faces and 1D WP forseparation of 1D signals. If we choose the same representa-tion space as for the source signals, we express eachcomponent of the observed data x as
xjknðtÞ ¼
Xl
fjknljjlðtÞ, (8)
where n represents sensor index. Let vectors fl and cl beconstructed from the lth coefficients of the mixtures and
ARTICLE IN PRESSI. Kopriva, D. Sersic / Neurocomputing 71 (2008) 1642–16551644
Author's personal copy
sources, respectively. From (1) and (8) using the orthogon-ality property of the functions jj(t) we obtain
f l ¼ Acl . (9)
If additive noise is present this relation holds approxi-mately, i.e. expanding the noise term in the samerepresentation n
jkmðtÞ ¼
Ple
jkmljjlðtÞ would add up to (9)
a vector el of the decomposition coefficients of the noise[25]. Thus, when the noise is present Eq. (9) becomesflEAcl. From (1) and (9) we see the same relation betweensignals in the original domain and WP representationdomain. Inserting (9) into (8) and using (7) we obtain:
xjkðtÞ ¼ As
jkðtÞ, (10)
as it was announced by (6) where no assumption ofmultiscale decomposition has been made. For eachcomponent xn of the observed data vector x, the WPtransform creates a tree with the nodes that correspond tothe sub-bands of the appropriate scale. In order to selectthe sub-band with least dependent components sk, wemeasure MI between the same nodes in the WP trees. Forthis purpose we use the small cumulant approximation ofthe Kullback–Leibler divergence, which is an exactmeasure of MI, obtained under weak correlation and weaknon-Gaussianity assumptions [9]:
Ij
k xjk1;x
jk2; . . . ;x
jkN
� ��
1
4
X0pnolpN
nal
cum2ðxjkn;x
jklÞ
þ1
12
X0pnolpN
nal
cum2ðxjkn; x
jkn;x
jklÞ þ cum2ðx
jkn;x
jkl ; x
jklÞ
� �
þ1
48
X0pnolpL
nal
cum2ðxjkn;x
jkn; x
jkn;x
jklÞ
�
þcum2ðxjkn;x
jkn;x
jkl ; x
jklÞ þ cum2ðx
jkn;x
jkl ;x
jkl ; x
jklÞ
�, ð11Þ
where cum() in (11) denotes second, third or fourth-ordercross-cumulants [6,32]. Approximation of the joint MI bythe sum of pair-wise MI is commonly used in the ICAcommunity in order to simplify computational complexityof the linear instantaneous ICA algorithms [17]. We havedemonstrated in Appendix A as follows: (i) that cumulantbased approximation (11) of the MI is consistent as andcomputationally more efficient than entropy based approx-imation of MI [18]; (ii) although MI is estimated formeasured signals only index of the sub-band with leastdependent components of the measured signals willcorrespond with the index of the sub-band with leastdependent components of the sources. However, becausewe use small cumulant-based approximation (11) of theMI, it means that sometimes instead of selecting a sub-band with least dependent sub-components of the sourcesignals we shall only get close to this sub-band. Byconsistency, we assume that for mutually independentprocesses approximation (11) of the MI converges towardzero when the sample size grows toward infinity and that
its value is increased when dependence level between theprocesses is increased. Once the sub-band with the leastdependent components is selected, we obtain eitherestimation of the inverse of the basis matrix W orestimation of the basis matrix A by applying standardICA algorithms on the model (10). Reconstructed sourcesignals y are obtained by applying W on the original data x
as it is given by (3). Alternatively, mixed signals can bereconstructed through the synthesis part of the WPtransform, where sub-bands with a high level of MI areremoved from the reconstruction. This is demonstrated inexperiments 1 and 2 in the next section. We summarizemultiscale SDICA BSS algorithm in the following foursteps:
1. Perform multiscale WP decomposition of each compo-nent of the multivariate data x. Wavelet tree will beassociated to each component of x (Eq. (6)–(10).
2. Select sub-band with the least dependent components byestimating MI between the same nodes (sub-bands) inthe wavelet trees (Eq. (10) and (11).
3. Learn basis matrix A or its inverse W by executingstandard ICA algorithm for linear instantaneous pro-blem on the selected sub-band (Eq. (10)).
4. Obtain recovered sources y by applying W on datavector x (Eq. (3)).
4. Simulation examples
In this section we examine performance of derived WPSDICA algorithm using the same examples as in [40]. Thecode for MATLAB implementation of the WP SDICAalgorithm can be found and downloaded from [26]. Inorder to quantify separation quality, we also use theAmari’s performance index Perr [1] that measures closenessof the global separation matrix Q ¼WA to the generalpermutation matrix P ¼ KI, where K represents diagonalmatrix and I represents identity matrix. Perr is calculated as
Perr ¼1
NðN � 1Þ
XN
i¼1
XN
j¼1
qij
�� ��maxk qik
�� ��� 1
!(
þXN
j¼1
qji
�� ��maxk qki
�� ��� 1
!, ð12Þ
where qij ¼ [Q]ij and 0pPerrp2. Perr approaches zero whenQ approaches P. We compare WP SDICA algorithm withadaptive prefiltering algorithms [40,41], innovations ap-proach [21,24] and high pass prefiltering approach [14,15].We have implemented the high pass prefiltering by meansof 1D wavelet transform with the one decomposition level.It realizes coarse approximation and details by means ofthe two half band filters. All the experiments have beenconducted in MATLABr environment on a 3GHz PCmachine with dual core microprocessor and 2GB ofinternal memory. In a case of all tested algorithms, wehave used the FastICA algorithm in a symmetric mode
ARTICLE IN PRESSI. Kopriva, D. Sersic / Neurocomputing 71 (2008) 1642–1655 1645
Author's personal copy
with tanh nonlinearity [22]. We would like to comment thatchoice of the nonlinearity for the FastICA algorithm is notcritical as it would be for the ICA algorithm based on theclassical score function [4,43]. We have also tested FastICAalgorithm with cubic nonlinearity as well as JADEalgorithm [8] that is known to be ‘‘nonlinearity free’’.The obtained results were only slightly different than thoseobtained by FastICA algorithm with tanh nonlinearity.The reason why FastICA and any standard ICA algorithmwill fail when applied directly on measured data is violationof statistical independence assumption.
4.1. Experiment 1: WP-based SDICA with artificially
generated data
In this experiment, we have used the four independentsub-components: amplitude-modulated signal, a squarewave, a high-frequency noise signal and a speech signal.Each signal has 10,000 samples. Each original signal si
contains one of the above independent signals togetherwith a sine wave with the same frequency, O ¼ 0.41 rad,but different phases for different sources, which representsdependent sub-component. Fig. 1 shows these signals aswell as their magnitude spectra. Fig. 2 shows waveforms(left) and corresponding power spectrums (right) of theobserved or mixed signals xi. We have used the samemixing matrix as in [40]. Fig. 3 shows waveforms of thetrue dependent source signals (left) and waveforms of theestimated dependent source signals (right) reconstructed bydirect application of the FastICA algorithm on theobserved signals xi. It is evident that separation quality is
poor. The value of the corresponding Amari’s error isPerr ¼ 0.4728, which is significantly greater than zero.Evidently, the presence of the dependent sub-componentsprevented ICA algorithm from learning de-mixing matrixW correctly. Fig. 4 shows equivalent results when describedWP-based SDICA has been applied on observed data. Wehave used 1D shift-invariant WPs with five decompositionlevels. Regarding the type of the wavelet our choice wassymmlets with eight vanishing moments. The amount ofthe MI equals to 9.3059e-8, where MI has been normalizedwith respect to the maximal value in the whole WP tree. Inorder to derive a truly unsupervised approach to blindseparation of statistically dependent sources, we need acriterion for the selection of decomposition level. One wayto automate algorithm completely is to monitor the rate ofchange of the MI, as decomposition level is increasing andstop when rate of change is small. For example, bycarefully looking Table 1 we observed that change of MI inthe part of the wavelet tree that spreads from the high passpart of the decomposition level 1 is insignificant incomparison with the MI at the high pass part ofdecomposition level 1. This issue is commented in moredetails at the end of the section. Table 1 shows normalizedMI among the corresponding nodes of the wavelet treesobtained after WP-based decomposition of the observeddata. The observed data were reconstructed using thesynthesis part of the WP tree. Only the nodes withnormalized MIo0.05 have been retained, i.e. the nodeswith higher level of MI were eliminated. In this way thedependent sub-component has been treated as interfererand eliminated from measurements before separation. We
ARTICLE IN PRESS
Fig. 1. The source independent sub-components (top four) and the waveform of dependent sub-component (the bottom one), as well as their magnitude
spectra. Only 400 samples are plotted for illustration. The sources are si ¼ si,I+si,D and si,D are sinusoid waves with the same frequency but different phases.
I. Kopriva, D. Sersic / Neurocomputing 71 (2008) 1642–16551646
Author's personal copy
then applied FastICA algorithm on the reconstructedobserved data in order to learn de-mixing matrix W. TheAmari’s error was Perr ¼ 0.0051, that is 92.7 times less thanthe value when the FastICA algorithm has been applieddirectly on measured data. The whole process, analysis andsynthesis part, took approximately 9.5 s.
The competing innovations approach has been testedwith a 10th order of the AR model. The prediction filter
was very efficient in eliminating the dependent sine wavecomponent. The FastICA algorithm has been applied oninnovations to learn the de-mixing matrix W. The Amari’serror was Perr ¼ 0.0885 that is 5.3 times less than the valuewhen the FastICA algorithm has been applied directly onmeasured data. The whole process took approximately0.3 s. When HP filter has been applied to solve describedBSS problem value of the Amari’s error was Perr ¼ 0.0823,
ARTICLE IN PRESS
Fig. 2. The observed signals (left) and their magnitude spectra (right) in the experiment 1.
Fig. 3. The source independent sub-components with added dependent component (left) and reconstructed waveforms by direct application of FastICA
method on the observed data (right).
I. Kopriva, D. Sersic / Neurocomputing 71 (2008) 1642–1655 1647
Author's personal copy
that is 5.7 times less than the value when the FastICAalgorithm has been applied directly on measured data. Thewhole process took approximately 0.22 s. The HP filterused in the experiment has been obtained from the WPs. Itis interesting to note that use of simple differentiator, i.e.FIR filter with coefficients [1 �1], that works very well withface images gives very poor performance in the presentedcase of 1D signals. For such a filter value of Amari’s errorwas Perr ¼ 0.5646.
Thus, it is evident that both HP and innovations arereasonably accurate and yet computationally very efficient.On the other hand, WP approach offers an order ofmagnitude better accuracy with an increased computa-tional time that is still of the order of a few seconds.However, if we change the frequency of the sine dependentcomponents to O ¼ 2.14 rad, the WP approach achievesthe value of Amari’s error Perr ¼ 0.0069, the innovationsapproach Perr ¼ 0.01271 and the HP filtering approachPerr ¼ 0.4423. The reason while HP filtering approachfailed has been of course that interference is in a high passpart of the spectrum. Without a priori information aboutlocation of the interference HP approach fails. By thissimple experiment we have demonstrated the meaning andimportance of adaptability of the proposed WP SDICAapproach.
The adaptive prefilter algorithm has been applied onmixed data with adaptive filter learning gain Zh ¼ 0.1, de-mixing matrix learning gain ZW ¼ 0.15 and filter orderequal to 17. These values were suggested in [40]. We haveused the code provided on the web-site [42]. The algorithmfailed to remove sine dependent component. When applied
in independence enhanced mode [40], which assumesavailability of the subset of source signals, the adaptivefilter strongly suppressed sine dependent component andachieved value of the Amari’s error was Perr ¼ 0.0046.However, we find assumption about availability of thesubset of source signals unfair in the context of blindscenario.In [16] additive information cost function was proposed
for the selection of the best basis of WPs. Table 1 shows thevalues of normalized MI for all nodes in the WP tree. It isimportant to note that MI is not an additive function thatfits in the best basis scheme. Nevertheless, we can compareMI measure in the parent and child nodes and decidewhether the further splitting was fruitful. If the child nodesdo not differ significantly among each other in the MI, wecan stop the further splitting in that direction.
4.2. Experiment 2: WP-based SDICA for separating images
of human faces
This example has been chosen from the same reason as in[40]. Namely, it has been already found out that faceimages, represented as 1D signals, are dependent and hardto separate [21]. We demonstrate here the capability of theWP SDICA algorithm to successfully separate four imagesof human faces and compare its performance with adaptiveprefiltering algorithm [40], innovations approach and HPprefiltering approach. We have used differentiator with theimpulse response h ¼ [1 �1] to implement HP filter forwhich it is known to works well with the face images.Regarding the type of the wavelet, we have also used
ARTICLE IN PRESS
Fig. 4. Reconstructed source independent sub-components (left) and their magnitude spectra (right) obtained by WP-based SDICA. It can be seen from
the magnitude spectra that dependent component with OE0.41 rad is removed from the reconstructed sources.
I. Kopriva, D. Sersic / Neurocomputing 71 (2008) 1642–16551648
Author's personal copy
symmlets with eight vanishing moments. The four originalimages are shown in Fig. 5. It has been already demon-strated that both innovations and HP prefiltering aresuccessful in separating faces from their mixtures [21,40].However, let us assume the existence of the externalbackground Gaussian noise that is treated in BSS scenarioas another source signal [12]. It obviously represents a
dependent sub-component. We have added it such that anaverage SNR with respect to source images was 29.1 dB.The noise corrupted source images are mixed with the 4� 4random mixing matrix. Mixed images are shown inFig. 6(a). It can be observed that external backgroundnoise is not visible. Fig. 6(b) shows separated imagesobtained after applications of the FastICA algorithm onthe mixed images. The quality of the separated images isevidently very poor, that is also confirmed by the value ofthe Amari’s error Perr ¼ 0.7382. If we now apply the WPSDICA algorithm, innovations approach and HP prefilter-ing approach to these noisy data, we respectively obtain thefollowing values of the Amari’s error Perr ¼ 0.0344, 0.2761and 0.1685. Performance of WP SDICA algorithm hasbeen improved 21.4 times, performance of the innovationsapproach has been improved only 2.7 times and perfor-mance of the HP prefiltering approach has been improved4.38 times. Fig. 6(c) shows separated image by WP SDICAalgorithm, while Fig. 6(d) shows separated images byinnovations approach. Normalized mutual informationamong the corresponding nodes of the wavelet treesobtained after WP based decomposition of the observeddata is shown in Table 2. The observed data werereconstructed using the synthesis part of the wavelet packettree. Only the nodes with normalized mutual informationless than 0.05 have been retained. The reason whyinnovations failed is that dependent white Gaussian noiseterm is not predictable. The HP prefiletring approach faileddue to the wideband character of the external noise term,i.e. the HP filtering kept half of its power. The WP SDICAalgorithm with two decomposition levels took approxi-mately 91 s. The innovations approach took approximately19 s while HP prefiltering approach took approximately12 s. When WP transform was implemented in decimatedversion obtained Amari’s error for two decompositionlevels was Perr ¼ 0.04, which is slightly worse than for non-decimated version while the execution time was equal to14 s. This experiment is additional evidence about therobustness of the WP SDICA algorithm as well as itsflexibility in achieving tradeoff between accuracy andcomputational efficiency.The adaptive prefilter algorithm has been additionally
tested on a case without external noise with the parametervalues as proposed in [40] and using code provided byZhang et al. [42]. The filter order was 11, filter learning
ARTICLE IN PRESS
Fig. 5. Separating mixtures of images of human faces: original images.
Table 1
Values of the normalized mutual information between the same nodes in
the wavelet packet trees of the mixed signals in experiment 1
Node that corresponds to the minimum MI is marked with an asterisk,
while nodes that represent sub-bands with a high level of mutual
information are marked with two asterisks. The italic nodes denote the
sub-bands from which the partial reconstruction has been done.
I. Kopriva, D. Sersic / Neurocomputing 71 (2008) 1642–1655 1649
Author's personal copy
gain 0.04 and de-mixing matrix learning gain in theinstantaneous ICA stage 0.08. Initial value of de-mixingmatrix has been obtained by the FastICA algorithm. Thebest result obtained after 265 iterations was Perr ¼ 0.0134.However after 500 iterations algorithm stopped at similarvalue Perr ¼ 0.0124, but the error has been oscillatingduring the adaptation process. The value of the Amari’serror vs. iteration index is shown in Fig. 7(a), whichillustrates oscillating character of the learning process. Wehave decreased the prefilter and de-mixing matrix learninggains 10 times to 0.004 and 0.008. Fig. 7(b) shows behaviorof the Amari’s error vs. iteration index for 5000 iterations.Minimum of the Amari’s error Perr ¼ 0.0077 has been
achieved after 4567 iterations. Evidently, the learning ismore stable but small oscillations still remained. We foundvery difficult to determine when to stop the algorithm. Theoscillating character of the learning process may representserious difficulty in real world applications where Amari’serror is impossible to be calculated and monitored. Thereare no such difficulties associated with the WP SDICAapproach. We have also tested the [41] on the sameexample in a slightly modified mode, i.e. the prefilter hasbeen applied on observed data only and ICA algorithm hasbeen applied to filtered data in order to learn de-mixingmatrix. Obtained results are equivalent to those alreadyreported when filtering and instantaneous ICA are
ARTICLE IN PRESS
Fig. 6. Separating mixtures of images of human faces with an external Gaussian noise source. Average SNR value was equal to 29.1 dB. (a) Mixed images.
(b) Separation results obtained by FastICA algorithm. (c) Separation results obtained by WP SDICA algorithm. (d) Separation results obtained by
innovations based algorithm.
I. Kopriva, D. Sersic / Neurocomputing 71 (2008) 1642–16551650
Author's personal copy
executed together [40]. Regarding computational complex-ity of the [40] algorithm it took approximately 360 s for 500iterations. For an accuracy of Perr ¼ 0.0077 achieved after4567 iterations it took 3240 s. The same accuracy has beenachieved by WP SDICA algorithm in 91 s only.
5. Conclusion
We have formulated a new method for blind separationof statistically dependent sources. The method is basedupon assumption that wideband sources are dependent,but there exists sub-band where sources are independent.Adaptive sub-band decomposition is realized through WP
transform implemented in a form of iterative filter bank.The sub-band with least dependent components is selectedby measuring MI between the same nodes of the wavelettrees in the multiscale decomposition of the mixed signals.The computationally efficient small cumulant approxima-tion of the MI has been used for this purpose. Ifreconstruction of dependent sub-components is desired,the de-mixing matrix, learnt by applying the standard ICAalgorithm on selected sub-band mixed signals, is applied onoriginal mixed signals. This is illustrated by successfulseparation of images of human faces from the mixtures. Ifreconstruction of dependent sub-components is not de-sired, the learnt de-mixing matrix is applied on a cleaned
version of the mixed signals obtained through the synthesispart of the WP transform, where only sub-bands withnormalized MI less than predefined threshold retained.Demixing matrix is learnt by applying ICA on recon-structed mixed signals with eliminated dependent sub-components. This is illustrated by successful separation ofartificially generated 1D signals corrupted by a sine wave asdependent sub-component. In relation to the innovationsalgorithm proposed in [21] and high pass prefilteringalgorithm proposed in [14,15], it has been demonstratedthat WP SDICA algorithm exhibits consistent performancein terms of accuracy and robustness with respect to variousinterfering signals that act as dependent components, whilethe other two methods fail under certain circumstances.Yet, the WP SDICA algorithm keeps reasonably smallcomputational complexity. Unlike the adaptive prefilteringalgorithms proposed in [40,41], the WP SDICA algorithmdoes not have any stability problems associated with anadaptive learning process and is an order of magnitudecomputationally more efficient. Thus, we consider that WPSDICA algorithm may be quite useful in practicalapplications due to its robustness and computationalefficiency.
ARTICLE IN PRESS
Table 2
Values of the normalized mutual information between corresponding
nodes in the wavelet packet trees of the mixed signals in experiment 2
Level 0 Level 1 Level 2
3.1104e-2 3.6502e-2 1.00000**
3.4570e-2
2.5430e-2*
3.2171e-2
3.9920e-2 5.2337e-2**
3.8559e-2
2.9874e-2
3.3041e-2
2.8660e-2 2.9893e-2
2.8736e-2
2.7758e-2
3.1539e-2
3.0871e-2 2.8029e-2
2.9845e-2
3.0672e-2
3.1825e-2
For description of used notation see Table 1.
Fig. 7. Amari’s error vs. iteration index for algorithm (Zhang and Chan, 2006a) in blind separation of human faces. (a) Prefilter learning gain 0.04, de-
mixing matrix learning gain 0.08. (b) Prefilter learning gain 0.004, de-mixing matrix learning gain 0.008.
I. Kopriva, D. Sersic / Neurocomputing 71 (2008) 1642–1655 1651
Author's personal copy
Acknowledgment
The authors wish to thank the anonymous reviewers fortheir constructive comments that certainly improved thequality of the paper. Part of this work was supportedthrough grant 098-0982903-2558 funded by the Ministry ofScience, Education and Sport, Republic of Croatia.
Appendix A. Approximation of the MI and its use in sub-
band selection
Under weak correlation and weak non-Gaussian as-sumptions, it has been shown in [9], Eqs. (70), (66) and(68), that MI can be approximated via small cumulantapproximation of the Kullback–Leibler divergence(Gram–Charlier expansion of non-Gaussian distributionsaround normal distribution) as
I y1; y2; . . . ; yN
� ��
1
4
X1piojpN
iaj
c2ij þ1
2
XrX3
1
r!
XN
i1i2;...;ir¼1
c2i1i2;...;ir,
(A.1)
where c2i1i2;...;irdenotes square of the related rth order cross-
cumulant and i1i2,y,ir denotes partition of indices suchthat they are not all identical and cij denotes cross-correlation or second-order cross-cumulant. If we furtherassume that all the cumulants of the order higher than 4(weak non-Gaussian assumptions) are small, then (A.1) isreduced to
I y1; y2; . . . ; yN
� ��
1
4
X1piojpN
iaj
c2ij þ1
12
XN
i1i2i3¼1
c2i1i2i3
þ1
48
XN
i1i2i3i4¼1
c2i1i2i3i4. ðA:2Þ
We now apply another approximation commonly donein ICA community [17] and approximate joint indepen-dence I(y1,y2,y,yn) by the sum of pairwise independencesand arrive to
I y1; y2; . . . ; yN
� ��
X1piojpN
iaj
I yi; yj
� ��
1
4
X1piojpN
iaj
c2ij
þ1
12
X1piojpN
iaj
Xn;m
cumðyi; . . . ; yi|fflfflfflfflffl{zfflfflfflfflffl}n times
; yj ; . . . ; yj|fflfflfflfflffl{zfflfflfflfflffl}m times
Þ
0B@
1CA
2
þ1
48
X1piojpN
iaj
Xp;q
cumðyi; . . . ; yi|fflfflfflfflffl{zfflfflfflfflffl}p times
; yj ; . . . ; yj|fflfflfflfflffl{zfflfflfflfflffl}q times
Þ
0B@
1CA
2
, ðA:3Þ
where in (A.3) n,mA{1,2}, p,qA{1,2,3}, n+m ¼ 3 andp+q ¼ 4. As it has been proven in [9] it applies for sucha measure I(y1,y2,y,yN)X0. It is equal to zero when
yn
� N
n¼1are mutually statistically independent or pairwise
statistically independent. This pairwise approximation ofMI has been also used in [40], Eq. (19), as a cost functionfor filter adaptation. It follows that above approximationrepresents consistent measure of statistical (in)dependence.From the practical reasons we shall rewrite Eq. (A.3) interms of explicit use of second-, third- and fourth-ordercross-cumulants
I c y1; y2; . . . ; yN
� ��
1
4
X1piojpN
iaj
cum2ðyi; yjÞ
þ1
12
X1piojpN
iaj
cum2ðyi; yi; yjÞ þ cum2ðyi; yj ; yjÞ
� �
þ1
48
X1piojpN
iaj
cum2ðyi; yi; yi; yjÞ þ cum2ðyi; yi; yj ; yjÞ
�
þcum2ðyi; yj ; yj ; yjÞ
�, ðA:4Þ
where subscript c denotes cumulant-based approximation.Alternative to this approximation is to employ the entropy-based approximation that follows from the definition of theMI [18]:
IðyÞ ¼XN
n¼1
HðynÞ �HðyÞ, (A.5)
where H denotes the entropy, HðynÞ ¼ �E½logðpynðynÞÞ�
and pynis the density function of yn. It is necessary to
approximate marginal entropies, joint entropy and relateddensity functions:
IH ðyÞ ¼XN
n¼1
HðynÞ � HðyÞ, (A.6)
where subscript H denotes entropy-based approximation.We follow approach exposed in [3]:
HðynÞ ¼ �1
T
XT
t¼1
log pynðynðtÞÞ,
pynðynðtÞÞ ¼
1
T
XT
m¼1
Kh ynðtÞ � ynðmÞ� �
,
HðyÞ ¼ �1
T
XT
t¼1
log py yðtÞð Þ,
pyðyðtÞÞ ¼1
T
XT
m¼1
YNn¼1
Kh ynðtÞ � ymðtÞ� �
, (A.7)
where T denotes the sample size. In (A.7) Kh representskernel. Various kernels might be used for density estima-tion [11,35,37]. The commonly used kernel is Gaussiankernel
KhðynðtÞ � ynðmÞÞ ¼ exp �ynðtÞ � ynðmÞ� �2
h
!, (A.8)
ARTICLE IN PRESSI. Kopriva, D. Sersic / Neurocomputing 71 (2008) 1642–16551652
Author's personal copy
where h denotes the bin width of the kernel. Recommenda-tion for choosing h is [37]:
h ¼ 0:6sffiffiffiffiT5p , (A.9)
where s denotes sample standard deviation of yn. The jointdensity estimator in (A.7) is expressed using a product of1D kernels. This simplifies computational complexity, butrepresents a restriction since joint density does not factorizeexcept if marginal variables are independent. However, ithas been proven in [3] that this choice improves the bias ofthe estimator of the MI that asymptotically convergestoward zero with a speed 1/T. Therefore, MI estimator(A.6) satisfies IH ðyÞX0 with the equality to zero only whenvariables of y are statistically independent.
We now demonstrate consistency of the two MIestimators I cðyÞ and IH ðyÞ. Fig. A.1 shows MI as afunction of the sample size for two Laplacian distributedprocesses, where ‘x’ denotes entropy-based approximationand ‘o’ denotes cumulant-based approximation. It isevident that both approximations converge toward zeroas sample size increases. Fig. A.2 shows approximations ofthe MI for two partially dependent processes s1 ¼ s1
u+csd
and s2 ¼ s2u+csd where s1
u and s2u were two zero mean
uniformly distributed processes and dependent componentsd being normally distributed process with zero mean andunit variance. The scale factor c has been varied between0.1 and 1 with the step equal to 0.1. It is evident fromFig. A.2 that both I c and IH are consistent, i.e.I cðHÞ s1ðc1Þ; s2ðc1Þð Þ4I cðHÞ s1ðc2Þ; s2ðc2Þð Þ for c14c2. In thecase of dependent processes both approximations are lessaccurate, IH due to the already discussed reasons related tothe use of 1D kernel in the joint density estimator (A.7) andI c due to the fact that only cumulants up to the order fourwere used in the approximation (A.4). Nevertheless,demonstrated consistency makes both measures suitablefor detection of the sub-band with the least dependentcomponents. Due to the huge difference in computationalcomplexity, we have selected cumulant-based approxima-tion I c for the use in sub-band detection criterion. For thesake of illustration, Fig. A.3 shows estimated computationtime for MATLAB implementation of both approxima-
tions as a function of the sample size for the exampleanalyzed in Fig. A.1. Note that small inconsistency withthe cumulant-based measure is due to the poor resolutioncapability of the MATLAB functions clock and etime, onthe milliseconds scale, that were used for measuring thecomputation time.After demonstrating consistency of the MI estimators
I cðyÞ and IH ðyÞ, we now want to prove that index of thesub-band with the statistically least dependent componentsof the measured signals, that is found at the minimum ofthe MI of the measured signals, corresponds with the indexof the sub-band with the statistically least dependentcomponents of the source signals. For that purpose, werewrite sub-band decomposition representation, Eq. (6), ofthe wideband BSS model:
xkðtÞ ¼ AskðtÞ. (A.10)
Index of the sub-band with the statistically leastdependent components of the measured signals k* isobtained as
k� ¼ argmink
IðxkÞ. (A.11)
ARTICLE IN PRESS
Fig. A1. MI approximation as a function of the sample size for two Laplacian distributed processes. ‘x’ denotes entropy-based approximation while ‘o’
denotes cumulant-based approximation.
Fig. A2. MI approximation as a function of dependence factor. Two
uniformly distributed processes were used as independent processes.
Normally distributed process was added as a dependent process with the
scale factor that influenced dependence level. ‘x’ denotes entropy-based
approximation while ‘o’ denotes cumulant based approximation.
I. Kopriva, D. Sersic / Neurocomputing 71 (2008) 1642–1655 1653
Author's personal copy
Consequently, we want to prove
k� ¼ argmink
IðxkÞ ¼ argmink
IðskÞ. (A.12)
Let us now express the MI in a form of Kullbackdivergence:
IðxkÞ ¼ D pðxkÞYNn¼1
pnðxknÞ
����� !
¼
ZpðxkÞ log
pðxkÞQNn¼1pnðxknÞ
¼XN
n¼1
HðxknÞ �HðxkÞ, ðA:13Þ
where H(xkn) and H(xk) respectively represent marginaland joint differential entropy [18, pp. 224–233]. MI definedas Kullback divergence is a convex function with theproperty I(xk)X0 and equality with zero only if xkn areindependent, i.e. pðxkÞ ¼
QNn¼1pnðxknÞ, see theorem 9.6.1 in
[18]. Also from the corollary defined by Eq. (9.59) in [18], itfollows:
HðxkÞpXN
n¼1
HðxknÞ, (A.14)
with equality only if xkn are independent. Because we arelooking for the sub-band with least dependent componentsof the measured signals, sum of the marginal entropies insuch a case will be closest to the joint entropy. Therefore,from (A.13) and (A.14) it follows:
k� ¼ argmink
IðxkÞ ¼ argmink
XN
n¼1
HðxknÞ � argmaxk
HðxkÞ.
(A.15)
It follows from (A.14) and (A.15) that minimizing theMI is equivalent to the simultaneous minimization of thesum of marginal entropies and maximization of the jointentropy. Thus, maximization of joint entropy impliessimultaneous minimization of the sum of marginalentropies driving MI toward zero. We now express jointentropy H(xk) in terms of the joint entropy H(sk) using therule for the entropy of the linear transformation xk ¼ Ask
[18] (theorem 9.6.4, i.e. Eq. (9.67)):
HðxkÞ ¼ HðskÞ þ log jdet Aj. (A.16)
Taking account that second term in (A.16) is invariantwith respect to the sub-band index k it follows:
argmaxk
HðxkÞ ¼ argmaxk
HðskÞ. (A.17)
Due to already discussed reasons maximization of H(sk)will simultaneously imply minimization of
PNn¼1HðsknÞ.
Hence
k� ¼ argmink
IðxkÞ
¼ argmink
XN
n¼1
HðxknÞ � argmaxk
HðxkÞ
¼ argmink
XN
n¼1
HðsknÞ � argmaxk
HðskÞ
¼ argmink
IðskÞ, ðA:18Þ
which proves (A.12). However, because we use smallcumulant-based approximation (A.4) of the MI (A.5), theequality (A.12) holds approximately. It means that some-times instead of selecting a sub-band with least dependentsub-components of the source signals we shall only getclose to it.
References
[1] S. Amari, S.A. Cichocki, A.H.H. Young, A new learning algorithm
for blind signal separation, MIT Press, Cambridge, MA, 1996, pp.
757–763.
[2] S. Araki, S. Makino, R. Aichner, T. Nishikawa, H. Saruwatari,
Subband-based blind separation for convolutive mixtures of speech,
IEICE Trans Fundament E88-A (2005) 3593–3603.
[3] S. Archard, D.T. Pham, Ch. Jutten, Criteria based on mutual
information minimization for blind source separation in post
nonlinear mixtures, Signal Process 85 (2005) 965–975.
[4] A.J. Bell, T.J. Sejnowski, An information-maximization approach to
blind separation and blind deconvolution, Neural Comput. 7 (1995)
1129–1159.
[5] R. Boscolo, H. Pan, Independent component analysis based on
nonparametric density estimation, IEEE Trans. Neural Networks 15
(2004) 55–65.
ARTICLE IN PRESS
Fig. A3. Computation time for two approximations of the MI as a function of the sample size. Two Laplacian distributed processes were used as
