Patent Grant
7519512
(1) Field of the Invention
(2) Description of the Art
This invention relates to dynamic blind signal separation, that is to say blind signal separation which can be used inter alia in a non-stationary environment, and to a method, an apparatus and a computer program for implementing it.
In recent years, blind signal (or source) separation (BSS) has emerged as an important field in signal processing. It is now successfully exploited in applications such as biomedicine, financial analysis, speech recognition and telecommunications.
The objective of BSS is to recover or separate a number of statistically independent signals given only observations of mixtures of those sources detected using sensors, together with noise and unwanted interference. The term ‘blind’ is used to indicate that no prior information (e.g. waveform, frequency, bandwidth, location) is available concerning the individual sources or the manner in which they have been combined.
Although the term ‘blind’ is often associated with these techniques there are several assumptions that need to be satisfied such that signal separation can be achieved. The mixing conditions, i.e. the combination process of the signals at the sensors, must be linear and time-invariant. The term time-invariant means that the spatial statistics (combinations) of the signals received at the sensor array do not change with time. For this invention, another assumption for the combination process is that the signals arrive at the same time at the sensors. This is referred to as the instantaneous mixing case.
Many BSS methods have been developed for the problem when the combination (mixing) process is linear, instantaneous and time-invariant. These include:
Published International Application No. WO 03/028550 A2 discloses a specific application of BLISS to monitoring fetal ECG signals using abdominal electrodes applied externally to a mother. Electrode signals are filtered and subjected to a BSS technique based on Independent Component Analysis (ICA), I. J. Clarke, “Direct Exploitation of non-Gaussianity as a Discriminant”, EUSIPCO '98, Rhodes, Greece, 8-11 Sep., 1998.
Published International Patent Application No. WO 01/17109 discloses BSS by windowing data, applying a discrete Fourier Transform and cross-correlating the result. A gradient descent method is then used to define a finite impulse response filter which will separate mixed signals. U.S. Pat. No. 6,321,200 (Casey) discloses extraction of features from mixed signals by filterbank analysis, windowing resulting data to produce multidimensional observation matrices, reducing matrix dimensionality and processing by independent component analysis. Published International Patent Application No. WO 03/090127 discloses BSS using fourth-order cumulants. BSS of disjoint orthogonal signals is disclosed by A. Jourjine et al, Proc. 2000 IEEE Conference on Acoustics, Speech and Signal Processing, 2000 ICASSP'00, 5-9 Jun. 2000, Istanbul Turkey, Vol. 5, pp 2985-2988.
When non-stationary conditions apply, the accuracy with which some prior art BSS methods can detect and separate signals is restricted because it is normally not an objective of the prior art to deal with non-stationarity. This is because only a limited amount of data can be used reliably in the separation process, and statistically based estimates of the parameters may not be accurate. Additionally, processing of large data sets (e.g. over hours) is not possible. These limitations have led to the development of BSS methods for the non-stationary problem. Such methods include:
Particle filtering is a relatively new area of signal processing and it provides a robust method for tracking signals. However, its computational cost renders its use impractical for real time BSS.
In contrast to this, the EASI and NGA (with a Kalman filter) algorithms are much more computationally efficient. These techniques rely on stochastic and natural gradient approaches respectively. They tend to have slow convergence and in a non-stationary environment they are only suitable for tracking a slowly changing system.
Alternatively, a more robust technique uses the ‘block’ algorithms (such as BLISS and JADE described earlier) but in a moving window approach. The term window is used for a section or ‘block’ of data. It is assumed that the data window advances in time to cover all recorded data as a series of sections and the signals are separated at each window location. These block algorithms process blocks or windows of data in isolation from one another. They may be used in blind demixing, and they assume that:
A major problem with using a moving window approach is called ‘signal swapping’. This arises from unmixed signals being in different orders in successive windows without the orders being discernible: i.e. a first window might be processed to separate six signals numbered 1 to 6 in the order 123456, and processing of the immediately following window might yield these signals in a different unknown order such as 413562. The resulting observed signals therefore contain parts assembled perhaps randomly from different source signals. It may not be convenient, practical or possible to reorder the signals correctly to allow continuous signal monitoring.
The blind demixing algorithm overcomes this problem by assuming prior knowledge of signal statistics to reorder signals. However in many applications signal statistics change from window to window and they may also be unknown.
Moreover, as in particle filtering, the computational cost of a moving window technique can render its use impractical for real time BSS processing.
It is an object of the invention to provide an alternative method of signal separation.
The present invention provides a method for dynamic blind signal separation including processing signals associated with windows of data characterised in that the method includes:
The invention provides the advantage that it is normally capable of processing data faster than the prior art of block algorithms, because initialisation avoids the need to carry out a full analysis of each data window (other than the first such) and it results in computation being reduced. This is important in applications such as fetal heartbeat monitoring in particular, where processing in as near to real time as possible is desirable to enable clinical intervention in the case of an emergency. Moreover, as will be described later, separate initialisation of orthogonality and independence avoids a problem associated with two stage BSS, in which signals become normalised and consequently cannot be updated appropriately. The invention is applicable to monitoring systems changing sufficiently slowly with time such that leading window results are useable in processing of following windows.
Processing of the approximate results may incorporate updating orthogonality using small updates to produce decorrelation in a second order statistics procedure. Updating orthogonality may be implemented by a technique referred to as Jacobi and involving diagonalisation of a symmetric matrix by determining and applying rotations iteratively until off-diagonal elements of the matrix become substantially equal to zero.
The method may include a second stage of initialisation using results obtained for each leading window to initialise independence of decorrelated signals associated with the respective following window, this second stage using independent component analysis (ICA) to apply small rotation updates to initialised signals in a higher than second order statistics procedure to produce signal independence and separation. The higher than second order statistics procedure may be at least one of a third order and a fourth order statistics procedure.
The method may be implemented in an acquisition phase in which signals are separated and desired signals are identified among the separated signals, and in a subsequent phase in which only desired signals are processed to separation.
The signals associated with windows may be statistical measures of data in the windows.
The method may comprise an acquisition stage of processing a first leading window of data to obtain first results, and a subsequent stage of processing following windows by iteratively updating immediately preceding results using subsequent data snapshots to produce snapshot results and combining the snapshot results with the immediately preceding results, the immediately preceding results being those obtained in a respective immediately preceding update if any and being the first results otherwise. Prior to combining the snapshot results with the immediately preceding results, the immediately preceding results may be weighted with a forget factor to implement an exponentially fading following window.
The first results may comprise a mean vector of signal samples and a covariance matrix of a data matrix of the first leading window in each case, and decorrelation and normalisation of the data matrix to obtain signal vectors from which to obtain their moment as a fourth order tensor.
The snapshot results may comprise a mean snapshot vector and a snapshot covariance matrix, a decorrelated and normalised snapshot equivalent providing signal vectors from which to obtain their moment as a fourth order tensor update, and prior to combining the snapshot results with the immediately preceding results, the snapshot results may be weighted by a forget factor p and the immediately preceding results may be weighted by a further forget factor (1−p) to produce an exponentially fading window, where 0<p<1.
In another aspect, the present invention provides computer apparatus for dynamic blind signal separation programmed to process signals associated with windows of data characterised in that the computer apparatus is also programmed to:
The computer apparatus may be programmed to update orthogonality of the approximate results using small updates to produce decorrelation in a second order statistics procedure. It may be programmed to update orthogonality by a technique referred to as Jacobi and involving diagonalisation of a symmetric matrix by determining and applying rotations iteratively until off-diagonal elements of the matrix become substantially equal to zero.
The computer apparatus may be programmed to implement a second stage of initialisation using results obtained for each leading window to initialise independence of decorrelated data associated with the respective following window, this second stage involving ICA to apply small rotation updates to Initialised data in a higher than second order statistics procedure to produce signal independence and separation. The higher than second order statistics procedure may be at least one of a third order and a fourth order statistics procedure.
The computer apparatus may be programmed to implement an acquisition phase in which signals are separated and desired signals are identified among the separated signals, and a subsequent phase in which only desired signals are processed to separation.
The signals associated with windows may be statistical measures of data in the windows.
The computer apparatus may be programmed to implement an acquisition stage of processing a first leading window of data to obtain first results, and a subsequent stage of processing following windows by iteratively updating immediately preceding results using subsequent data snapshots to produce snapshot results and combining the snapshot results with the immediately preceding results, the immediately preceding results being those obtained in a respective immediately preceding update if any and being the first results otherwise. It may be programmed to implement an exponentially fading following window by weighting the immediately preceding results with a forget factor prior to combining the snapshot results with the immediately preceding results. The first results may comprise a mean vector of signal samples and a covadance matrix of a data matrix of the first leading window in each case, and decorrelation and normalisation of the data matrix to obtain signal vectors from which to obtain their moment as a fourth order tensor.
The snapshot results may comprise a mean snapshot vector and a snapshot covariance matrix, and the computer apparatus may be programmed to produce a decorrelated and normalised snapshot equivalent and to obtain therefrom signal vectors and their moment as a fourth order tensor update, to weight the snapshot results by a forget factor p and the immediately preceding results by a further forget factor (1−p) to implement an exponentially fading window, where 0<p<1, and to implement such weighting prior to combining the snapshot results with the immediately preceding results.
In a further aspect, the invention provides computer software for dynamic blind signal separation including processing signals associated with windows of data characterised in that the software has instructions for controlling computer apparatus to:
The computer software may include instructions for processing the approximate results using small updates to update orthogonality and produce decorrelation in a second order statistics procedure. It may include instructions for updating orthogonality by a technique referred to as Jacobi and involving diagonalisation of a symmetric matrix by determining and applying rotations iteratively until off-diagonal elements of the matrix become substantially equal to zero.
The computer software may include instructions for implementing a second stage of initialisation using results obtained for each leading window to initialise independence of decorrelated data associated with the respective following window, this second stage using ICA to apply small rotation updates to initialised data in a higher than second order statistics procedure to produce signal independence and separation. The higher than second order statistics procedure may be at least one of a third order and a fourth order statistics procedure.
The computer software may include instructions for implementing an acquisition phase in which signals are separated and desired signals are identified among the separated signals, and a subsequent phase in which only desired signals are processed to separation.
The signals associated with windows may be statistical measures of data in the windows.
The computer software may include instructions for implementing an acquisition stage of processing a first leading window of data to obtain first results, and a subsequent stage of processing following windows by iteratively updating immediately preceding results using subsequent data snapshots to produce snapshot results and combining the snapshot results with the immediately preceding results, the immediately preceding results being those obtained in a respective immediately preceding update if any and being the first results otherwise. It may include instructions for implementing an exponentially fading following window by weighting the immediately preceding results with a forget factor prior to combining the snapshot results with the immediately preceding results.
The first results may comprise a mean vector of-signal samples and a covariance matrix of a data matrix of the first leading window in each case, and decorrelation and normalisation of the data matrix to obtain signal vectors from which to obtain their moment as a fourth order tensor. The snapshot results may comprise a mean snapshot vector and a snapshot covariance matrix, and the computer software may include instructions for producing a decorrelated and normalised snapshot equivalent and for obtaining therefrom signal vectors and their moment as a fourth order tensor update, and for weighting the snapshot results by a forget factor p and the immediately preceding results by a further forget factor (1−p) to implement an exponentially fading window, where 0<p<1, such weighting being prior to combining the snapshot results with the immediately preceding results.
In order that the invention might be more fully understood, an embodiment thereof will now be described, with reference to the accompanying diagrams, in which:
Referring to
The sensors 14 supply input signals to a data acquisition stage 16 of the system 10 in which the Input signals are sampled in sets at a rate of approximately 500 sets of samples/second. The samples are then digitised and recorded. Each set of samples consists of a respective sample from each sensor 14, and the samples in each set are recorded substantially simultaneously. After recordal, these signals pass to a signal separation stage 17 where they are processed to yield separated or unmixed signals corresponding to source signals 15. The separated signals pass to an analysis stage 18, where they are analysed to distinguish fetal ECG signal, maternal ECG signal and noise.
The structure of signals recorded in fetal ECG monitoring are complex, as the required fetal ECG signal is subject to interference signals such as the much stronger maternal ECG signal and muscle noise. Typical prior art BSS methods are restricted to a relatively short section or time window of data over which prior art BSS assumptions hold good. There is however a well known and long felt want for signals to be separated and monitored over relatively longer periods of time, for which comparable prior art methods are not adequate.
Data from the acquisition stage 16 is processed as a series of what are referred to as data “windows” consisting of digital signals derived by digitising signals from each sensor 14, a typical signal may consist of 2,000-3,000 signal samples. Each window may partially overlap its preceding window, i.e. sets of samples may be common to both as described later. It is assumed that the data window advances through a succession of window locations in a block of data, signals are separated for each window location and eventually the entire data block is processed. A first data window is processed in the signal separation stage 17 by a block mode method such as BLISS or JADE. Before describing this in more detail, a signal model for each window will be established and BSS used therewith.
In fetal ECG analysis, the number of signal samples per set is normally taken to be equal to the number of signal electrodes, i.e. ignoring ground or reference electrodes.
For BSS, an appropriate (noise-free) data model for data window number k having m digital signals from respective signal electrodes, each signal with n samples, is given by:
Xk=AkSk (1)
where: Xk is an (m×n) matrix of data in window number k, and has rows which are respective signals (each expressed as a set of m data samples); Sk is an (m×n) data matrix of similar kind to Xk but representing unmixed signals, i.e. signals which have been separated and corresponding to independent source signals 15; and Ak is an (m×m) mixing matrix which represents a transformation by which the source signals 15 are combined to form mixed signals received by electrodes and corresponding to Xk. The mixing model in each window expressed by Equation (1) is assumed to be linear, instantaneous and time-invariant. Data in each window is assumed to have a zero mean (or indeed the data is transformed to zero mean by subtracting its mean from each data sample), as this is convenient and does not affect the required result.
The signal separation problem to be solved in BSS is to recover or unmix signals (up to a permutation) represented by the data matrix Sk: this can be done by determining an inverse mixing matrix Ak−1, the inverse of Ak, and premultiplying both sides of Equation (1) by it, i.e.:
Ak−1Xk=Ak−1AkSk=Sk (2)
Equation (2) for the kth data window can be written for the immediately following or (k+1)th window as:
Ak+1−1Xk+1=Sk+1 (3)
To determine the inverse mixing matrix Ak−1, a BSS algorithm may use two stages. The first stage, which only uses second order statistics, is to decorrelate and normalise the signals. The second stage, which normally relies on fourth order statistics, makes the decorrelated and normalised signals (obtained from the first stage) independent.
In order for these two processes to be carried out efficiency and effectively in each window, the tracking method has several important features. The first two features include initialisation and update processes.
The expression “initialisation” means that information about a mixing matrix for an immediately preceding data window is used to initialise signals from a current window. This initialisation has the effect of rendering initialised signals nearly separated, so that a complete computation to unmix signals is not required. Although in this example the signals themselves are initialised, it is not in fact essential to do so: it is also possible to initialise intermediate quantities instead, i.e. parameters derived from data instead, e.g. the data's statistical measures such as its covariance matrix (as will be described later). To determine the mixing matrix of the current window, initialised signals will only require updating by a smaller amount than would be the case without initialisation. Moreover, as will be described later in more detail, initialisation in combination with small updates can prevent the problem of signal swapping.
It has been considered that the signals in the current window could be initialised using Ak−1, such that a first approximation S′k+1 to Sk+1 is provided by multiplying Xk+1 by Ak−1 instead of Ak+1−1, i.e.:
Ak−1Xk+1=S′k+1 (4)
The estimated S′k+1 provided by Equation (4) will be closer to the desired value of Sk+1 than Xk+1. However, in the two stage BSS approach (mentioned above), estimated signals are normalised to have the same power. This has a consequence that, in a first stage (decorrelation), signals cannot be updated. This will be shown later. Even if the normalisation process were to be ‘undone’, updating the decorrelation process would then undo the effect of the previous stage. A feature of the example of the invention to be described is that it avoids this problem by updating the corresponding second and higher order components of the demixing matrix separately.
In the present example, an initialisation process is implemented in two stages as will be described in more detail later: in the first of these stages, the orthogonality of the signals is initialised and then the signals' orthogonality is updated using a decorrelation method.
This makes the signals orthogonal. Then, in the second stage, the independence of the signals is initialised and the signals' independence is updated using a higher order statistics stage. This makes the signals independent.
In the first of the stages referred to above, one option is to carry out a decorrelation and normalising method on the (m×n) matrix Xk. This may be referred to as a singular value decomposition (SVD) when the singular values are arranged in decreasing magnitude order (M. Moonen and B. De Moor, “SVD and Signal Processing, III Algorithms, Architectures and Applications”, ASIN: 0444821074, Elsevier 1995). However, ordering of signals in this manner gives difficulty, because signal powers can change from window to window resulting in signal reordering. Hence signal swapping can occur, i.e. the process generates apparent signals which are actually segments of different signals joined together. The terminology of decorrelation and normalisation is preferred here to SVD, as signals are not ordered according to their power.
A second order decomposition of the data matrix Xk can be expressed as:
Xk=UkΣkVkT (5)
where Uk is an (m×m) unitary matrix containing eigenvectors of the data matrix's covariance matrix XkXkT; Vk is an (n×n) unitary matrix containing eigenvectors of XkTXk; superscript index T in XkT or VkT denotes a transpose of matrix Xk or Vk: and Σk denotes an (m×n) diagonal matrix having singular values (square roots of eigenvalues) on its diagonal and zero values in all off-diagonal positions. These singular values are not necessarily arranged in order of decreasing magnitude and they indicate the relative contributions of associated eigenvectors to data to which both correspond.
The columns of Uk are respective orthonormal singular vectors: they are spatial vectors indicating mixing of the source signals 15. The rows of VkT are orthonormal singular vectors which are estimates of the source signals 15, but in general these estimates are not fully separated versions of source signals because they require a further relative rotation to become so. Decorrelation fails to separate the source signals fully because it is a second order decorrelation by which decomposed vectors are made mutually orthogonal, i.e. uncorrelated with one another. In many applications of the invention, the spatial characteristics (locations) of the signals will be similar, and so a solution that makes them dissimilar will not succeed in separating them.
This failure to separate signals is illustrated in
Plot (b) is a scatter plot of two signals which are mixtures of the source signals, i.e. after mixing but before processing to unmix them. It shows that the mixing process corresponds to a complicated combination of stretching, shearing and rotation. Plot (c) is a scatter plot of two signals derived by decorrelation and normalisation of the plot (b) signals. Such decorrelation and normalisation of signals is obtainable by application of a processing method which in statistical terms is restricted to second order: this method provides estimated signals from which stretching and shearing have been removed. However, comparison of plots (a) and (c) shows that estimated signals in (c) are rotated relative to (a): i.e. the second order method has failed to remove a rotation. The estimated signals are therefore related to the source signals 15 by an unknown rotation matrix, and it is necessary to determine this matrix to rotate the estimated signals in plot (c) to obtain appropriately separated signals. Plot (d) shows two signals obtained by rotating those in plot (c). Processing which produces only a rotation does not produce shearing or stretching, and thus does not counteract results achieved by the second order method. The rotation is chosen to align edges of scatter plot (d) with co-ordinate axes (not shown). The correct rotation is determined as will now be described using higher order statistics, where “higher order” means higher than second order.
In signal separation, determination of the unknown rotation matrix may satisfy a statistical independence condition. Mathematically, statistical independence of m signals X1 to Xm can be expressed as an ability to factorise the signals' joint probability function p(x1 . . . xm) into a product of the signals' joint probability functions p(xi), i.e.:
This process uses higher order statistics (HOS). The use of HOS to separate unknown and independent signals is often referred to as independent component analysis (ICA, which normally involves second and higher order statistics), and this is the second stage in the separation process in this example.
From Equation (5), and only using data with non-zero singular values, Xk can be expressed as Ũk{tilde over (Σ)}k{tilde over (V)}kT, where Ũk is an (m×m) matrix consisting of the first m columns of Uk; {tilde over (V)}kT is an (m×n) matrix consisting of the first m rows of VkT; and {tilde over (Σ)}k is an (m×m) matrix containing the singular values of the kth window data matrix Xk. From
Xk=Ũk{tilde over (Σ)}kRkTRk{tilde over (V)}kT (7)
In Equation (7), the product Ũk{tilde over (Σ)}kRkT is an estimate of the mixing matrix Ak and the product Rk{tilde over (V)}kT is a matrix having rows which are respective estimated signals. The estimated signals in Rk{tilde over (V)}kT are orthonormal, i.e. orthogonal and normalised to have the same power, but have relative powers stored in the estimated mixing matrix Ak. This normalisation aspect is inherent in the second order stage approach to the blind signal separation process.
Rk is calculated (as will be described later) by using an iterative sequence of pairwise rotations, each of which is designed to maximise the statistical independence of a respective pair of signals (taken from {tilde over (V)}kT) as measured by a corresponding pairwise cost function. In this maximisation process, the BLISS and JADE algorithms use fourth order statistics for the higher order statistics stage. In the present embodiment, the BLISS higher order statistics stage is used.
For each window, the main objective is to calculate parameters contained in Equation (7). If each window were to be processed in isolation as in the prior art, then results for each window would be determined separately: the computational cost of using k windows would then be k times that for a single window. This is impractical if the number of windows is large, e.g. >20 windows with 2,500 samples per signal per window. However, as previously described it has been found that it is possible to reduce the processing required by using results obtained for each window in initialisation of the separation process for a subsequent window. This initialisation converts data into an intermediate form which is closer to the desired result and reduces processing. It results in signals which are closer to being separated and a full computation of signal statistics is not required. The mixing matrix Ak need not be stationary to permit this, but should be changing relatively slowly with time for the purposes of this invention.
Having discussed the BSS problem, it can now be shown that initialising signals with Ak is not feasible. For example, rewriting Equation (7) by substituting index (k+1) for k, multiplying both sides by (Ũk{tilde over (Σ)}kRkT)−1, the inverse of the estimated mixing matrix Ak for the kth window, and replacing inverses of matrices by transposes:
(Ũk{tilde over (Σ)}kRkT)−1Xk+1=(Rk{tilde over (Σ)}k−1ŨkT)Xk+1=(Rk{tilde over (Σ)}k−1ŨkT)(Ũk+1{tilde over (Σ)}k+1Rk+1TRk+1{tilde over (V)}k+1T) (8)
In Equation (8), (Ũk+1{tilde over (Σ)}k+1Rk+1TRk+1{tilde over (V)}k+1T) is unknown and has to be estimated from the data. Ũk and Rk are unitary matrices, and therefore Ũk−1=ŨkT, ŨkTŨk=I and ŨkŨkT=I, where I denotes the identity matrix. If Ak is slowly changing, then ŨkTŨk+1≈I, RkRk+1T≈I and {tilde over (Σ)}k−1{tilde over (Σ)}k+1≈I. Using these three approximations:
(Ũk{tilde over (Σ)}kRkT)−1Xk+1≈Rk+1{tilde over (V)}k+1T (9)
On the right of Approximation (9), the product Rk+1{tilde over (V)}k+1T is a matrix having rows which are respective estimated signals: these signals are initialised and close to being orthonormal, i.e. to orthogonality and equality of signal powers, and can be processed by moving window ICA using a decorrelation method. However, they were obtained using a decorrelation and normalisation method in Equation (5), and this discriminates between unnormalised signals on the basis of signal power. Such discrimination makes Equation (5)'s decorrelation process much more difficult. A data window initialisation process is preferred (as will be described) which initialises the orthogonality of the incoming signals, i.e. those in the (k+1)th data window, but which avoids normalisation of these signals.
The present example is implemented in two separate stages, each involving a respective initialisation process: in the first of these stages, the orthogonality of the signals is initialised by applying the kth window matrix Ũk to them, and then their orthogonality is updated to reflect data in the (k+1)th window using a decorrelation method. This makes the signals orthogonal. Then, in the second stage, the independence of the signals is initialised by applying the kth window rotation matrix Rk to them, and then their independence is updated using higher order statistics. This makes the signals independent.
The matrix Rk is input at 24 and applied to the orthonormal signals at 25 to initialise them using Equation (26). The initialised orthonormal signals then undergo separation by rotation at 27 by ICA with small angle updates using higher order statistics and Equation (29). Here again, “small” means more than one possible updating angle is available and the smaller/smallest angle is chosen. Steps 25 and 27 are collectively a higher order stage of processing in statistical terms producing separated signals at an output 28.
In many decorrelation routines, signals are made orthonormal by using an iterative sequence of pairwise rotations. Each of these rotations is designed to maximise the orthogonality of a given pair of signals as measured by a corresponding pairwise (second order) cost function. In the Jacobi method, this can be accomplished by reducing a symmetric matrix to diagonal form using Given's rotation as described in the Haykin reference mentioned above. For example, in a two signal case the objective of using the Jacobi method is to diagonalise a (2×2) covariance matrix X12 of two signals x1 and x2 expressed as (n·x1) vectors. This matrix is given by Haykin as:
This can be implemented by iteratively forming a (2×2) matrix B by rotating the Equation (10) covariance matrix using a rotation matrix Q to force off-diagonal elements of B to zero, i.e.:
At each iteration of Given's rotation, the rotation matrix Q is determined which minimises off-diagonal elements of the symmetric matrix X12. When the matrix X12 becomes fully diagonal, i.e. when all its off-diagonal matrix elements have become zero, the signals it represents have become mutually orthogonal. This can be achieved by simply setting b12 to zero in Equation (11), and determining a value of Θ that satisfies this. B is a symmetric matrix so b12≈b21 and it is only necessary to determine Θ for b12=0. By forcing b12 and b21 to be equal to 0, signals x1 and x2 are made orthogonal and hence B becomes a diagonal matrix. Thus, the diagonalisation of B can be a necessary condition for making vectors orthogonal.
Using Equation (7), a covariance matrix XkXkT is formed for Xk, i.e.:
XkXkT=(Ũk{tilde over (Σ)}kRkTRk{tilde over (V)}kT){tilde over (V)}kRkTRk{tilde over (Σ)}kTŨkT (12)
Because Rk is a unitary matrix, and because estimates of signals contained in the rows of {tilde over (V)}k are orthogonal to one another:
{tilde over (V)}kT{tilde over (V)}k=RkTRk=I (13)
Then XkXkT=˜k{tilde over (Σ)}k{tilde over (Σ)}kT˜kT=˜k{tilde over (α)}k˜kT (14)
where {tilde over (α)}k is an (m×m) diagonal matrix whose diagonal elements are eigenvalues of XkT.
Writing (k+1) for k in Equation (14), for the (k+1)th window, the covariance matrix Xk+1Xk+1T is expressed as:
Xk+1Xk+1T=Ũk+1{tilde over (α)}k+1Ũk+1T (15)
In order to avoid computing decorrelation afresh, the covariance matrix Ũk+1{tilde over (α)}k+1Ũk+1T is initialised by premultiplying by ŨkT and postmultiplying by Ũk, which were obtained for the immediately preceding or kth window, i.e.:
ŨkT(Xk+1Xk+1T)Ũk=ŨkT(Ũk+1{tilde over (α)}k+1Ũk+1T)Ũk (16)
Equation (16) represents Ũk initialising orthogonality for information or signals or vectors in the covariance matrix Xk+1Xk+1T of the data matrix Xk+1 using results obtained for the kth data window. Although signals represented in the covariance matrix Xk+1Xk+1T are not expressly those of the data matrix Xk+1, they are derived from those of the data matrix Xk+1. The expression “signals associated with data windows” and similar expressions as used herein therefore mean signals of the data matrix Xk+1 or signals derived from them such as those of the covariance matrix Xk+1Xk+1T.
The covariance matrix Xk+1Xk+1T is a second order measure of the statistics of the data matrix Xk+1. Its initialisation is at 22 in
As previously mentioned, the covariance matrix being diagonal is a condition for conferring orthogonality upon vectors having covariance represented by that matrix. Initialisation at 22 is followed by a further improvement in orthogonality brought about using information from the (k+1)th data window as described below.
Now in a slowly changing system, ŨkTŨk+1≈I and Ũk+1TŨk≈I, and so Equation (16) can be rewritten as:
ŨkT(Ũk+1{tilde over (α)}k+1Ũk+1)Ũk≈{tilde over (α)}k+1 (17)
For an (m×m) symmetric matrix, and without initialisation of data, a Jacobi diagonalisation method is iterative with of order m3 operations required for each sweep consisting of a single update of all the relevant signals represented as rows of the matrix Xk+1Xk+1T. Many sweeps may be required. However, the matrix (ŨkTXk+1Xk+1TŨk), initialised at 22 in accordance with the invention and given by Equations (16) and (17), is already close to representing a diagonal matrix: the Jacobi diagonalisation method therefore requires only a few iterations to converge. Thus, the computational cost is reduced in accordance with the invention by using an transformation to initialise signal orthogonality in the (k+1)th window, i.e. the current window.
The next step is to complete diagonalisation for the (k+1)th window using in this example the Jacobi method at 23. It employs matrices ŨzT and Ũz which are eigenvectors estimated when applying Jacobi to the matrix (ŨkTXk+1Xk+1TŨk) which has been initialised as regards orthogonality. ŨzT and Ũz are rotation matrices implementing small rotations that diagonalise the initialised matrix (ŨkTXk+1Xk+1TŨk), i.e. they should make vectors in this matrix orthogonal by updating as follows, i.e.:
ŨzT(ŨkTXk+1Xk+1TŨk)Ũz={tilde over (α)}k+1 (18)
In Equation (18), Ũz is computed using Jacobi diagonalisation of the orthogonality-initialised covariance matrix in Equation (16): a unitary matrix is required to diagonalise this matrix. As before, a transform in the form of a rotation implemented by a matrix Q is sought such that QT(Xk+1Xk+1T)Q is diagonal, where Q=ŨkŨz. This uses the fact that Ũk initialises data and Ũz is based on small updates. They can be used in the form of their product as they are both are unitary and their product is then also unitary. Estimates of eigenvectors Ũk+1 for the (k+1)th window are expressed by:
Ũk+1=ŨkŨz (19)
and the estimated signals ({circumflex over (V)}k+1T) are defined in Equation (24). This all occurs in stages 22 and 23.
For the higher order stage 25/27, signals obtained from Equation (24) are initialised in Equation (26) using Rk. This occurs in stage 25, shown in
Writing k+1 for k in Equation (5):
Xk+1=Ũk+1{tilde over (Σ)}k+1{tilde over (V)}k+1T (20)
Premultiplying Equation (20) by Ũk+1T:
Ũk+1TXk+1=Ũk+1TŨk+1{tilde over (Σ)}k+1{tilde over (V)}k+1T={tilde over (Σ)}k+1{tilde over (V)}k+1T (21)
Premultiplying Equation (21) by {tilde over (Σ)}k+1−1:
{tilde over (Σ)}k+1TŨk+1TXk+1={tilde over (Σ)}k+1−1{tilde over (Σ)}k+1{tilde over (V)}k+1T={tilde over (V)}k+1T (22)
where {tilde over (Σ)}k+1 is a diagonal matrix having elements which are the square roots of the eigenvalues contained in the matrix {tilde over (α)}k+1, i.e.
{tilde over (Σ)}k+1=√{square root over ({tilde over (α)}k+1)} (23)
Orthonormal signal vectors contained in the matrix {tilde over (V)}k+1T (for use in the higher order statistics stage 25/27 of the method 20) can be expressed by eliminating {tilde over (Σ)}k+1 from Equation (22) using Equation (23), i.e.:
Vk+1T={tilde over (α)}k+1−1/2Ũk+1TXk+1 (24)
The vectors or signals {tilde over (V)}k+1T are decorrelated signals output from the second order stage 22/23 in
For the (k+1)th window, m orthonormal signal vectors each with non-zero eigenvalues and zero mean have been obtained as {tilde over (V)}k+1T from the Jacobi decorrelation stage 23. They are now expressed as a product of an (m×m) unitary rotation matrix Rk+1T and signal estimates Ŝk+1 which are unmixed versions of the source signals 15.
{tilde over (V)}k+1T=Rk+1TŜk+1 (25)
The estimates Ŝk+1 are related to the source signals 15 by an arbitrary permutation to reflect possible signal swapping and a scale factor. It is necessary to estimate a rotation matrix Rk+1 which is the transpose and inverse of Rk+1T in order to multiply both sides of Equation (25) by it: this eliminates Rk+1T from the right hand side of Equation (25) and consequently recovers the required estimates Ŝk+1. If the mixing matrix is changing slowly, then by using information from the preceding or kth window:
Rk{tilde over (V)}k+1T=RkRk+1TŜk+1 (26)
Equation (26) expresses decorrelated signals output from step 23 having their independence initialised at 25 by multiplication by Rk obtained from the kth data window. This initialises fourth order information in the decorrelated signals. To improve signal independence and separate the signals, they are updated in step 27 on the basis of the (k+1)th data window as follows:
RkRk+1T≈I for slow change, and so:
Rk{tilde over (V)}k+1T≈Ŝk+1 (27)
The higher order statistics stage 25/27 in a block based algorithm such as BLISS is implemented in accordance with the invention using the initialised matrix Rk{tilde over (V)}k+1T of signals from Equation (27). This stage is iterative with approximately the following number of operations per sweep, a sweep being a single update of all signal pairs in the matrix Rk{tilde over (V)}k+1T.
No. of operations per sweep ≈20 m(m−1)n/2 (28)
where the term m(m−1)/2 denotes the number of pairwise operations and n denotes the number of signal samples in the (k+1)th window or elements in the associated data matrix Xk+1. Since the signals defined in Equation (27) are close to being separated, only a few iterations of sweeps in the higher order statistics stage 25/27 will be required for convergence of the signals to separation. Separation is judged to have occurred when angles obtained from the BLISS estimation process are small. Initialisation in accordance with the invention reduces the number of computation steps required to unmix signals. Signal independence is improved by updating in step 25, and unmixed or separated estimates Ŝk+1 of original signal sources are given by:
Ŝk+1=Rz(Rk{tilde over (V)}k+1T) (29)
where Rz is a unitary rotation matrix implementing a small rotation or update. Rz is computed using the initialised vectors (Rk{tilde over (V)}k+1T) in the higher order stage 25/27. Similarly to Rk, Rz is calculated by using an iterative sequence of pairwise rotations, each of which is designed to maximise the statistical independence of a respective pair of signals taken from {tilde over (V)}k+1T, independence being measured by a corresponding pairwise cost function. A unitary rotation matrix Rk+1 for the (k+1)th window is required to rotate signal vectors obtained by decorrelation as described earlier. This uses the fact that Rk initialises data and Rz is a small update to Rk: they can be expressed as a product as they are both are unitary and their product is also unitary. The rotation matrix Rk+1 is therefore expressed as a product of an initialisation matrix and an update matrix, i.e.:
Rk+1=RzRk (30)
The initialisation arid update processes described with reference to Equations (19) to (30) can reduce the computation load associated with separating signals. Additionally, they can be used to prevent a known problem referred to as signal swapping as will be described later.
To summarise, Equation (16) represents initialisation by Ũk of orthogonalisation of information (vectors or signals) in the data matrix Xk+1. This corresponds to step 22 of
Attention will now turn to the problem of signal swapping, which is illustrated in
The phenomenon of signal swapping can potentially occur in the second order stage 22/23 (Equations (17) to (22)) and/or the higher order stage 25/27. It is counteracted by applying only small rotation updates in the initialisation process previously described. 10 Here (as mentioned earlier) “small” means there is more than one possible rotation and the smaller/smallest is selected. For the purposes of the discussion to follow, x and y are vectors and each is a respective row of {tilde over (V)}kT representing all samples of a signal from a respective sensor extending the full length of a data window. After x and y are initialised in the second order stage 22/23, they are close to being orthogonal to one another. This is shown in
To make the vectors x and y more accurately orthogonal, a relative rotation angle θ is introduced between them. This rotation is implemented by a decorrelation method such as Jacobi diagonalisation of a symmetric matrix as described earlier. For signals which are already close to orthogonal, two solutions are obtained by this method, i.e. θ≈0 and θ≈90 degrees. To avoid signal swapping the smaller value of θ≈0 should be chosen. If the signals are indeed rotated relative to one another by the smaller value of θ, this corresponds to a Givens rotation matrix given by:
When θ≈0 then a pair of initialised (nearly orthogonal) vectors x and y can be updated using:
In this case, the updated signals are in the same order and signal swapping has been avoided. This is in direct contrast to the case of θ≈90, for which a Givens rotation matrix for rotating one vector relative to another is:
Updating a pair of initialised vectors x and y using Equation (33) gives:
In this case, the updated signals have swapped as indicated by the exchange of positions of xT and −yT, and an inversion also occurs indicated by the negative sign of −yT.
In the Jacobi decorrelation method, the smaller of the two angles is therefore chosen to achieve convergence to a solution. Thus, signal swapping is prevented in accordance with the invention when using an initialisation process for the second order stage and Jacobi with smaller angle updates. For validity, this approach relies on the system under consideration being slowly changing and resulting in the value of θ being small.
Signal swapping can potentially occur in the higher order statistics stage also, as again there will be more than one solution or rotation angle to choose from. After Rk is applied to data using Equation (27), estimated signal vectors will be related to corresponding source signals by a small rotation. In the BLISS or higher order statistics stage, since fourth order statistics are used, there are now four angles that can be chosen to achieve signal separation. This is illustrated in
As has been said, BLISS is based on fourth order statistics. It is also possible to use a technique based on third order statistics instead of BLISS, as described in “An Improved Cumulant Based Method for Independent Component Analysis”, T. Blaschke and L. Wiskott, Proc. Int. Conf. on Artificial Neural Networks, ICANN'02, 2002. It is also possible to use a combination of third and fourth order statistics.
The invention as described above was used for fetal ECG analysis using the data shown analysed in
The initialisation and update processes also allow desired signals to be targeted alone. This technique avoids redundant computation incurred by tracking unwanted signals. It exploits the fact that signals of interest (e.g. maternal ECG and fetal ECG) can be classified in an acquisition stage. For fetal ECG analysis, this classification process could either depend on the choice of the user (clinician or doctor) or on a pattern classification technique. Pattern recognition methods for ECG signals are described by J. Pan and W. J. Tompkins in “A Real-Time QRS Detection Algorithm”, IEEE Transactions on Biomedical Engineering 32(3): 230-236, 1985.
Once the desired signals are identified, this information can be used to focus or target processing on them. After classifying signals of interest in the acquisition stage 16, as compared to the invention without such targeting, the targeted technique only differs in restricting updating to desired signals in the signal separation stage 17.
Targeting of desired signals at 73 and 77 in
Before describing the targeting method in more detail, it is re-emphasised that a higher order stage 76 to 78 can be implemented using an iterative sequence of pairwise rotations. Each of these rotations is designed to maximise the statistical independence of two signals forming pair of signals, statistical independence being measured by an associated pairwise cost function. With m signals v1 to vm, there are mC2 or m(m−1)/2 signal pairs. After initialisation of these signals to form orthonormal vectors v′1 to v′m with zero mean at 76, the total number of vector pairs are expressed as:
In (35), it is assumed that the orthonormal vectors v′1 to v′m are initialised using higher order information obtained using an immediately preceding data window as described previously. All of these pairs would have to be updated in an iterative sequence in the absence of targeting, which means that m(m−1)/2 pairs would need updating. In contrast to this, the targeting approach requires only to update pairs that include the desired vectors. This is possible as:
If pairs involving v′1 and v′2 are targeted and no others, only the first two rows of (35) need updating. This involves only (m−1)+(m−2) signal pairs, as opposed to m(m−1)/2 pairs for a full update. In this ‘targeted’ process, only the maternal and fetal signals will be separated. Separation of other signals is not accomplished, as in this example they are not required.
As described above, the targeting principle can also be applied to the second order stage 72-74. After an acquisition stage and by initialising second order information as in Equation (11) to produce an initialised covariance matrix, rows of this matrix that correspond to desired signals will be updated, but not other rows.
Referring now to
Application of the invention to produce the results shown in
After the acquisition stage, only desired signals are processed to full separation. Once the order of the separated signals has been determined in the acquisition stage, as shown in the triangle of signals in (35), the signals can be re-ordered such that the desired signals occupy the top rows of the triangle: i.e. signals 115 and 121 in
In a second example, the targeted method of the invention was applied to a singleton fetal ECG recording, i.e. an ECG of a single fetus. The recording employed twelve signal electrodes and associated ground and reference electrodes and lasted for three minutes. In
Another method for dealing with signal swapping is to use overlapping windows and a cross-correlation of overlapping estimated signals which result. For example, let overlapping sections of estimated signals in windows k and k+1 be f samples long and be represented by (m×f) matrices Ŝok and Ŝok respectively. The rows of these matrices contain the estimated signals and the overlapping sections will contain similar signals. The matrix Ŝok is related to Ŝok+1 by:
Ŝok+1≈DPŜok (36)
where D denotes an (m×m) matrix, whose diagonal elements are 1 or −1 (providing for signal inversions) and P denotes an (m×m) permutation matrix (providing for signal swapping). These conditions arise in adjacent windows as the estimated signals can be ordered differently and inversions are also possible. The overlapping sections of the estimated signals are ordered until this correlation matrix resembles a diagonal matrix. In this case, the signals will be in the same order. This technique is described in by G. Spence in “A Joint Estimation Approach for Periodic Signal Analysis”, PhD thesis, Queen's University Belfast, June 2003.
Referring now to
The lead box 173 contains processing circuitry for signals from the signal electrodes 171a etc. are shown inset at 173a: circuitry is shown for two signal electrodes 171a and 171b, dots D indicate like circuitry for other signal electrodes. The circuitry comprises for each signal electrode 171a etc. a respective low-noise differential amplifier 175 with an analogue high-pass filter and an analogue low-pass anti-aliasing filter both incorporated in a filter unit 176. The filters 176 are connected to a common multi-channel analogue to digital converter (ADC) 177. The amplifiers 175 subtract a signal from the reference electrode R from signals from respective signal electrodes 171a etc., and amplify the resulting difference signals. The difference signals are converted to digital signals by the ADC 177, and are then recorded and processed using the computer 174 to separate signals as described earlier.
Referring now to
Twelve signals 201 to 212 (marked on left) in
In
It is possible in principle to envisage implementing signal separation without initialisation and update processes, but with the above ordering technique to correct the problem of signal swapping. However, this has disadvantages. It would sacrifice the ability of the invention to exploit the computationally efficient targeted approach. Additionally, to provide reliable estimates in the cross-correlation technique, a large window overlap would have to be used. This means, that in order to process the data, a larger number of windows is required: e.g. with 50% Window overlap the data is processed twice (and takes twice as long) compared to processing once (i.e. with no overlap). This is disadvantageous because it lengthens the time between detection of mixed signals by sensors 14 and generation of unmixed signals. In the case of fetal ECG in particular, this is important because it may be necessary for a clinician to intervene rapidly if a heartbeat abnormality is detected, and it may be too late to save the fetus if processing time is too long.
The foregoing example of the invention may also be applied to complex-valued data. The principles of initialising the second and higher order stages and using small updates remain unchanged; however complex-valued versions of the second order stage (Jacobi) and higher order (ICA) stage are required. The use of ICA for complex-valued data is well researched and many decorrelation and ICA techniques have been developed for this problem. For more information on a complex-valued Jacobi refer to: S. Haykin, “Adaptive filter theory, second edition”, Prentice Hall, ISBN 0-13-013236-5, 1991.
In addition, more information about complex-valued ICA methods can be found in:
P. Comon, “Independent component analysis, A new concept”, Signal Processing 36, pp. 287-314, 1994.
I. J. Clarke, “Exploiting non-Gaussianity for signal separation”, MILCOM 2000. 21st Century Military Communications Conference Proceedings, Volume: 2, 22-25 Oct. 2000.
In the last of the above three references, Clarke showed that an extended version of BLISS can deal with complex signals by treating their in-phase and quadrature components as individual real signals. By using Jacobi and this extended version of BLISS small updates of the initialised signals are again possible. This is because the smaller angle is chosen in the Jacobi technique and the rotation of the initialised signals in the higher order stage will again depend on fourth order statistics.
The equations for a complex-valued version of the invention, in the (k+1)th window, can be summarised in four stages as follows. The superscript index H denotes the hermitian of a matrix or vector. These equations can also be used to represent the real-valued case.
1. Initialisation in second order stage:
ŨkH(Xk+1Xk+1H)Ũk=ŨkH(Ũk+1{tilde over (α)}k+1Ũk+1H)Ũk (37)
2. Update using complex-valued Jacobi:
ŨzH(ŨkHXk+1Xk+1HŨk)Ũz={tilde over (α)}k+1 (38)
Ũk+1=ŨkŨz (39)
{tilde over (V)}k+1H={tilde over (α)}k+1−1/2Ũk+1HXk+1 (40)
3. Initialisation in higher order stage:
Rk{tilde over (V)}k+1H=RkRk+1HSk+1 (41)
4. Update using complex-valued ICA
Ŝk+1=Rz(Rk{tilde over (V)}k+1H) (42)
Rk+1=RzRk (43)
Equations (39) and (40) represent updated eigenvectors and orthonormal signals respectively. Equations (42) and (43) represent updated independent signals and the higher rotation matrix respectively. The four stage process is summarised as described earlier with reference to
Equations given in the foregoing description for calculating intermediate quantities and results can clearly be evaluated by an appropriate computer program embodied in a carrier medium and running on a conventional computer system. Such a program and system are straightforward for a skilled programmer to implement without requiring invention, and will therefore not be described further.
A further embodiment of the invention will now be described. When processing data from sensors in real-time, only one or a few data snapshots may be available to update a data matrix at a time: here a “snapshot” is defined as a set of simultaneous data samples, a respective single sample of a signal from each of the sensors, i.e. a column of the data matrix. In such cases, any moving data window approach as described earlier will require a large amount of window overlap, i.e. the position of the window can only be moved a limited amount between successive processing stages. Thus a large number of windows may be required to cover the data: for each window, this corresponds to a significant degree of and perhaps a large amount of redundancy as the window's statistics (although initialised to the desired solution) will be re-computed. One approach for dealing with this for a sliding window version of the invention (the example described earlier) is to subtract the statistics of the oldest snapshot(s) from the overall statistics. The statistics of the latest added snapshot(s) are then computed and added to the difference between the overall statistics and those of the oldest snapshot(s). There are also two ways of implementing the BLISS algorithm: one version operates directly on sensor signals or digitised equivalents, and the other operates directly on a fourth order tensor derived from those signals. Sets of results from these versions will be equivalent and may be used for the moving window approach of the preceding embodiment.
The redundancy problem may also be avoided by using a recursive update approach: in this approach signal separation is implemented using weighted measures of statistics for old data computed earlier (and not recomputed) and newly computed statistics obtained from most recently added data. This approach is not conceptually dissimilar to the moving window approach previously described. It once more uses initialisation and small update processes, and again these are applied in the second and higher order statistics stages. However, instead of using a moving window approach with a rectangular window, an exponential fading window is used. This allows the invention to process the latest available data snapshot(s) recursively while avoiding the problem of recomputing overlapping regions of windows. It then follows that redundancy is minimal, because the statistics of overlapping window sections are not recomputed.
In this embodiment, the recursive update approach is divided into two parts. The first part uses a rectangular window of data as in the previous embodiment, and this is referred to as the acquisition stage. In the acquisition stage signals are separated as previously described. The first, second and higher order statistics of data are calculated in the acquisition stage for use in a second stage, in which these statistics are recursively updated (see Tables 2 to 4). In the acquisition stage which spans n samples, the data is defined as an (m×n) data matrix X0=[x1 x2 . . . xm]T with rows x1 to xm each representing a respective sensor output signal. X0 is therefore represented by m vectors x1 to xm each with n elements. The subscript index 0 of X0 and other parameters in Table 1 indicates t=0 for the acquisition stage. Fourth order moments which form a tensor are set out in Table 1 below: these provide a measure of the relationship between four vectors of arbitrary kind e.g. yi, yj, yk and yl such that tensor elements Mijkl are given by:
where yi(τ) denotes the τth sample of vector yi: here τ is an integer time index, and expresses a sample time when multiplied by a time interval Ts between successive samples.
Each element Mijkl of M0 is defined as:
where vi(τ) is a sample of the ith signal vector v1 at time τ, τis in units of a sampling time interval, and n is the number of samples per vector as previously indicated.
Using the definitions in Table 1, the second or recursive update stage of this example may be defined mathematically as follows. For convenience of explanation this is described in terms of successive updates each made using data comprising the preceding data window with the contributions of previous snapshots diminished by weighting and addition of a single snapshot of new data per update; this means that progressively older data is progressively more diminished. Each snapshot has an index number t, where t=1 to T in units of a sampling time interval, and the first or acquisition stage is denoted by t=0. This embodiment of the invention is not however limited to a single snapshot per update.
For convenience, the data snapshot (hitherto xt) will now be defined as dt, and the decorrelated and normalised equivalent of the snapshot obtained from the second order statistics stage (hitherto vt) will now be defined as zt.
After acquisition and for each snapshot dt, a respective mean vector vt is recursively updated using Equation (50) below and subtracted from the data snapshot using Equation (51) below. In the following Tables 2-4, processes in the recursive update stage are defined. For cases in which each update involves more than one snapshot, the mean is calculated over all snapshots in the update.
Equations (50) and (51) provide for first order statistics to be generated from an exponentially fading window of data by the use of p and (1−p), which are referred to in signal processing as “forget factors”. The value of p lies between 0 and 1, and so p and (1−p) progressively reduce the effect of older data. In this example p is likely to lie between 0 and 0.1. Exponential fading occurs because (1−p) and p are applied on each update, so the more updates that follow a snapshot the more its effect diminished: e.g. after u updates the wth updating snapshot is diminished by p(1−p)u−w+1, snapshots in the acquisition stage window being ignored as regards the value of w. The value of the forget factor p may be set in relation to the acquisition data window length. It also provides the advantage that p may be set adaptively to reflect confidence in data: e.g. the value of p may be adjusted upwards or downwards during processing depending on whether the data's noise component is reducing or increasing.
The next part of the processing procedure in the recursive update stage employs Equations (52) to (56) in Table 3 below. In this part, second order statistics expressed by the data matrix's covariance matrix χt are recursively updated using Equation (52): here the covariance matrix χt for time t is produced by adding (1−p) times the covariance matrix χt−1 for time (t−1) to p times the covariance matrix {tilde over (d)}t{tilde over (d)}tT of the new data snapshot {tilde over (d)}t to provide an exponentially fading window once more. For the first update, t−1=0 indicates that the relevant preceding results are from the acquisition stage described with reference to Table 1. The covariance matrix χt has orthogonal components which are initialised in Equation (53) by premultiplication by Ut−1T and postmultiplication by Ut−1, as described in the preceding embodiment (see Equation (16)). Here and subsequently xr. (r=1, 2 . . . ) will be used to denote the r mode product, such that for the second order case XXT=UλUT=λ×1U×2U, where U, λ and T are as defined earlier. The expression λ×1U×2U therefore means the matrix U is applied to the first dimension of λ and then to the second dimension of λ. As in the preceding embodiment, the objective in the present embodiment is to diagonalise the second order data covariance matrix, albeit this matrix is generated from a recursive update in Equation (52) instead of by processing a whole data window. In Equation (54) the matrix Ut−1TχtUt−1 is diagonalised using small updates and the results of updating are defined in Equations (55) and (56). Equation (56) produces a second order estimate zt of a single snapshot, because in this embodiment each update involves a single snapshot. If each update is implemented with multiple snapshots, the estimate will be for multiple snapshots also.
The next (higher order statistics) part of this embodiment employs Equations (57) to (62) in Table 4 below. The current estimate of the snapshot from Equation (56) in the second order stage is passed to the higher order stage. The fourth order moment tensor of this snapshot is formed in Equation (57), and it is recursively updated. This tensor is approximately diagonalised by application of Rt−1 in Equation (59), and small updates Rz are then used to diagonalise it more fully in Equation (60). The diagonalisation of this fourth order tensor is a requirement for signal separation in this example. Again, these can be calculated using the BLISS algorithm. The rotation matrix is determined using Equation (61), and an estimate of the signals (for time t) is given by Equation (62). These are the same types of processes as in the earlier embodiment, except that now fourth order statistics are generated from a recursive update.
m×m×m×m
m×m×m×m
In practice, recursive updating may be implemented using more than one snapshot per update, in which case this embodiment yields as many estimated signal snapshots as there are snapshots per update.
In this embodiment, Equations (50), (52) and (58) represent statistics that are generated from signals that span an exponentially fading window, i.e. using a forget factor to bias processing. For the moving window approach in the first example described with reference to
Analysis of electroencephalogram (EEG) measurements by a clinician is a well established technique for diagnosing brain dysfunction such as epilepsy, brain tumour and brain injury. However, in practice, this analysis can be very difficult as EEG measurements produce weak electrical signals which are susceptible to severe corruption by unwanted electrical interference (artefact): such interference arises from sources such as eye movement and/or blink (ocular artefact) and muscle activity of the patient and from mains power supply interference.
A BSS approach may be used to suppress unwanted interference in order to isolate signals of interest in EEG data. When desired EEG signals are repetitive (‘pseudo-continuous’) a tracking BSS method can be used to track them, and thus allow long-term monitoring. Clinical EEG data were obtained from a patient suffering from repetitive epileptic seizures. These data yielded repetitive EEG signals, i.e. EEG voltage against time: they were obtained conventionally by applying eight scalp electrodes to a patient and monitoring electrode voltages with EEG equipment. The EEG signals consisted of repeated epileptic seizure signals adulterated by interference. Epileptic seizures in the signals could be inferred but not reliably discerned from interference; for example, counting signal maxima for any of the eight signals did not yield the correct number of seizures. The recursive tracking embodiment of the invention was applied to the EEG signals to provide an estimated signal. An acquisition stage was used that spanned five seconds and each update was produced using a respective set of 128 data snapshots with a sample rate of 256 Hz, i.e. 0.5 second to produce 128 data samples. The estimated signal did not suffer from signal swapping and it showed repetitive epileptic seizures indicated by easily discernable spikes of greater amplitude than interference. The signal to interference and noise ratio was in the range 3 to 6 as the spikes varied in height.
| Number | Date | Country | Kind |
|---|---|---|---|
| 0326539.4 | Nov 2003 | GB | national |
| Filing Document | Filing Date | Country | Kind | 371c Date |
|---|---|---|---|---|
| PCT/GB2004/004714 | 11/9/2004 | WO | 00 | 5/11/2006 |
| Publishing Document | Publishing Date | Country | Kind |
|---|---|---|---|
| WO2005/052848 | 6/9/2005 | WO | A |
| Number | Name | Date | Kind |
|---|---|---|---|
| 5776073 | Garfield et al. | Jul 1998 | A |
| 5959966 | Torkkola | Sep 1999 | A |
| 6321200 | Casey | Nov 2001 | B1 |
| 6751498 | Greenberg et al. | Jun 2004 | B1 |
| 20010044719 | Casey | Nov 2001 | A1 |
| Number | Date | Country |
|---|---|---|
| WO 0117109 | Mar 2001 | WO |
| WO 0127874 | Apr 2001 | WO |
| WO 02088907 | Nov 2002 | WO |
| WO 03028550 | Apr 2003 | WO |
| WO 03090127 | Oct 2003 | WO |
| Number | Date | Country | |
|---|---|---|---|
| 20070088226 A1 | Apr 2007 | US |
