# Extraction of Correlated Sparse Sources from Signal Mixtures

Authors: M. S. Woolfson; C. Bigan; J. A. Crowe; B. R. Hayes-Gill
Journal: ISRN Signal Processing
 (2013)

---

## Abstract

A blind source separation method is described to extract sources from data mixtures where the underlying sources are sparse and correlated.  The approach used is to detect and analyze segments of time where one source exists on its own.  The method does not assume independence of sources and probability density functions are not assumed for any of the sources.  A comparison is made between the proposed method and the Fast-ICA and Clusterwise PCA methods. It is shown that the proposed method works best for cases where the underlying sources are strongly correlated because Fast-ICA assumes zero correlation between sources and Clusterwise PCA can be sensitive to overlap between sources.  However, for cases of sources that are sparse and weakly correlated with each other, there is a tendency for Fast-ICA and Clusterwise PCA to have better performances than the proposed method, the reason being that these methods appear to be more robust to changes in input parameters to the algorithms.  In addition, because of the deflationary nature of the proposed method, there is a tendency for estimates to be more affected by noise than Fast-ICA when the number of sources increases.  The paper concludes with a discussion concerning potential applications for the proposed method.


---

## Body

## 1. Introduction

One general problem in signal processing is the extraction of individual source signals{sj[n]} from measurements {zi[n]} that are linear combinations of these sources:
(1)zi[n]=∑j=1NAijsj[n](i=1,2,…,S),
where {Aij} are the mixing coefficients, S is the number of sets of measurement data, and there are N underlying sources. In the case where both the sources and mixing coefficients are unknown, this problem comes under the heading of blind source separation (BSS). There are many applications in this area [1–8].There are various approaches to extracting the underlying sources{si}. For example, in Independent Component Analysis (ICA), the aim is to determine a transformation of the data that maximizes the negentropy [9]. In the approach designed in [10], the aim is to minimize the cross-cumulants.Now in some applications, the underlying sources can be approximated as sparse; that is, each source has nonzero values for a finite segment of time, with the other sources having zero values. In practice, this definition can be approximated by each source being dominant over other sources for at least one segment of time. Various specialized BSS methods have been derived for the case of sparse sources [11–18]. In one approach [12], the sparsity of a signal {x[k]} is defined over K data points by
(2)S(x)=K∑k=1K(x[k])2∑k=1K|x[k]|.
Assuming two sources, the aim is to find a transformation on the data {z1[k],z2[k]} such that the sparsity is maximized. The extension to three or more sources is also discussed. Another approach is based on the fact that regions of the phase plot of the data {zi[k]}, where one source exists on its own, are represented by points that form a straight line. These directions are found using clustering techniques on the phase plot. In [11], Minor Component Analysis (MCA) is applied to the clusters to find elements of the separating matrix B so that the sources can be directly estimated from the data from
(3)s=Bz.
This method will be referred to as the Clusterwise PCA method.In [17] a method to separate out sparse sources from signal mixtures was presented where the mixing coefficients are unknown as well as the sources. This standard Blind Source Separation (BSS) problem was addressed by analysis of localized segments of the phase space trajectory where one source dominates. The method was found to have a comparable performance to the standard Fast-Independent Component Analysis (Fast-ICA) [9] method for mixtures of sparse sources. Another assumption made by the proposed method was that the underlying sources are uncorrelated, an assumption that is also made by the Fast-ICA method [9] and some other, although not all, BSS methods. However, as the method proposed in [17] analyzes local features of phase space, it makes no assumption about the probability distributions of the underlying sources unlike some BSS methods which make such assumptions; hence if the underlying sources had Gaussian distributions this would be no problem for the method in [17], whilst an approach such as Fast-ICA may encounter difficulty in this situation. In this paper, the method presented in [17] is extended to the case where the underlying sources are sparse and correlated. It will be shown that the “prewhitening” stage that is used by many standard BSS and ICA methods is not required here and that the proposed method can work better than the Fast-ICA and Clusterwise PCA methods, for cases where the underlying sources are sparse and strongly correlated. The proposed method can be more sensitive to noise than two other methods tested and can also be more sensitive to changes in input parameters. In particular, for sources that are weakly correlated or uncorrelated, the proposed method offers no advantages compared with Clusterwise PCA and Fast-ICA.


## 2. Theory

### 2.1. Minimum Heading Change (MHC) Method

To motivate the approach used in this paper, an example will be shown where the two sources before mixing are shown in Figure1.Figure 1
Two source signals.It can be seen that source A is zero for sample numbers larger than 19, whilst source B is zero for sample numbers less than 8. Hence, these two sources can be considered as sparse. Now suppose that the two sources are mixed linearly so that the mixed data looks like Figure2.Figure 2
Mixed data.It can be seen that because of the sparsity of the two underlying sources, there are segments of each mixture which depend on one source only; these segments are indicated in Figure2.One can identify segments of data wheres1 exists on its own by plotting data input 2, z2 against data input 1, z1. In Figure 3 this is done for the example in Figure 2; in the following discussion, we associate source B in Figure 1 with s1 and source A with s2.Figure 3
Plot ofz2 against z1.In [17], the estimation of sparse sources from mixtures, looking for straight line segments in the phase plots, was described. A summary of the method is described here.Suppose that there are two mixtures of two sources. Then, from (1),
(4)z1=A11s1+A12s2,z2=A21s1+A22s2,
where the explicit dependences of the measurements and the sources on sample number have been dropped for notational convenience. Source s1 is associated with the direction (A11,A21) in the phase plot with source s2 being associated with the direction (A12,A22). Now for the specific case of mixtures of sources that are sparse, these directions are visible in the phase plot; for example, AB in Figure 3 corresponds to source A on its own in Figure 2, whilst CD (shown dashed in Figure 3) corresponds to source B in Figure 2; this information can be used to estimate the underlying sources.Suppose that sources1 is sparse and that sources s1 and s2 can be associated with unit vectors in the phase plot R^1 and R^2, respectively. One can think of the plot of z2 against z1 as analogous to a radar plot and vectors R^1 and R^2 as “heading” vectors. Any point in the phase plot can be associated with a vector:
(5)z=s1R1+s2R2=∑i=12siRiR^i,
where the vector magnitude Ri(i=1,2) depends on the mixing coefficients and Ri=RiR^i. Note that R^1 and R^2 are not necessarily orthogonal.Now suppose that three consecutive points in the phase plot are detected that are collinear corresponding toR^1. One can make an estimate of the component s1,s~1, by resolving the vector z in the direction of  R^1:
(6)s~1=z·R^1=∑i=12siRiR^i·R^1=s1R1+s2R2R^2·R^1.
The second term in (6) is an error term as it depends on s2.Now as sources1 is associated with a direction R^1 in the phase plot, this source can be represented in the plot as a vector s~1 given by
(7)s~1=s~1R^1.
This source is associated with a constant direction or heading. As s~1 varies with time, so does the magnitude of the vector s~1.Substituting fors~1 from (6) into (7), this vector can be written as
(8)s~1=s1R1R^1+s2R2(R^1·R^2)R^1.
Now from (5), the vector representing the phase plot is given by
(9)z=s1R1R^1+s2R2R^2.
Comparing (8) and (9), it can be seen that the contribution from s2 has been reduced in (8) because in general |R^1·R^2|<1.  It can also be seen in (8) and (9) that the terms in s1 are identical. Hence, one can subtract (8) from (9) as follows:
(10)s~2=z-s~1=s2R2[R^2-(R^1·R^2)R^1].
and one now has an estimate of the vector representing source s2 in the phase plot. The magnitude of this vector is equal to the underlying source s2 to that within a multiplicative constant.Referring again to the analogy with radar plots, one is looking for data segments where the “heading” does not change over three consecutive data points; in practice, because of the effects of noise, one would be looking for aminimum change in heading over these three points; in [17], this method is referred to as the minimum heading change (MHC) method and more details regarding the algorithm for this method can be found in this reference.The MHC method is based on computing a “velocity” vectorv[n], which is given by
(11)v[n]=z[n]-z[n-1]
at each sample point n.In [17] the algorithm used to detect the principal directions {R^i} is described in detail. Also in [17] an analysis was presented for the general N source case where it was shown that the N sources could be estimated as follows:
(12)s~1=C11s1+∑j=2NC1jsj,(13)s~2=C22s2+∑j=3NC2jsj,⋮(14)s~N-1=CN-1,N-1sN-1+CN-1,NsN,(15)s~N=CNNsN,
where Cij is related to the cross-correlation coefficient between source i,si, and source j,sj. In fact, when sources i and j are uncorrelated and prewhitening of data has taken place, it can be shown that Cij = 0 for all i andj such that i≠j; in this particular case, it can be seen from (12) to (15) that
(16)s~i=Ciisi,
and each source is extracted to be within a multiplicative constant.In order to avoid the effects of spurious noise, let us write theN components of the velocity vector as
(17)v[n]=(v1[n],v2[n],…,vN[n]).
A velocity vector is accepted at sample point n if
(18)Vmax[n]≥vth·vmax,
where
(19)Vmax[n]=max{|v1[n]|,|v2[n]|,…,|vN[n]|}0<vth<1,(20)vmax=max{|v[1]|,|v[2]|,…,|v[M]|}
is the maximum value of the magnitude of the velocity vector over all M sample points.


### 2.2. Extension to the MHC Method

In this paper, the algorithm described in [17] is extended to extracting individual sparse sources {s1,s2,…,sN-1} from (12) to (15) for the general case where the sources are correlated. The sparsity of these sources is used to achieve this.


#### 2.2.1. Mixtures of Two Sources

First, we will look at the case of two mixtures of two sources. From (12) and (13) with N=2, it can be seen that, after applying MHC, the two estimated sources are given by
(21)s~1=C11s1+C12s2,(22)s~2=C22s2.
Now, if  s~2 is plotted against  s~1, the following vector in two-dimensional phase space results:
(23)z=(s~1,s~2).Substituting fors~1 from (21) and s~2 from (22) into (23),
(24)z=(C11s1+C12s2,C22s2)=(C11,0)s1+(C12,C22)s2.
Hence,
(25)z=s1R1+s2R2,
where the vectors R1 and R2, corresponding to the underlying sources 1 and 2, respectively, are given by
(26)R1=(C11,0),R2=(C12,C22).
Equation (25) can be rewritten as
(27)z=s1R1R^1+s2R2R^2,
where R^1 and R^2 are unit vectors in the directions of  R1 and R2, respectively.The question now to be asked is whether one can use the fact thats~2 is proportional to s2 only to somehow subtract the contribution from this source to s~1. Now s2 is assumed sparse, so that s1 will exist on its own for a segment of time which will manifest itself in the phase plot as a straight line over at least 3 points. If one assumes that s1 is also sparse, then there will be a straight line segment corresponding to s2 as well.To illustrate this approach, the MHC has been applied to the data mixture in Figure2.The phase plot fors~2 against s~1 for the estimated sources is given in Figure 4.Figure 4
Phase plot ofs~2 against s~1.Two sets of line segments can be seen where the “heading” is constant over three or more data points: between A and B this segment corresponds to vectorR1 in (26) corresponding to samples between 20 and 27 in Figure 1 where estimated source s~1 (corresponding to Source B in Figure 1) exists on its own; this straight line has a zero slope in phase space as estimated source s~2, corresponding to source A in Figure 1, is zero in this data segment.  The trajectory between C and D in Figure 4 corresponds to samples 0 to 7 in Figure 1, where source s2 (corresponding to source A in Figure 1) is on its own (R2 in (26)); from (21) and (22), this source contributes to both s~1 and s~2 so this line has a slope in phase space that is nonzero and of finite value.Suppose that the MHC method is applied with the data inputs, this time beings~1 and s~2. Let us assume that the direction CD is detected at a certain sample point say at n=M; this direction can be associated with the unit vector R^2 in (26).If one takes the scalar product of the vectorz in (27) and R^2 we obtain an “estimate” of the underlying source s2 from
(28)s~2(1)=z·R^2,
where, from (26),
(29)R^2=1R2(C12,C22).
Substituting for z from (27) into (28),
(30)s~2(1)=s1R1(R^1·R^2)+s2R2.s~2(1) is not an exact estimate of s2; this is obtained from (22), but s~2(1) will later be used to estimate source s1.Now the estimate of source 1 is given by (21):
(31)s~1=C11s1+C12s2.
Comparing (30) with (31), an improved estimate of source s1 can be obtained from
(32)s~1(1)=s~1-C12R2s~2(1)=C′s1,
where
(33)C′=(C11-C12R2(R1R^1·R^2)).
Note that the subtraction in (32) can be carried out as we already know the components of R2=(C12,C22) from the detection step and R2=C122+C222; prior knowledge of the coefficient C11 is not required. Another way of expressing (32) is to note that, from (29), C12 is the first component of  R^2:
(34)C12=R^2(1).
Using (28) and (34), (32) can be rewritten:
(35)s~1(1)=s~1-(z·R^2)R^2(1).
From (32), it can be seen that source s1 has now been estimated to be within a scaling constant.This procedure is like operating the MHC in “reverse” starting from the last estimated source,s~2, and then obtaining an improved estimate s~1(1) of source s~1 so this processing step is referred to as “reverse iteration.” The first processing step described in Section 2.1 and [17] will be referred to as the “forward iteration.”In the Reverse processing step, it is necessary to detect segment CD in Figure4; to avoid detecting segment AB, an additional test is carried out. The change in heading is calculated as described in [17], and, in addition, the correlation coefficient between s~1 and s~2 over P data points is computed. In order to guard against detecting segments like AB where a source exists in one component only, the minimum heading change is accepted only if the magnitude of the correlation coefficient is no smaller than a minimum threshold value Corrmin. Once the point where this occurs is found, the heading vector is used in (35) to subtract off the contribution from s2 to the overall phase plot. However, to avoid chance correlations by noise, the correlation is only calculated if the velocity components in phase space have a magnitude that is larger than a minimum threshold, Tmin, given by
(36)Tmin=λthvmax,
where vmax is the maximum velocity component, defined in (20), for vector z in the reverse iteration, and λth is a chosen threshold. Condition (36) is used instead of (18) in the reverse iteration. The following expression is used to calculate the correlation function:
(37)Corr(x,y)=∑i=1P(xi-x-)(yi-y-)∑i=1P(xi-x-)2∑j=1P(yj-y-)2,
where x- and y- are the mean values of x and y defined by
(38)x-=1P∑i=1Pxi,y-=1P∑i=1Pyi.
It has been confirmed in simulations that if this correlation test is not carried out, then the method breaks down as the wrong source can be detected in the reverse iteration. No subtraction takes place if condition (36) is not satisfied.It should be noted that, in many BSS methods, the first step is whitening. In this case, uncorrelated sources have characteristic directions in phase space that are orthogonal. In the MHC method, if whitening is used, uncorrelated components are extracted exactly in forward iteration only asCij=0 for all i and  j such that i≠j in (12) to (15); in this case, there is no need for the reverse iteration to be carried out; in [17], the Gram-Schmidt method was applied to perform the whitening. However, such a procedure, whilst useful, is not necessary for the MHC with reverse iteration, as whitening is not a required preprocessing procedure for this method. In particular, the MHC does not make any assumptions as to how many sources, if any, are correlated. In addition, centring, that is the subtraction of DC components from the data is also not necessary and can lead to undesirable effects where previously uncorrelated sources are assumed to be correlated. However, the problem with nonwhitening is that uncorrelated components cannot be individually estimated in the forward iteration. In this case, one is relying on the reverse iteration to achieve this and this could lead to incomplete separation in the reverse direction for inappropriate choice of input parameters. However, for sources that are correlated strongly enough, there is not much point in prewhitening. In subsequent implementation of the MHC method with reverse iteration, the data will not be whitened initially.


#### 2.2.2. Extension toN Sources

The reverse processing step can be generalized to theN-source case, assuming that one has N data inputs. In particular, it can be shown that individual sources can be extracted even if there is correlation between pairs of sources. A partial analysis for the N-source case is given in Appendix A, with a pseudocode for the complete algorithm being given in Appendix B.


## 2.1. Minimum Heading Change (MHC) Method

To motivate the approach used in this paper, an example will be shown where the two sources before mixing are shown in Figure1.Figure 1
Two source signals.It can be seen that source A is zero for sample numbers larger than 19, whilst source B is zero for sample numbers less than 8. Hence, these two sources can be considered as sparse. Now suppose that the two sources are mixed linearly so that the mixed data looks like Figure2.Figure 2
Mixed data.It can be seen that because of the sparsity of the two underlying sources, there are segments of each mixture which depend on one source only; these segments are indicated in Figure2.One can identify segments of data wheres1 exists on its own by plotting data input 2, z2 against data input 1, z1. In Figure 3 this is done for the example in Figure 2; in the following discussion, we associate source B in Figure 1 with s1 and source A with s2.Figure 3
Plot ofz2 against z1.In [17], the estimation of sparse sources from mixtures, looking for straight line segments in the phase plots, was described. A summary of the method is described here.Suppose that there are two mixtures of two sources. Then, from (1),
(4)z1=A11s1+A12s2,z2=A21s1+A22s2,
where the explicit dependences of the measurements and the sources on sample number have been dropped for notational convenience. Source s1 is associated with the direction (A11,A21) in the phase plot with source s2 being associated with the direction (A12,A22). Now for the specific case of mixtures of sources that are sparse, these directions are visible in the phase plot; for example, AB in Figure 3 corresponds to source A on its own in Figure 2, whilst CD (shown dashed in Figure 3) corresponds to source B in Figure 2; this information can be used to estimate the underlying sources.Suppose that sources1 is sparse and that sources s1 and s2 can be associated with unit vectors in the phase plot R^1 and R^2, respectively. One can think of the plot of z2 against z1 as analogous to a radar plot and vectors R^1 and R^2 as “heading” vectors. Any point in the phase plot can be associated with a vector:
(5)z=s1R1+s2R2=∑i=12siRiR^i,
where the vector magnitude Ri(i=1,2) depends on the mixing coefficients and Ri=RiR^i. Note that R^1 and R^2 are not necessarily orthogonal.Now suppose that three consecutive points in the phase plot are detected that are collinear corresponding toR^1. One can make an estimate of the component s1,s~1, by resolving the vector z in the direction of  R^1:
(6)s~1=z·R^1=∑i=12siRiR^i·R^1=s1R1+s2R2R^2·R^1.
The second term in (6) is an error term as it depends on s2.Now as sources1 is associated with a direction R^1 in the phase plot, this source can be represented in the plot as a vector s~1 given by
(7)s~1=s~1R^1.
This source is associated with a constant direction or heading. As s~1 varies with time, so does the magnitude of the vector s~1.Substituting fors~1 from (6) into (7), this vector can be written as
(8)s~1=s1R1R^1+s2R2(R^1·R^2)R^1.
Now from (5), the vector representing the phase plot is given by
(9)z=s1R1R^1+s2R2R^2.
Comparing (8) and (9), it can be seen that the contribution from s2 has been reduced in (8) because in general |R^1·R^2|<1.  It can also be seen in (8) and (9) that the terms in s1 are identical. Hence, one can subtract (8) from (9) as follows:
(10)s~2=z-s~1=s2R2[R^2-(R^1·R^2)R^1].
and one now has an estimate of the vector representing source s2 in the phase plot. The magnitude of this vector is equal to the underlying source s2 to that within a multiplicative constant.Referring again to the analogy with radar plots, one is looking for data segments where the “heading” does not change over three consecutive data points; in practice, because of the effects of noise, one would be looking for aminimum change in heading over these three points; in [17], this method is referred to as the minimum heading change (MHC) method and more details regarding the algorithm for this method can be found in this reference.The MHC method is based on computing a “velocity” vectorv[n], which is given by
(11)v[n]=z[n]-z[n-1]
at each sample point n.In [17] the algorithm used to detect the principal directions {R^i} is described in detail. Also in [17] an analysis was presented for the general N source case where it was shown that the N sources could be estimated as follows:
(12)s~1=C11s1+∑j=2NC1jsj,(13)s~2=C22s2+∑j=3NC2jsj,⋮(14)s~N-1=CN-1,N-1sN-1+CN-1,NsN,(15)s~N=CNNsN,
where Cij is related to the cross-correlation coefficient between source i,si, and source j,sj. In fact, when sources i and j are uncorrelated and prewhitening of data has taken place, it can be shown that Cij = 0 for all i andj such that i≠j; in this particular case, it can be seen from (12) to (15) that
(16)s~i=Ciisi,
and each source is extracted to be within a multiplicative constant.In order to avoid the effects of spurious noise, let us write theN components of the velocity vector as
(17)v[n]=(v1[n],v2[n],…,vN[n]).
A velocity vector is accepted at sample point n if
(18)Vmax[n]≥vth·vmax,
where
(19)Vmax[n]=max{|v1[n]|,|v2[n]|,…,|vN[n]|}0<vth<1,(20)vmax=max{|v[1]|,|v[2]|,…,|v[M]|}
is the maximum value of the magnitude of the velocity vector over all M sample points.


## 2.2. Extension to the MHC Method

In this paper, the algorithm described in [17] is extended to extracting individual sparse sources {s1,s2,…,sN-1} from (12) to (15) for the general case where the sources are correlated. The sparsity of these sources is used to achieve this.


### 2.2.1. Mixtures of Two Sources

First, we will look at the case of two mixtures of two sources. From (12) and (13) with N=2, it can be seen that, after applying MHC, the two estimated sources are given by
(21)s~1=C11s1+C12s2,(22)s~2=C22s2.
Now, if  s~2 is plotted against  s~1, the following vector in two-dimensional phase space results:
(23)z=(s~1,s~2).Substituting fors~1 from (21) and s~2 from (22) into (23),
(24)z=(C11s1+C12s2,C22s2)=(C11,0)s1+(C12,C22)s2.
Hence,
(25)z=s1R1+s2R2,
where the vectors R1 and R2, corresponding to the underlying sources 1 and 2, respectively, are given by
(26)R1=(C11,0),R2=(C12,C22).
Equation (25) can be rewritten as
(27)z=s1R1R^1+s2R2R^2,
where R^1 and R^2 are unit vectors in the directions of  R1 and R2, respectively.The question now to be asked is whether one can use the fact thats~2 is proportional to s2 only to somehow subtract the contribution from this source to s~1. Now s2 is assumed sparse, so that s1 will exist on its own for a segment of time which will manifest itself in the phase plot as a straight line over at least 3 points. If one assumes that s1 is also sparse, then there will be a straight line segment corresponding to s2 as well.To illustrate this approach, the MHC has been applied to the data mixture in Figure2.The phase plot fors~2 against s~1 for the estimated sources is given in Figure 4.Figure 4
Phase plot ofs~2 against s~1.Two sets of line segments can be seen where the “heading” is constant over three or more data points: between A and B this segment corresponds to vectorR1 in (26) corresponding to samples between 20 and 27 in Figure 1 where estimated source s~1 (corresponding to Source B in Figure 1) exists on its own; this straight line has a zero slope in phase space as estimated source s~2, corresponding to source A in Figure 1, is zero in this data segment.  The trajectory between C and D in Figure 4 corresponds to samples 0 to 7 in Figure 1, where source s2 (corresponding to source A in Figure 1) is on its own (R2 in (26)); from (21) and (22), this source contributes to both s~1 and s~2 so this line has a slope in phase space that is nonzero and of finite value.Suppose that the MHC method is applied with the data inputs, this time beings~1 and s~2. Let us assume that the direction CD is detected at a certain sample point say at n=M; this direction can be associated with the unit vector R^2 in (26).If one takes the scalar product of the vectorz in (27) and R^2 we obtain an “estimate” of the underlying source s2 from
(28)s~2(1)=z·R^2,
where, from (26),
(29)R^2=1R2(C12,C22).
Substituting for z from (27) into (28),
(30)s~2(1)=s1R1(R^1·R^2)+s2R2.s~2(1) is not an exact estimate of s2; this is obtained from (22), but s~2(1) will later be used to estimate source s1.Now the estimate of source 1 is given by (21):
(31)s~1=C11s1+C12s2.
Comparing (30) with (31), an improved estimate of source s1 can be obtained from
(32)s~1(1)=s~1-C12R2s~2(1)=C′s1,
where
(33)C′=(C11-C12R2(R1R^1·R^2)).
Note that the subtraction in (32) can be carried out as we already know the components of R2=(C12,C22) from the detection step and R2=C122+C222; prior knowledge of the coefficient C11 is not required. Another way of expressing (32) is to note that, from (29), C12 is the first component of  R^2:
(34)C12=R^2(1).
Using (28) and (34), (32) can be rewritten:
(35)s~1(1)=s~1-(z·R^2)R^2(1).
From (32), it can be seen that source s1 has now been estimated to be within a scaling constant.This procedure is like operating the MHC in “reverse” starting from the last estimated source,s~2, and then obtaining an improved estimate s~1(1) of source s~1 so this processing step is referred to as “reverse iteration.” The first processing step described in Section 2.1 and [17] will be referred to as the “forward iteration.”In the Reverse processing step, it is necessary to detect segment CD in Figure4; to avoid detecting segment AB, an additional test is carried out. The change in heading is calculated as described in [17], and, in addition, the correlation coefficient between s~1 and s~2 over P data points is computed. In order to guard against detecting segments like AB where a source exists in one component only, the minimum heading change is accepted only if the magnitude of the correlation coefficient is no smaller than a minimum threshold value Corrmin. Once the point where this occurs is found, the heading vector is used in (35) to subtract off the contribution from s2 to the overall phase plot. However, to avoid chance correlations by noise, the correlation is only calculated if the velocity components in phase space have a magnitude that is larger than a minimum threshold, Tmin, given by
(36)Tmin=λthvmax,
where vmax is the maximum velocity component, defined in (20), for vector z in the reverse iteration, and λth is a chosen threshold. Condition (36) is used instead of (18) in the reverse iteration. The following expression is used to calculate the correlation function:
(37)Corr(x,y)=∑i=1P(xi-x-)(yi-y-)∑i=1P(xi-x-)2∑j=1P(yj-y-)2,
where x- and y- are the mean values of x and y defined by
(38)x-=1P∑i=1Pxi,y-=1P∑i=1Pyi.
It has been confirmed in simulations that if this correlation test is not carried out, then the method breaks down as the wrong source can be detected in the reverse iteration. No subtraction takes place if condition (36) is not satisfied.It should be noted that, in many BSS methods, the first step is whitening. In this case, uncorrelated sources have characteristic directions in phase space that are orthogonal. In the MHC method, if whitening is used, uncorrelated components are extracted exactly in forward iteration only asCij=0 for all i and  j such that i≠j in (12) to (15); in this case, there is no need for the reverse iteration to be carried out; in [17], the Gram-Schmidt method was applied to perform the whitening. However, such a procedure, whilst useful, is not necessary for the MHC with reverse iteration, as whitening is not a required preprocessing procedure for this method. In particular, the MHC does not make any assumptions as to how many sources, if any, are correlated. In addition, centring, that is the subtraction of DC components from the data is also not necessary and can lead to undesirable effects where previously uncorrelated sources are assumed to be correlated. However, the problem with nonwhitening is that uncorrelated components cannot be individually estimated in the forward iteration. In this case, one is relying on the reverse iteration to achieve this and this could lead to incomplete separation in the reverse direction for inappropriate choice of input parameters. However, for sources that are correlated strongly enough, there is not much point in prewhitening. In subsequent implementation of the MHC method with reverse iteration, the data will not be whitened initially.


### 2.2.2. Extension toN Sources

The reverse processing step can be generalized to theN-source case, assuming that one has N data inputs. In particular, it can be shown that individual sources can be extracted even if there is correlation between pairs of sources. A partial analysis for the N-source case is given in Appendix A, with a pseudocode for the complete algorithm being given in Appendix B.


## 2.2.1. Mixtures of Two Sources

First, we will look at the case of two mixtures of two sources. From (12) and (13) with N=2, it can be seen that, after applying MHC, the two estimated sources are given by
(21)s~1=C11s1+C12s2,(22)s~2=C22s2.
Now, if  s~2 is plotted against  s~1, the following vector in two-dimensional phase space results:
(23)z=(s~1,s~2).Substituting fors~1 from (21) and s~2 from (22) into (23),
(24)z=(C11s1+C12s2,C22s2)=(C11,0)s1+(C12,C22)s2.
Hence,
(25)z=s1R1+s2R2,
where the vectors R1 and R2, corresponding to the underlying sources 1 and 2, respectively, are given by
(26)R1=(C11,0),R2=(C12,C22).
Equation (25) can be rewritten as
(27)z=s1R1R^1+s2R2R^2,
where R^1 and R^2 are unit vectors in the directions of  R1 and R2, respectively.The question now to be asked is whether one can use the fact thats~2 is proportional to s2 only to somehow subtract the contribution from this source to s~1. Now s2 is assumed sparse, so that s1 will exist on its own for a segment of time which will manifest itself in the phase plot as a straight line over at least 3 points. If one assumes that s1 is also sparse, then there will be a straight line segment corresponding to s2 as well.To illustrate this approach, the MHC has been applied to the data mixture in Figure2.The phase plot fors~2 against s~1 for the estimated sources is given in Figure 4.Figure 4
Phase plot ofs~2 against s~1.Two sets of line segments can be seen where the “heading” is constant over three or more data points: between A and B this segment corresponds to vectorR1 in (26) corresponding to samples between 20 and 27 in Figure 1 where estimated source s~1 (corresponding to Source B in Figure 1) exists on its own; this straight line has a zero slope in phase space as estimated source s~2, corresponding to source A in Figure 1, is zero in this data segment.  The trajectory between C and D in Figure 4 corresponds to samples 0 to 7 in Figure 1, where source s2 (corresponding to source A in Figure 1) is on its own (R2 in (26)); from (21) and (22), this source contributes to both s~1 and s~2 so this line has a slope in phase space that is nonzero and of finite value.Suppose that the MHC method is applied with the data inputs, this time beings~1 and s~2. Let us assume that the direction CD is detected at a certain sample point say at n=M; this direction can be associated with the unit vector R^2 in (26).If one takes the scalar product of the vectorz in (27) and R^2 we obtain an “estimate” of the underlying source s2 from
(28)s~2(1)=z·R^2,
where, from (26),
(29)R^2=1R2(C12,C22).
Substituting for z from (27) into (28),
(30)s~2(1)=s1R1(R^1·R^2)+s2R2.s~2(1) is not an exact estimate of s2; this is obtained from (22), but s~2(1) will later be used to estimate source s1.Now the estimate of source 1 is given by (21):
(31)s~1=C11s1+C12s2.
Comparing (30) with (31), an improved estimate of source s1 can be obtained from
(32)s~1(1)=s~1-C12R2s~2(1)=C′s1,
where
(33)C′=(C11-C12R2(R1R^1·R^2)).
Note that the subtraction in (32) can be carried out as we already know the components of R2=(C12,C22) from the detection step and R2=C122+C222; prior knowledge of the coefficient C11 is not required. Another way of expressing (32) is to note that, from (29), C12 is the first component of  R^2:
(34)C12=R^2(1).
Using (28) and (34), (32) can be rewritten:
(35)s~1(1)=s~1-(z·R^2)R^2(1).
From (32), it can be seen that source s1 has now been estimated to be within a scaling constant.This procedure is like operating the MHC in “reverse” starting from the last estimated source,s~2, and then obtaining an improved estimate s~1(1) of source s~1 so this processing step is referred to as “reverse iteration.” The first processing step described in Section 2.1 and [17] will be referred to as the “forward iteration.”In the Reverse processing step, it is necessary to detect segment CD in Figure4; to avoid detecting segment AB, an additional test is carried out. The change in heading is calculated as described in [17], and, in addition, the correlation coefficient between s~1 and s~2 over P data points is computed. In order to guard against detecting segments like AB where a source exists in one component only, the minimum heading change is accepted only if the magnitude of the correlation coefficient is no smaller than a minimum threshold value Corrmin. Once the point where this occurs is found, the heading vector is used in (35) to subtract off the contribution from s2 to the overall phase plot. However, to avoid chance correlations by noise, the correlation is only calculated if the velocity components in phase space have a magnitude that is larger than a minimum threshold, Tmin, given by
(36)Tmin=λthvmax,
where vmax is the maximum velocity component, defined in (20), for vector z in the reverse iteration, and λth is a chosen threshold. Condition (36) is used instead of (18) in the reverse iteration. The following expression is used to calculate the correlation function:
(37)Corr(x,y)=∑i=1P(xi-x-)(yi-y-)∑i=1P(xi-x-)2∑j=1P(yj-y-)2,
where x- and y- are the mean values of x and y defined by
(38)x-=1P∑i=1Pxi,y-=1P∑i=1Pyi.
It has been confirmed in simulations that if this correlation test is not carried out, then the method breaks down as the wrong source can be detected in the reverse iteration. No subtraction takes place if condition (36) is not satisfied.It should be noted that, in many BSS methods, the first step is whitening. In this case, uncorrelated sources have characteristic directions in phase space that are orthogonal. In the MHC method, if whitening is used, uncorrelated components are extracted exactly in forward iteration only asCij=0 for all i and  j such that i≠j in (12) to (15); in this case, there is no need for the reverse iteration to be carried out; in [17], the Gram-Schmidt method was applied to perform the whitening. However, such a procedure, whilst useful, is not necessary for the MHC with reverse iteration, as whitening is not a required preprocessing procedure for this method. In particular, the MHC does not make any assumptions as to how many sources, if any, are correlated. In addition, centring, that is the subtraction of DC components from the data is also not necessary and can lead to undesirable effects where previously uncorrelated sources are assumed to be correlated. However, the problem with nonwhitening is that uncorrelated components cannot be individually estimated in the forward iteration. In this case, one is relying on the reverse iteration to achieve this and this could lead to incomplete separation in the reverse direction for inappropriate choice of input parameters. However, for sources that are correlated strongly enough, there is not much point in prewhitening. In subsequent implementation of the MHC method with reverse iteration, the data will not be whitened initially.


## 2.2.2. Extension toN Sources

The reverse processing step can be generalized to theN-source case, assuming that one has N data inputs. In particular, it can be shown that individual sources can be extracted even if there is correlation between pairs of sources. A partial analysis for the N-source case is given in Appendix A, with a pseudocode for the complete algorithm being given in Appendix B.


## 3. Results

### 3.1. Two Sources, Low Sparsity

The application of this procedure to the extraction of sources from the mixed signals in Figure2 is now described; this is an example of two sources which are sparse over relatively few samples. The outputs from the “forward” iteration, s~1 and s~2, are taken and the reverse iteration of the MHC is applied to extract estimates s~1(1) and  s~2(1)  of the underlying sources s1 and s2, respectively. In this simulation, vth=0.2 (18), λth=10-5 (36), and the minimum allowed value of Corr(x,y) in (37) is Corrmin=0.9, with P=3 in (37); all these parameters were chosen by trial and error.The figure of merit used is to correlate each estimate with each underlying source and make the association that has the largest correlation value—the CORRCOEF function in MATLAB is used to compute this parameter. The correlation coefficients, over the whole signals, between the estimated and actual sources for each methods are shown in Table1. For comparison, the corresponding calculations are carried out using (i) the MHC method implementing just forward iteration with the Gram-Schmidt whitening [17], (ii) the Fast-ICA software Version 2.5 using the deflation method with the Gaussian nonlinearity [19] and (iii) the Clusterwise PCA method [11]. The MHC, with reverse iteration, has estimated both sources to be within a scaling constant. Exact estimation of the underlying sources has not taken place with the Fast-ICA and Clusterwise PCA methods. It should be noted that the cross-correlation between the actual sources in Figure 1 is 0.361 and is not negligible, which explains why the Fast- ICA method cannot extract any of the sources exactly for this example. The Clusterwise PCA method does not extract the sources exactly because the sources are sparse for only a few data points, so do not form obvious clusters of points. As predicted using (10), when using the MHC with just forward iteration [17], only one of the sources can be estimated exactly.Table 1
Correlation coefficients between estimated and actual sources.
Method
Source A
Source B
Clusterwise PCA [11]
0.9994
0.6310
Fast-ICA [19]
0.9961
0.8970
MHC (Forward) [17]
1
0.9854
MHC (Reverse)
1
1The robustness to noise for the MHC method is now investigated for the mixed data shown in Figure2. The noise, which is modeled as Gaussian, is added to the samples of each data input after the mixing has taken place.At theqth Monte Carlo run, the absolute value of the correlation coefficient is calculated between each source and estimate, over the whole data segment. An association between source and estimate is found by looking for the estimate that has highest correlation with that source.Calculations are carried out for 2000 Monte Carlo runs and a noise standard deviation of 0.005 is used; this value corresponds to 1% of the peak value ofz1 in Figure 2 and 0.5% of the magnitude of the peak value of z2. The absolute value of the correlation coefficients between the estimated and actual sources, averaged over all Monte Carlo runs, is found to be 0.991 and 0.817 for sources A and B, respectively, meaning that there is a tendency for the estimates of source A to be superior to those for source B. This can be explained as being due to the deflationary nature of the algorithm, where the various sources are estimated and subtracted one by one. For this particular set of simulations, it is found that source B tends to be detected first in the forward iteration. For those Monte Carlo runs where source B is detected first in the forward iteration, we can identify s1 with source B and s2 with source A, and so, neglecting noise, the estimates are given from (21) and (22) by
(39)s~B=C11sB+C12sA,(40)s~A=C22sA,
where the subscripts “A” and “B” correspond to sources A and B, respectively. It can be seen from (40) that s~A, corresponding to source A, is estimated to be within a scaling constant. From (32), (39), and (40), source sB is estimated in the reverse iteration according to
(41)s~B(1)=s~B-C12RAs~A(1),
where
(42)RA=(C12,C22).
From (39), (40), and (41), it can be seen that the reverse iteration step requires an estimate of the heading vector RA in (42) corresponding to source A; this estimate will be affected by noise. The estimate s~B(1) in (41) will depend on this noisy estimate of RA as well as on the estimates s~A and s~B in the forward iteration, which will both also be affected by noise; thus these combined errors will cause the estimate,  s~B(1), of source B in the reverse iteration to be more affected by noise than the estimate of source A, s~A, in the forward iteration. This sensitivity to noise will limit the number of underlying sources that could be extracted with this method because, due to the deflationary nature of the method, the later a source is detected in the reverse iteration, then the more it will be affected by noise.


### 3.2. Simulations for Two Sources: High Sparsity

In Section2, the proposed method was tested on mixtures of sources, where each source exists on its own for a relatively short segment of time compared with the whole data interval. In this section, the case where the underlying sources are sparser than the sources in Figure 1 are shown. A comparison is made between the performance of the MHC with reverse iteration and both the Clusterwise PCA and Fast-ICA methods. There are two datasets, each dataset consisting of a mixture of two periodic signals and each period consisting of Gaussian functions given as follows:
(43)si(t)=aiexp[-(t-t0i)22σi2],
where ai is the amplitude, σi is the standard deviation, and t0i is the time shift of the first period. This model was used in [17].This Gaussian is repeated with a periodTi. 10 seconds of data are analyzed. The sampling frequency is set at 250 Hz.The parameters used for the two source signals are as follows.Source 1:
a
1 various values,
t
01
=
0.1 s,
σ
1
=
12.5 ms,
T
1
=
1 s.
Source 2:
a
2 = various values,
t
02
=
0.22 s,
σ
2
=
6.25 ms,
T
2
=
0.688 s.We now look at the case for two signals each composed of a mixture of sources 1 and 2 above with the following amplitudes:measurement set 1:a1=6.5, a2=1,measurement set 2:a1=3, a2=1,N=2 in (1).The two mixed signals are shown in Figure 5.Simulated mixed signals.
(a)
(b)These signals could represent, for example, a simplified model of the ECG measurements taken from the abdomen of an expectant mother, where source 1 represents a contribution from the maternal ECG and source 2 a contribution from the fetal ECG. In this paper, this dataset is separated into two segments: segment 1 where there is no coincidence between the two sources, between 0.5 s and 2.5 s in Figure5, and segment 2 where coincidence between the sources occurs, between 7 s and 9 s in Figure 5—this could, for example, represent the maternal/fetal coincidence in abdominal ECG data.The following parameters in the MHC are chosen after some experimentation:Corrmin=0.999, λth=10-5 and vth=0.5. Also, in the reverse iteration, it is found that best results are obtained by using P=5 in (37).The Fast-ICA method is run using the deflation method with the Gaussian nonlinearity. The parameters used in the Clusterwise PCA method are the same as listed in [11].The various methods are assessed by computing the correlation coefficient between the estimated and actual sources using the same methodology used for the analysis in Table1.The results of applying the Clusterwise PCA, Fast-ICA, and MHC (reverse) methods to segments 1 and 2 are shown in Table2.Table 2
Correlation between estimate and source for (a) Segment 1 and (b) Segment 2.
(a)
Method
Source 1
Source 2
MHC (Reverse)
1
1
Clusterwise PCA
1
1
Fast-ICA
1
0.9993
(b)
Method
Source 1
Source 2
MHC (Reverse)
1
1
Clusterwise PCA
0.6492
0.9993
Fast-ICA
0.9479
0.9997In Table2, it can be seen that for segment 1, both the Clusterwise PCA and MHC reverse iteration methods are able to extract the sources exactly. The Fast-ICA method subtracts the average from the original data and this introduces correlation between the assumed sources and so the correlation coefficients between estimated and actual sources are less than 1 for one of the sources. For segment 2, only the MHC method with reverse iteration manages to extract the two sources to be within a scaling constant. The Fast-ICA cannot extract the sources exactly because there is now significant correlation between the two underlying sources due to the coincidence occurring between the two sources. The Clusterwise PCA method also breaks down when estimating source 1 because of the overlap between sources 1 and 2; to see this further, in Figure 6, the phase plot of the data for segment 2 is shown: 1 and 2 refer to the contributions to the phase plot for sources 1 and 2, respectively. “12” is the contribution to the phase plot from the overlapping of sources 1 and 2 that occurs in segment 2; this feature is of significant size compared with the features representing source 1 only and source 2 only. The Clusterwise PCA method attempts to group the plot into two “clusters.” The contribution 12 is clustered with 1, which degrades the estimate of this source. However, this contribution has less of an effect on the cluster corresponding to source 2, and hence this source is estimated much better, as observed in Table 2.Figure 6
Phase plot for segment 2. “1” and “2” refer to intervals in the phase plot corresponding to sources 1 and 2, respectively, and the label “12” refers to coincidence between the two sources.The robustness to the noise of the proposed method is now looked for segment 1 where there are no coincidences between the two underlying sources. The amplitude of source 2 for each set of data is equal to 1; the Gaussian noise is added with standard deviations varying from 0 to 30% of the amplitude of source 2. 2000 Monte Carlo runs are used.The absolute value of the correlation coefficient between each source and closest estimate is averaged over all Monte Carlo runs. A comparison is made between the performances of the proposed method, the Fast-ICA and the Clusterwise PCA methods. In Figure7(a), the mean correlation is plotted as a function of the standard deviation of noise for source 1, and in Figure 7(b) the corresponding plot is shown for source 2. Note that the range of noise standard deviations used in the plots are different as the extraction of the smaller source 2 breaks down at smaller noise levels.Mean correlation coefficient as a function of noise standard deviation: (a) source 1 (b) and source 2. Full line: Fast-ICA, dashed line, Clusterwise PCA; dotted line: MHC (reverse).
(a)
(b)In Figure7(a), it can be seen that MHC and Fast-ICA have comparable performances when estimating source 1, with the latter method having slightly more robustness to noise. The Clusterwise PCA method is more sensitive to noise, particularly for noise standard deviation larger than 0.2; the problem here is that if the noise is too large, then points in phase space will be associated with the wrong cluster. All the methods give a comparable performance when estimating source 1 for noise standard deviations up to 0.05, corresponding to 5% of the amplitude of source 2 in each set of data.In Figure7(b), when estimating source 2, Fast-ICA and Clusterwise PCA have indistinguishable performances for the range of standard deviations of noise used (0 to 5% of the amplitude of source 2). MHC is more sensitive to noise than the other two methods for noise standard deviations greater than about 1% of the amplitude of source 2. The reason for this is that Fast-ICA and Clusterwise PCA process the data globally, whilst MHC does this locally, so is more susceptible to noise.Now it is of interest to ask whether, when using the MHC, the reverse iteration method is more susceptible to noise than if the forward estimation only is used. In Figure8, the mean correlation coefficient between estimate and source for the forward MHC with no whitening is compared with MHC with the reverse iteration for source 2. We can see that reverse iteration does not contribute significantly more errors to the estimation of source 2. Hence most of the errors occurring in the MHC reverse method are occurring in the forward iteration.Figure 8
Extraction of source 2: full line: MHC with reverse iteration, dashed line: MHC with forward iteration only.For this particular simulation, the mixing matrix is given byA1=(6.5131) and the phase space plot for the mixed signals is shown in Figure 6. The directions corresponding to sources 1 and 2 are very clear in the directions (6.5,3) and (1,1). However, a much stricter test for the reverse MHC method is when these two directions are much closer. Simulations have also been carried out for the following mixing matrices:
(44)A2=(3.2131),A3=(3.02131)
for segments of time 0.5 to 2.5 seconds (no coincidence of sources 1 and 2) and 7 s to 9 s where there is coincidence. For 0.5 to 2.5 seconds, it is found that all three methods tested have a similar performance as when using A1, namely that Clusterwise PCA and MHC reverse extracts the sources exactly, and Fast-ICA extracts the sources with correlation coefficients 0.999 for source 1 and 1 for source 2. A stricter test of the performance of these methods when using mixing matrices A2 and A3 is to look at the segment of data, where there is coincidence of sources between 7 s and 9 s.  The performances of these methods for this data segment are shown in Tables 3 and 4.Table 3
Correlation between estimate and source for mixture matrixA2.
Method
Source 1
Source 2
MHC (Reverse)
1
1
Clusterwise PCA
0.7970
0.9894
Fast-ICA
0.9480
0.9997Table 4
Correlation between estimate and source for mixture matrixA3.
Method
Source 1
Source 2
MHC (Reverse)
1
1
Clusterwise PCA
0.7133
0.9992
Fast-ICA
0.9480
0.9997The MHC (reverse) method still gives perfect extraction of sources to be  within a scaling constant and has the best performance of the three methods tested for this particular dataset.


### 3.3. Two Sources, Varying Sparsity

It is of interest to investigate the performance of the proposed algorithm when processing mixtures of sparse sources with different “sparsity.” At one extreme, one has sources that completely overlap and where there are no segments where one source exists on its own. At the other extreme, one could have sources that do not overlap. An intermediate situation is where there is partial overlap between sources.In the next set of simulations, two mixtures of two sources are shown, where each source is generated from a uniform random number generator with values between −1 and +1; each source consists of 100 samples.The two sources are mixed using the randomly selected mixing matrix:(45)A=(0.5745807850.166294188-0.668082703-0.933171682).
Next, the second source is shifted by L samples with respect to the first source and the mixing matrix (45) applied to the two sources to produce mixtures with 100 + L samples. For each value of L, the MHC algorithm with reverse iteration is applied; the following parameters are found empirically to work for this data: vth=0.1 (18), λth=10-5 (36), and the minimum allowed value of  Corr(x,y) in (37) is Corrmin=0.9, with P=3 in (37).In Figure9, the correlation coefficient between estimated sources and closest actual sources, for sources 1 and 2, is plotted as a function of L for the Fast ICA, Clusterwise PCA, and MHC reverse methods.Mean correlation coefficient as a function of shift,L, between sources: (a) source 1 and (b) source 2. Full line: Fast-ICA; dashed line, Clusterwise PCA; dotted line: MHC (reverse).
(a)
(b)The following observations can be made.(i)
Clusterwise PCA extracts both sources with correlation coefficients larger than 0.98 forL≥70.(ii)
Fast-ICA should be expected to work best forL=0 and L around 100 because in these two extremes the sources are uncorrelated. At some intermediate values of L,where there is some correlation, it is expected that Fast-ICA should not work so well, and this is observed in this example for L between 30 and 50 for source 1. However, the estimates for source 2 have correlation coefficients larger than 0.98 for all values of L.(iii)
Fast-ICA has the best overall performance for all values ofL. However, if L≥3, then MHC (reverse) has the best performance as this method requires only three samples to work; Clusterwise PCA requires larger values of L to work because it performs a global clustering of the data.(iv)
The Fast-ICA algorithm is based on iteratively updating weights applied to the data. The initial values of these weights are chosen randomly; hence, one can obtain different results by rerunning the Fast-ICA algorithm and it is found that reinitializing two or more times can lead to better results for some values ofL. When using Fast ICA, the symmetric approach is found to work better than deflation for these signals. Changing the nonlinearity from gauss to pow 3 and sometimes changing the step size in the iteration of weights can lead to better results; in Figure 9, the best estimates, changing the input parameters, are used to construct the performance curves for Fast ICA. So, for this mixture of sources, the estimation using Fast-ICA can be very sensitive to the approaches used and to input parameters; however, for MHC (reverse), the same set of parameters works for all values of  L≥3.


### 3.4. 4-Lead ECG Data

In order to compare the performances on actual data, of the Fast-ICA, MHC, and Clusterwise PCA methods, data taken from the Daisy Database [20] will be analyzed. These data consists of ECG signals taken from an expectant mother. The data consists of 8 leads, 1 to 5 being abdominal and 6 to 8 thoracic. There are 2500 samples of data; there is some uncertainty about the sampling frequency, also pointed out in [10], but it is probably 250 Hz.These data were analysed in [17], where the MHC with forward iteration was applied to the data. Successful extractions of the fetal and maternal signals were achieved.In this paper, the MHC with reverse iteration is also applied. 4 leads are chosen for analysis: leads 1, 2, 3, and 4—this is a more challenging problem in extracting the fetal signal than when using 8 leads. In [17], all the data samples are analyzed; in this paper, data used is between sample numbers 1451 to 1751 where there is an overlap between a maternal and fetal complex as shown in Figure 10(a) where the maternal and fetal R-waves are marked by “M” and “F,” respectively.Abdominal and thoracic ECG data [20]. (a) Lead 1, (b) extraction of fetal ECG using Fast-ICA, (c) extraction of fetal ECG using MHC (forward iteration, Gram-Schmidt’s initialization), (d) Extraction of fetal ECG using MHC (Reverse iteration), (e) Clusterwise PCA.
(a)
(b)
(c)
(d)
(e)A comparison is made between the Fast-ICA (deflation, Gaussian), MHC forward iteration (with the Gram-Schmidt initialization), MHC with reverse iteration with no whitening of data, and Clusterwise PCA. The input parameters for the MHC calculations werevth=0.5 in (18), Corrmin=0.999, P=5 in (37), and  λth=10-5 in (36).For each method applied, the output that is closest to the expected fetal signal is chosen. The results for the Fast-ICA, MHC (forward), MHC (reverse), and Clusterwise PCA methods are shown in Figures10(b), 10(c), 10(d), and 10(e), respectively.It can be seen that all four methods give comparable results for the extraction of the fetal sources.Comparing Figures10(c) and 10(d), it is observed that the MHC with forward iteration gives similar results for the fetal source as the MHC with reverse iteration. This is due to the fact that despite the maternal/fetal coincidence in the data, there is still low correlation between the maternal and fetal sources. Hence when whitening is carried out, the vectors representing the maternal and fetal sources in the resulting phase plot are almost perpendicular and, from (16), this means that the extracted sources in the forward iteration are approximately a scaled version of the actual sources. The MHC with reverse iteration does not whiten the data initially but is still able to extract the fetal source with comparable quality to the other three methods applied to this data. It is interesting to note that noise is not having a deleterious effect on the quality of the estimated fetal sources in Figures 10(b), 10(c), 10(d), and 10(e); for this particular set of data, the anti-aliasing filter has already filtered off a significant amount of noise.


## 3.1. Two Sources, Low Sparsity

The application of this procedure to the extraction of sources from the mixed signals in Figure2 is now described; this is an example of two sources which are sparse over relatively few samples. The outputs from the “forward” iteration, s~1 and s~2, are taken and the reverse iteration of the MHC is applied to extract estimates s~1(1) and  s~2(1)  of the underlying sources s1 and s2, respectively. In this simulation, vth=0.2 (18), λth=10-5 (36), and the minimum allowed value of Corr(x,y) in (37) is Corrmin=0.9, with P=3 in (37); all these parameters were chosen by trial and error.The figure of merit used is to correlate each estimate with each underlying source and make the association that has the largest correlation value—the CORRCOEF function in MATLAB is used to compute this parameter. The correlation coefficients, over the whole signals, between the estimated and actual sources for each methods are shown in Table1. For comparison, the corresponding calculations are carried out using (i) the MHC method implementing just forward iteration with the Gram-Schmidt whitening [17], (ii) the Fast-ICA software Version 2.5 using the deflation method with the Gaussian nonlinearity [19] and (iii) the Clusterwise PCA method [11]. The MHC, with reverse iteration, has estimated both sources to be within a scaling constant. Exact estimation of the underlying sources has not taken place with the Fast-ICA and Clusterwise PCA methods. It should be noted that the cross-correlation between the actual sources in Figure 1 is 0.361 and is not negligible, which explains why the Fast- ICA method cannot extract any of the sources exactly for this example. The Clusterwise PCA method does not extract the sources exactly because the sources are sparse for only a few data points, so do not form obvious clusters of points. As predicted using (10), when using the MHC with just forward iteration [17], only one of the sources can be estimated exactly.Table 1
Correlation coefficients between estimated and actual sources.
Method
Source A
Source B
Clusterwise PCA [11]
0.9994
0.6310
Fast-ICA [19]
0.9961
0.8970
MHC (Forward) [17]
1
0.9854
MHC (Reverse)
1
1The robustness to noise for the MHC method is now investigated for the mixed data shown in Figure2. The noise, which is modeled as Gaussian, is added to the samples of each data input after the mixing has taken place.At theqth Monte Carlo run, the absolute value of the correlation coefficient is calculated between each source and estimate, over the whole data segment. An association between source and estimate is found by looking for the estimate that has highest correlation with that source.Calculations are carried out for 2000 Monte Carlo runs and a noise standard deviation of 0.005 is used; this value corresponds to 1% of the peak value ofz1 in Figure 2 and 0.5% of the magnitude of the peak value of z2. The absolute value of the correlation coefficients between the estimated and actual sources, averaged over all Monte Carlo runs, is found to be 0.991 and 0.817 for sources A and B, respectively, meaning that there is a tendency for the estimates of source A to be superior to those for source B. This can be explained as being due to the deflationary nature of the algorithm, where the various sources are estimated and subtracted one by one. For this particular set of simulations, it is found that source B tends to be detected first in the forward iteration. For those Monte Carlo runs where source B is detected first in the forward iteration, we can identify s1 with source B and s2 with source A, and so, neglecting noise, the estimates are given from (21) and (22) by
(39)s~B=C11sB+C12sA,(40)s~A=C22sA,
where the subscripts “A” and “B” correspond to sources A and B, respectively. It can be seen from (40) that s~A, corresponding to source A, is estimated to be within a scaling constant. From (32), (39), and (40), source sB is estimated in the reverse iteration according to
(41)s~B(1)=s~B-C12RAs~A(1),
where
(42)RA=(C12,C22).
From (39), (40), and (41), it can be seen that the reverse iteration step requires an estimate of the heading vector RA in (42) corresponding to source A; this estimate will be affected by noise. The estimate s~B(1) in (41) will depend on this noisy estimate of RA as well as on the estimates s~A and s~B in the forward iteration, which will both also be affected by noise; thus these combined errors will cause the estimate,  s~B(1), of source B in the reverse iteration to be more affected by noise than the estimate of source A, s~A, in the forward iteration. This sensitivity to noise will limit the number of underlying sources that could be extracted with this method because, due to the deflationary nature of the method, the later a source is detected in the reverse iteration, then the more it will be affected by noise.


## 3.2. Simulations for Two Sources: High Sparsity

In Section2, the proposed method was tested on mixtures of sources, where each source exists on its own for a relatively short segment of time compared with the whole data interval. In this section, the case where the underlying sources are sparser than the sources in Figure 1 are shown. A comparison is made between the performance of the MHC with reverse iteration and both the Clusterwise PCA and Fast-ICA methods. There are two datasets, each dataset consisting of a mixture of two periodic signals and each period consisting of Gaussian functions given as follows:
(43)si(t)=aiexp[-(t-t0i)22σi2],
where ai is the amplitude, σi is the standard deviation, and t0i is the time shift of the first period. This model was used in [17].This Gaussian is repeated with a periodTi. 10 seconds of data are analyzed. The sampling frequency is set at 250 Hz.The parameters used for the two source signals are as follows.Source 1:
a
1 various values,
t
01
=
0.1 s,
σ
1
=
12.5 ms,
T
1
=
1 s.
Source 2:
a
2 = various values,
t
02
=
0.22 s,
σ
2
=
6.25 ms,
T
2
=
0.688 s.We now look at the case for two signals each composed of a mixture of sources 1 and 2 above with the following amplitudes:measurement set 1:a1=6.5, a2=1,measurement set 2:a1=3, a2=1,N=2 in (1).The two mixed signals are shown in Figure 5.Simulated mixed signals.
(a)
(b)These signals could represent, for example, a simplified model of the ECG measurements taken from the abdomen of an expectant mother, where source 1 represents a contribution from the maternal ECG and source 2 a contribution from the fetal ECG. In this paper, this dataset is separated into two segments: segment 1 where there is no coincidence between the two sources, between 0.5 s and 2.5 s in Figure5, and segment 2 where coincidence between the sources occurs, between 7 s and 9 s in Figure 5—this could, for example, represent the maternal/fetal coincidence in abdominal ECG data.The following parameters in the MHC are chosen after some experimentation:Corrmin=0.999, λth=10-5 and vth=0.5. Also, in the reverse iteration, it is found that best results are obtained by using P=5 in (37).The Fast-ICA method is run using the deflation method with the Gaussian nonlinearity. The parameters used in the Clusterwise PCA method are the same as listed in [11].The various methods are assessed by computing the correlation coefficient between the estimated and actual sources using the same methodology used for the analysis in Table1.The results of applying the Clusterwise PCA, Fast-ICA, and MHC (reverse) methods to segments 1 and 2 are shown in Table2.Table 2
Correlation between estimate and source for (a) Segment 1 and (b) Segment 2.
(a)
Method
Source 1
Source 2
MHC (Reverse)
1
1
Clusterwise PCA
1
1
Fast-ICA
1
0.9993
(b)
Method
Source 1
Source 2
MHC (Reverse)
1
1
Clusterwise PCA
0.6492
0.9993
Fast-ICA
0.9479
0.9997In Table2, it can be seen that for segment 1, both the Clusterwise PCA and MHC reverse iteration methods are able to extract the sources exactly. The Fast-ICA method subtracts the average from the original data and this introduces correlation between the assumed sources and so the correlation coefficients between estimated and actual sources are less than 1 for one of the sources. For segment 2, only the MHC method with reverse iteration manages to extract the two sources to be within a scaling constant. The Fast-ICA cannot extract the sources exactly because there is now significant correlation between the two underlying sources due to the coincidence occurring between the two sources. The Clusterwise PCA method also breaks down when estimating source 1 because of the overlap between sources 1 and 2; to see this further, in Figure 6, the phase plot of the data for segment 2 is shown: 1 and 2 refer to the contributions to the phase plot for sources 1 and 2, respectively. “12” is the contribution to the phase plot from the overlapping of sources 1 and 2 that occurs in segment 2; this feature is of significant size compared with the features representing source 1 only and source 2 only. The Clusterwise PCA method attempts to group the plot into two “clusters.” The contribution 12 is clustered with 1, which degrades the estimate of this source. However, this contribution has less of an effect on the cluster corresponding to source 2, and hence this source is estimated much better, as observed in Table 2.Figure 6
Phase plot for segment 2. “1” and “2” refer to intervals in the phase plot corresponding to sources 1 and 2, respectively, and the label “12” refers to coincidence between the two sources.The robustness to the noise of the proposed method is now looked for segment 1 where there are no coincidences between the two underlying sources. The amplitude of source 2 for each set of data is equal to 1; the Gaussian noise is added with standard deviations varying from 0 to 30% of the amplitude of source 2. 2000 Monte Carlo runs are used.The absolute value of the correlation coefficient between each source and closest estimate is averaged over all Monte Carlo runs. A comparison is made between the performances of the proposed method, the Fast-ICA and the Clusterwise PCA methods. In Figure7(a), the mean correlation is plotted as a function of the standard deviation of noise for source 1, and in Figure 7(b) the corresponding plot is shown for source 2. Note that the range of noise standard deviations used in the plots are different as the extraction of the smaller source 2 breaks down at smaller noise levels.Mean correlation coefficient as a function of noise standard deviation: (a) source 1 (b) and source 2. Full line: Fast-ICA, dashed line, Clusterwise PCA; dotted line: MHC (reverse).
(a)
(b)In Figure7(a), it can be seen that MHC and Fast-ICA have comparable performances when estimating source 1, with the latter method having slightly more robustness to noise. The Clusterwise PCA method is more sensitive to noise, particularly for noise standard deviation larger than 0.2; the problem here is that if the noise is too large, then points in phase space will be associated with the wrong cluster. All the methods give a comparable performance when estimating source 1 for noise standard deviations up to 0.05, corresponding to 5% of the amplitude of source 2 in each set of data.In Figure7(b), when estimating source 2, Fast-ICA and Clusterwise PCA have indistinguishable performances for the range of standard deviations of noise used (0 to 5% of the amplitude of source 2). MHC is more sensitive to noise than the other two methods for noise standard deviations greater than about 1% of the amplitude of source 2. The reason for this is that Fast-ICA and Clusterwise PCA process the data globally, whilst MHC does this locally, so is more susceptible to noise.Now it is of interest to ask whether, when using the MHC, the reverse iteration method is more susceptible to noise than if the forward estimation only is used. In Figure8, the mean correlation coefficient between estimate and source for the forward MHC with no whitening is compared with MHC with the reverse iteration for source 2. We can see that reverse iteration does not contribute significantly more errors to the estimation of source 2. Hence most of the errors occurring in the MHC reverse method are occurring in the forward iteration.Figure 8
Extraction of source 2: full line: MHC with reverse iteration, dashed line: MHC with forward iteration only.For this particular simulation, the mixing matrix is given byA1=(6.5131) and the phase space plot for the mixed signals is shown in Figure 6. The directions corresponding to sources 1 and 2 are very clear in the directions (6.5,3) and (1,1). However, a much stricter test for the reverse MHC method is when these two directions are much closer. Simulations have also been carried out for the following mixing matrices:
(44)A2=(3.2131),A3=(3.02131)
for segments of time 0.5 to 2.5 seconds (no coincidence of sources 1 and 2) and 7 s to 9 s where there is coincidence. For 0.5 to 2.5 seconds, it is found that all three methods tested have a similar performance as when using A1, namely that Clusterwise PCA and MHC reverse extracts the sources exactly, and Fast-ICA extracts the sources with correlation coefficients 0.999 for source 1 and 1 for source 2. A stricter test of the performance of these methods when using mixing matrices A2 and A3 is to look at the segment of data, where there is coincidence of sources between 7 s and 9 s.  The performances of these methods for this data segment are shown in Tables 3 and 4.Table 3
Correlation between estimate and source for mixture matrixA2.
Method
Source 1
Source 2
MHC (Reverse)
1
1
Clusterwise PCA
0.7970
0.9894
Fast-ICA
0.9480
0.9997Table 4
Correlation between estimate and source for mixture matrixA3.
Method
Source 1
Source 2
MHC (Reverse)
1
1
Clusterwise PCA
0.7133
0.9992
Fast-ICA
0.9480
0.9997The MHC (reverse) method still gives perfect extraction of sources to be  within a scaling constant and has the best performance of the three methods tested for this particular dataset.


## 3.3. Two Sources, Varying Sparsity

It is of interest to investigate the performance of the proposed algorithm when processing mixtures of sparse sources with different “sparsity.” At one extreme, one has sources that completely overlap and where there are no segments where one source exists on its own. At the other extreme, one could have sources that do not overlap. An intermediate situation is where there is partial overlap between sources.In the next set of simulations, two mixtures of two sources are shown, where each source is generated from a uniform random number generator with values between −1 and +1; each source consists of 100 samples.The two sources are mixed using the randomly selected mixing matrix:(45)A=(0.5745807850.166294188-0.668082703-0.933171682).
Next, the second source is shifted by L samples with respect to the first source and the mixing matrix (45) applied to the two sources to produce mixtures with 100 + L samples. For each value of L, the MHC algorithm with reverse iteration is applied; the following parameters are found empirically to work for this data: vth=0.1 (18), λth=10-5 (36), and the minimum allowed value of  Corr(x,y) in (37) is Corrmin=0.9, with P=3 in (37).In Figure9, the correlation coefficient between estimated sources and closest actual sources, for sources 1 and 2, is plotted as a function of L for the Fast ICA, Clusterwise PCA, and MHC reverse methods.Mean correlation coefficient as a function of shift,L, between sources: (a) source 1 and (b) source 2. Full line: Fast-ICA; dashed line, Clusterwise PCA; dotted line: MHC (reverse).
(a)
(b)The following observations can be made.(i)
Clusterwise PCA extracts both sources with correlation coefficients larger than 0.98 forL≥70.(ii)
Fast-ICA should be expected to work best forL=0 and L around 100 because in these two extremes the sources are uncorrelated. At some intermediate values of L,where there is some correlation, it is expected that Fast-ICA should not work so well, and this is observed in this example for L between 30 and 50 for source 1. However, the estimates for source 2 have correlation coefficients larger than 0.98 for all values of L.(iii)
Fast-ICA has the best overall performance for all values ofL. However, if L≥3, then MHC (reverse) has the best performance as this method requires only three samples to work; Clusterwise PCA requires larger values of L to work because it performs a global clustering of the data.(iv)
The Fast-ICA algorithm is based on iteratively updating weights applied to the data. The initial values of these weights are chosen randomly; hence, one can obtain different results by rerunning the Fast-ICA algorithm and it is found that reinitializing two or more times can lead to better results for some values ofL. When using Fast ICA, the symmetric approach is found to work better than deflation for these signals. Changing the nonlinearity from gauss to pow 3 and sometimes changing the step size in the iteration of weights can lead to better results; in Figure 9, the best estimates, changing the input parameters, are used to construct the performance curves for Fast ICA. So, for this mixture of sources, the estimation using Fast-ICA can be very sensitive to the approaches used and to input parameters; however, for MHC (reverse), the same set of parameters works for all values of  L≥3.


## 3.4. 4-Lead ECG Data

In order to compare the performances on actual data, of the Fast-ICA, MHC, and Clusterwise PCA methods, data taken from the Daisy Database [20] will be analyzed. These data consists of ECG signals taken from an expectant mother. The data consists of 8 leads, 1 to 5 being abdominal and 6 to 8 thoracic. There are 2500 samples of data; there is some uncertainty about the sampling frequency, also pointed out in [10], but it is probably 250 Hz.These data were analysed in [17], where the MHC with forward iteration was applied to the data. Successful extractions of the fetal and maternal signals were achieved.In this paper, the MHC with reverse iteration is also applied. 4 leads are chosen for analysis: leads 1, 2, 3, and 4—this is a more challenging problem in extracting the fetal signal than when using 8 leads. In [17], all the data samples are analyzed; in this paper, data used is between sample numbers 1451 to 1751 where there is an overlap between a maternal and fetal complex as shown in Figure 10(a) where the maternal and fetal R-waves are marked by “M” and “F,” respectively.Abdominal and thoracic ECG data [20]. (a) Lead 1, (b) extraction of fetal ECG using Fast-ICA, (c) extraction of fetal ECG using MHC (forward iteration, Gram-Schmidt’s initialization), (d) Extraction of fetal ECG using MHC (Reverse iteration), (e) Clusterwise PCA.
(a)
(b)
(c)
(d)
(e)A comparison is made between the Fast-ICA (deflation, Gaussian), MHC forward iteration (with the Gram-Schmidt initialization), MHC with reverse iteration with no whitening of data, and Clusterwise PCA. The input parameters for the MHC calculations werevth=0.5 in (18), Corrmin=0.999, P=5 in (37), and  λth=10-5 in (36).For each method applied, the output that is closest to the expected fetal signal is chosen. The results for the Fast-ICA, MHC (forward), MHC (reverse), and Clusterwise PCA methods are shown in Figures10(b), 10(c), 10(d), and 10(e), respectively.It can be seen that all four methods give comparable results for the extraction of the fetal sources.Comparing Figures10(c) and 10(d), it is observed that the MHC with forward iteration gives similar results for the fetal source as the MHC with reverse iteration. This is due to the fact that despite the maternal/fetal coincidence in the data, there is still low correlation between the maternal and fetal sources. Hence when whitening is carried out, the vectors representing the maternal and fetal sources in the resulting phase plot are almost perpendicular and, from (16), this means that the extracted sources in the forward iteration are approximately a scaled version of the actual sources. The MHC with reverse iteration does not whiten the data initially but is still able to extract the fetal source with comparable quality to the other three methods applied to this data. It is interesting to note that noise is not having a deleterious effect on the quality of the estimated fetal sources in Figures 10(b), 10(c), 10(d), and 10(e); for this particular set of data, the anti-aliasing filter has already filtered off a significant amount of noise.


## 4. Discussion and Conclusions

The MHC method looks locally for threeconsecutive points that are the most linear of all points in phase space; hence just one direction corresponding to one source is estimated at any given time. Conventional clustering methods estimate either the mixing or the separating matrices and some algorithms can be extended to the underdetermined case, which cannot be done at present with the MHC.  The MHC is a more straightforward method to implement with no iterative clustering techniques; however this technique requires the sources to be strictly sparse whilst clustering techniques can deal with sources that are approximately sparse where one looks for clustering of points in phase space around a straight line [11]. As the MHC analyzes data in phase space locally, it would be more sensitive to the effects of noise than the  clustering techniques. For the MHC to work, the underlying sources need to be sparse for only a few (down to three data points) which can be seen as an advantage of this method over methods that look for the clustering of data points along a line in phase space, where it is assumed that the underlying sources are sparse over several segments of the mixed data.The MHC reverse iteration method can be sensitive to the choice of input parametersvth in (18), λth in (36), Corrmin, and P in (37). On the other hand, both the Fast-ICA and Clusterwise PCA algorithms have been found to be relatively insensitive to the choice of the input parameters. This sensitivity to input parameters implies that, in its present form, the reverse MHC method should be used to analyse signals offline, rather than online, where the parameters can be varied to obtain optimal estimation of the underlying sources.It should be noted that MHC relies onall sources being sparse. The more the number of sources, the less likely this assumption will be obeyed. In addition, as discussed in Section 3.1, the sensitivity to the noise of the method limits the number of sources that could be extracted. One potential application of the MHC method is the analysis of electron paramagnetic resonance (EPR) data [5] to detect the presence of free radicals in chemical and biological samples where there are two mixtures of two sources to be analyzed. In [5] Fast-ICA is applied to two sets of measured spectra to extract the spectra corresponding to the individual radicals. The results found by authors are promising. However, the underlying source spectra (e.g., for OH and O2) consist of overlapping peaks, and hence there is some correlation between the underlying sources. In addition, the underlying spectra can be modeled as sparse as there are distinct peaks in the data. It would be of interest to investigate whether the MHC method could also be used in this application and whether it would give any improved results compared with using Fast-ICA.There are BSS methods existing in the literature that derive methods to estimate the underlying sources where the underlying sources are assumed to be correlated. For example, in [21], a method based on using precoders is described; the algorithm to estimate the sources is more complicated than the method proposed in this paper, but the former method does not assume that the underlying sources are sparse. In [22] it is assumed that the underlying sources are nonnegative, and, as for [21], the sources are not assumed to be sparse. It would be of interest to compare the method proposed in this paper with the algorithms in [21, 22] for mixtures of sparse, correlated sources.



---

*Source: 218651-2013-02-13.xml*