Abstract:
We describe a method of blind source separation for use, for example, in a listening or hearing aid. The method processes input data from multiple microphones each receiving a mixed signal from multiple audio sources, performing independent component analysis (ICA) on the data in the time-frequency domain based on an estimation of a spectrogram of each acoustic source. The spectrograms of the sources are determined from non-negative matrix factorization (NMF) models of each source, the NMF model representing time-frequency variations in the output of an acoustic source in the time-frequency domain. The NMF and ICA models are jointly optimized, thus automatically resolving an inter-frequency permutation ambiguity.

Description:
RELATED APPLICATIONS 
     This application is a continuation of U.S. patent application Ser. No. 14/678,419, filed 3 Apr. 2015, which is incorporated herein in its entirety. 
    
    
     FIELD OF THE INVENTION 
     This invention relates to methods, apparatus and computer program code for blind source separation, for example to assist listeners with hearing loss in distinguishing between multiple different simultaneous speakers. 
     BACKGROUND TO THE INVENTION 
     Many people&#39;s ability to understand speech is dramatically reduced in noisy environments such as restaurants and meeting rooms. This is especially true of the over 50s, and it is one of the first signs of age-related hearing loss, which severely curtails people&#39;s ability to interact in normal social situations. This can lead to a sense of isolation that has a profound influence on general lifestyle, and many studies suggest that it may contribute to dementia. 
     With this type of hearing loss, listeners do not necessarily suffer any degradation to their threshold of hearing; they could understand the speech perfectly well in the absence of the interfering noise. Consequently, many sufferers may be unaware that they have a hearing problem and conventional hearing aids are not very effective. Also, sufferers become more dependent on lip reading, so any technical solution should preferably have a low latency to keep lip sync. 
     Techniques are known for blind source separation, using independent component analysis (ICA) in combination with a microphone array. In broad terms such techniques effectively act as a beamformer, steering nulls towards unwanted sources. In more detail there in an assumption that the signal sources are statistically independent, and signal outputs are generated which are as independent as possible. The input signals are split into frequency bands and each frequency band is treated independently, and then the results at different frequencies are aligned. This may be done by assuming that when a source is producing power at one frequency it is probably also active at other frequencies. However this approach can suffer from problems if power at one frequency may be correlated with the absence of power at another frequency, for example the voiced and unvoiced parts of speech. This can lead to frequency bands from different sources being swapped in the output channels. 
     It is also known to employ non-negative matrix factorisation (NMF) to distinguish between speakers. In broad terms this works by learning a dictionary of spectra for each speaker. However whilst this can distinguish between characteristically different voices, such as male and female speakers, it has difficulty with finer distinctions. In addition this approach can introduce substantial latency into the processed signal, making it unsuitable for real time use and thus unsuitable, for example, to assist in listening to a conversation. 
     Accordingly there is a need for improved techniques for blind source separation. There is a further need for better techniques for assisting listeners with hearing loss. 
     SUMMARY OF THE INVENTION 
     According to the present invention there is therefore provided a method of blind source separation, the method comprising: inputting acoustic data from a plurality of acoustic sensors, said acoustic data comprising acoustic signals combined from a plurality of acoustic sources; converting said acoustic data to time-frequency domain data, wherein said time-frequency domain data is represented by an observation matrix X f  for each of a plurality of frequencies f; performing an independent component analysis (ICA) on said observation matrix X f  to determine a demixing matrix W f  for each said frequency such that an estimate Y f  of the acoustic signals from said source at said frequencies f is determined by X f  W f ; wherein said ICA is performed based on an estimation of a spectrogram of each said acoustic source; wherein said spectrogram of each said acoustic source is determined from a model of the corresponding acoustic source, the model representing time-frequency variations in a signal output of the acoustic source in the time-frequency domain. 
     In broad terms embodiments of the method employ a model that captures variations in both time and frequency (across multiple frequency bands), and then the independent component analysis effectively learns a decomposition which fits this model, for each source aiming to fit a spectrogram of the source. In this way the inter-frequency permutations of the demixing matrix are automatically resolved. Preferred embodiments of the model represent the behaviour of an acoustic source in statistical terms. 
     In principle the acoustic source model may be any representation which spans multiple frequencies (noting that the model may, for example, define that some frequencies have zero power). For example PCA (principle component analysis) or SVD (singular value decomposition) may be employed. However it is particularly advantageous to use an NMF model as this better corresponds to the physical process of adding together sounds. In principle a single component NMF model could be used, but preferably the NMF model has two or more components. 
     In preferred embodiments of the method the (NMF) model and independent component analysis (ICA) are jointly and iteratively improved. That is ICA is used to estimate the signals from the acoustic sources and then the (NMF) model is updated using these estimated signals to provide updated spectrograms, which are in turn used to once again update the ICA. In this way the ICA and NMF are co-optimised. 
     In some approaches the ICA and NMF models are updated alternately, but this is not essential—for example several updates may be performed to, say, the NMF model and then the ICA model may be updated, for example to more accurately align the permutations amongst frequency bands to sources. In practice, however, it has been found that interleaving updates of different types tends to approach the target joint optimisation faster, in part because the initial data tends to be noisy and thus does not benefit from an attempt to impose too much structure initially. 
     In some preferred implementations of the method the updating of the independent component analysis includes determining a permutation of elements of the demixing matrix over the acoustic sources prior to determining updated spectrograms (σ k ) for the acoustic sources. It is not essential to perform such a permutation but it can be helpful to avoid the method becoming trapped in local maxima. In broad terms the additional step of performing the permutation alignment based on the spectrogram helps to achieve a good fit of the NMF model more quickly (where frequency bands of different sources are cross-mixed the NMF model does not fit so well). 
     Preferably the updating of the ICA includes adjusting each of the demixing matrices W f  according to a gradient ascent (or descent) method, where the gradient is dependent upon both the estimate of the acoustic signals from the sources and the estimate of the spectrograms of the sources (from the NMF model). In broad terms this gradient search procedure aims to identify demixing matrices which make the output data channels (Y f ) look independent given (i.e. in) the NMF model representation. 
     In embodiments the NMF model is defined by latent variables U (a frequency-dependent spectral dictionary for each source) and V (time-dependent dictionary activations) for each acoustic source (noting that U and V are tensors). The NMF model is updated by updating these latent variables. In embodiments this may be performed in two stages—identify the best dictionary given a set of activations; and identify the best set of activations given a dictionary. Preferably the dictionaries and activations are jointly optimised with the demixing matrix for each frequency although, as previously noted, the update steps need not be performed alternately. 
     In embodiments the NMF model factorises time-frequency dependent variances to a power λ(λ≠2), σ k   λ  rather than σ k   2  because the ICA effectively performs a spatial rotation to decouple the signal components and with Gaussian data the squared power results in a circular cost contour which does not distinguish rotations. In some preferred embodiments λ=1 as this provides a good balance between some cost for rotation and avoiding being trapped in local maxima. 
     In broad terms embodiments of the method allocate audio sources to channels. However if too many channels are available the method may split a source over multiple channels, and fragmenting real sources in this way is undesirable. In embodiments, therefore, the method may preprocess the acoustic data to reduce an effective number of acoustic sensors to a target number of virtual sensors, in effect reducing the dimensionality of the data to match the actual number of sources. This may either be done based upon knowledge of the target number of sources or based on some heuristic or assumption about the number of likely sources. The reduction in the dimensionality of the data may be performed by discarding data from some of the acoustic sensors or, more preferably, may employ principal component analysis with the aim of retaining as much of the energy and shape of the original data as possible. 
     Preferably the method also compensates for a scaling ambiguity in the demixing matrices. In broad terms, whilst the independent component analysis may make the outputs of the procedure (source estimates) substantially independent, the individual frequency bands may be subject to an arbitrary scaling. This can be compensated for by establishing what a particular source would have sounded like at one or more of the acoustic sensors (or even at a virtual, reference acoustic sensor constructed from the original microphones). This may be performed by, for example, using one of the acoustic sensors (microphones) as a reference or, for a stereo output signal, using two reference microphones. The scaling ambiguity may be corrected by selecting time-frequency components for one of the output signals and by calculating the inverse of the estimated demixing matrix, since the combination of the demixing matrix and its inverse should reproduce what the source would have sounded like on (all the acoustic sensors). By employing this approach (eq(25) below) a user may determine what the selected source would have sounded like at each microphone. The user, or the procedure, may select a microphone to “listen to” (for example by selecting a row of the output data Y f  (k) corresponding to source estimate k), or for stereo may select two of the microphones to listen to. 
     Embodiments of the method perform the blind source separation blockwise on successive blocks of time series acoustic data. However the labelling of sources k may change from one block to the next (for example if W, U, V are initialised randomly rather than based on a previous frame). It is therefore helpful to be able to identify which real source corresponds to which source label k in each block, to thereby partially or wholly remove a source permutation ambiguity. The skilled person will appreciate that a number of different techniques may be employed to achieve this. For example in some applications a desired target source may have a substantially defined or fixed spatial relationship to the acoustic sensors (microphone array), for example when it is desired to target speech output from a driver or passenger of a car. In another approach a loudest source (that is the direction from which there is most audio power) may be assumed to be the target source—for example where it is desired to distinguish between a speaker and background from an air conditioning unit in, say, a video conference setting. Alternatively characteristics of a source (for example from the NMF model) may be used to identify a target source. 
     In one preferred approach the system or a user may select a target direction for a target source to be selected. Then the procedure may identify the source which best matches this selected direction. This may be performed by, for example, selecting the source with the highest phase correlation between the microphone array phase response from the selected direction and the corresponding portion (row) of the set of demixing matrices. 
     The skilled person will appreciate that embodiments of the procedure directly produce source estimates (Y) or a selected source estimate (Y(k)), albeit in the time-frequency domain. Such a source estimate may be used in the time-frequency domain or may be converted back to the time domain. Alternatively, however, the demixing matrices W may be converted from the time-frequency domain to the time domain, for example using an inverse Fourier transform, to determine a time domain demixing filter. 
     The calculations to implement the above described blind source separation technique can be relatively time consuming (for example taking around 1 second on a current laptop), but it is desirable for embodiments of the method, in particular when used as a hearing aid, to be able to operate in substantially real time. To achieve this the procedure may be operated at intervals on sections of captured acoustic data to establish coefficients for a demixing filter, that is to determine the demixing matrices W f . These coefficients may then be downloaded at intervals to a configurable filter operating in real time, in the time domain, on the acoustic data from the acoustic sensors. 
     In a related aspect, therefore, the invention provides apparatus to improve audibility of an audio signal by blind source separation, the apparatus comprising: a set of microphones to receive signals from a plurality of audio sources disposed around the microphones; and an audio signal processor coupled to said microphones, and configured to providing a demixed audio signal output; the audio signal processor comprising: at least one analog-to-digital converter to digitise signals from said microphone to provide digital time-domain signals; a time-to-frequency domain converter to convert said digital time domain signals to the time-frequency domain; a blind source separation module, to perform audio signal demixing in said time-frequency domain to determine a demixing matrix for at least one of said audio sources; and a digital filter to filter said digital time-domain signals in the time domain in accordance with filter coefficients determined by said demixing matrix, wherein said filter coefficients are determined asynchronously in said time-frequency domain; and wherein said audio signal processor is further configured to process said demixing matrix to select one or more said audio sources responsive to a phase correlation determined from said demixing matrix. 
     In embodiments the apparatus may be configured to resolve a scaling ambiguity in the demixing matrix as previously described and/or to reduce dimensionality of the input audio signal prior to demixing. Preferably the blind source separation module is configured to perform joint ICA and NMF processing to implement the audio signal demixing. 
     A demixed audio signal output, typically from an automatically or manually selected source, may be output and/or used in many ways according to the application. For example where the system is used as a listening aid the audio output may be provided to headphones, earbuds or the like, or to a conventional hearing aid. Alternatively the audio output may be provided to other electronic apparatus such as a video conferencing system or fixed line or mobile phone (for example, with an in-vehicle communications system). 
     The skilled person will appreciate that in the above described apparatus the audio signal processor may be implemented in hardware, firmware, software or a combination of these; on a dedicated digital signal processor, or on a general purpose computing system such as a laptop, tablet, smartphone or the like. Similarly the blind source separation module may comprise hardware (dedicated electronic circuitry), firmware/software, or a combination of the two. 
     In a further aspect the invention provides a method of blind source separation, the method comprising: processing an observation matrix X f  representing observations of signals at a plurality of frequencies f from a plurality of acoustic sources using a demixing matrix W f  for each of said frequencies to determine an estimate of demixed signals from said acoustic sources Y f  for each of said frequencies, the processing comprising iteratively updating Y f  from X f  W f ; wherein said processing is performed based on a probability distribution p(Y tkf ; σ tkf ) for Y dependent upon 
               1     σ   tkf   2       ⁢     ⅇ     -              Y   tkf          λ       σ   tkf   λ                 
where t indexes time intervals and k indexes said acoustic sources or acoustic sensors sensing said acoustic sources; and wherein σ tkf  are variances inferred from a non-negative matrix factorisation (NMF) model where
 
               σ   tkf   λ     =       ∑   l     ⁢       V   ltk     ⁢       U   lfk     .               
where l indexes non-negative components of said NMF model, U and V are latent variables of said NMF model, and λ is a parameter greater than zero.
 
     As previously described, in broad terms the NMF model imposes structure on the variances, noting that σ k  is an approximation of the spectrogram of source k across frequencies, in particular because the dictionaries and activations defined by the latent variables U, V form a sparse representation of the data. In broad terms the NMF model (or potentially some other model) represents the spectrogram of source k as the time varying sum of a set of dictionary components. The ICA model expresses the probability of the source estimates given the variances imposed by the NMF model. 
     The signal processing determines a set of demixing matrices Wand values for the latent variables U, V which (preferably) optimally fit the above equations. Thus, as previously described, embodiments of the procedure iteratively update U, V and W, improving each tensor given the other two. Preferred implementations of the procedure also apply an optimal permutation to W, as this can in effect provide a relatively large step in improving Was compared with gradient ascent. 
     Thus in embodiments of the processing the following steps are applied, potentially more than once each, and potentially in any order: update W using permutation and/or scaling; update W using a gradient-based update; update U; update V. The updates of W may be used to determine updated source estimates (for updating U, V). The updates of U and V may be used to determine updated source spectrogram estimates (for updating W). Optionally U and V may be initialised based on prior information, for example a previously learnt, stored or downloaded dictionary. 
     The skilled person will appreciate that embodiments of the above described methods may be implemented locally, for example on a general purpose or dedicated computer or signal processor, phone, or other consumer computing device; or may be partly or wholly implemented remotely, in the cloud, for example using the communication facilities of a laptop, phone or the like. 
     The invention further provides processor control code to implement the above-described apparatus and methods, for example on a general purpose computer system or on a digital signal processor (DSP). The code is provided on a non-transitory physical data carrier such as a disk, CD- or DVD-ROM, programmed memory such as non-volatile memory (eg Flash) or read-only memory (Firmware). Code (and/or data) to implement embodiments of the invention may comprise source, object or executable code in a conventional programming language (interpreted or compiled) such as C, or assembly code, or code for a hardware description language. As the skilled person will appreciate such code and/or data may be distributed between a plurality of coupled components in communication with one another. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       These and other aspects of the invention will now be further described, by way of example only, with reference to the accompanying figures in which: 
         FIG. 1  shows an example acoustic environment to illustrate the operation of a system according to an embodiment of the invention; 
         FIG. 2  shows the architecture of apparatus to improve audibility of an audio signal by blind source separation; 
         FIGS. 3 a  and 3 b    show, respectively, an example spatial filter for the apparatus of  FIG. 2 , and an STFT (short time fourier transform) implementation of time-frequency/frequency-time domain conversions for the system of  FIG. 2 ; 
         FIG. 4  shows modules of a frequency domain, filter-determining system for the apparatus of  FIG. 2 ; 
         FIG. 5  shows a flow diagram of a procedure for blind source separation according to an embodiment of the invention; and 
         FIG. 6  shows a general purpose computing system programmed to implement the procedure of  FIG. 5 . 
     
    
    
     DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 
     Broadly speaking we will describe techniques for blind source separation on the audio outputs of a small microphone array to separate a desired source from one or more interfering sources. In one application a user can listen to the desired source in real time over headphones or via a hearing aid. However the technology is not just applicable to listening aids and can be useful in any application where a sensor array is measuring a linear convolutive mixture of sources. In audio this includes applications such as teleconferencing and machine hearing. 
     By way of example, consider the acoustic scene of  FIG. 1 . This comprises four sources s 1 -s 4  with respective audio channels h 1 -h 4  to a microphone array  10  comprising (in this example) 8 microphones. The aim is to demix the microphone signals to make estimates of the original sources—that is to perform Blind Source Separation or Blind Signal Separation (BSS). We assume minimal information about the sources and the microphone locations. In some applications the microphone array may be placed on a table or chair in a social setting or meeting and embodiments of the systems we describe are used to separate a desired source, such as a person speaking, from undesired sounds such as other speakers and/or extraneous noise sources. 
     Using the multi-channel observations x, the task is to design a multi-channel linear filter w to create source estimates y. 
     
       
         
           
             
               
                 
                   
                     
                       y 
                       _ 
                     
                     
                       t 
                       ′ 
                     
                   
                   = 
                   
                     
                       ∑ 
                       τ 
                     
                     ⁢ 
                     
                       
                         w 
                         τ 
                       
                       ⁢ 
                       
                         
                           x 
                           _ 
                         
                         
                           
                             t 
                             ′ 
                           
                           - 
                           τ 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   1 
                   ) 
                 
               
             
           
         
       
     
     Given the lack of location information, rather than recover the original sources the objective is to recover the sources up to a permutation ambiguity P and an arbitrary linear transformation b k,τ ,
 
 y   P(k),t′ ≈Σ τ   b   k,τ   s   k,t′ −τ  (2)
 
where s k  labels source k.
 
STFT Framework
 
     Overlapped STFTs provide a mechanism for processing audio in the time-frequency domain. There are many ways of transforming time domain audio samples to and from the time-frequency domain. The NMF-ICA algorithm we describe can be applied inside any such framework; in embodiments we employ Short Time Fourier Transforms (STFT). Note that in multi-channel audio, the STFTs are applied to each channel separately. 
     Within this framework we define:
         K is the number of channels.   F is the number of STFT frequencies.   T is the number of STFT frames.       

     In the STFT domain, the source estimate convolution eq(1) becomes matrix multiplication. At each frequency we have the T×K observation matrix X ƒ , and an unknown demixing matrix W ƒ  such that the demixed output Y ƒ  is given by
 
 Y   ƒ   =X   ƒ   W   ƒ   (3)
 
     Here the demixed output Y ƒ  is also a T×K matrix where k labels sources, and the demixing matrix W ƒ  is a K×K matrix (Xε             T×K×F , Yε           T×K×F , Wε           K×K ).
     The equivalent objective to eq(2) is
 
 Y   ƒ   ≈S   ƒ   B   ƒ   P   (4)
 
where B ƒ  is an arbitrary diagonal scaling matrix, and P is a global permutation matrix. The task of Blind Source Separation is to use knowledge of the statistics of audio to estimate W ƒ  for each frequency.
 
Notation
                    means equal up to a constant offset (which can be ignored).   Σ a,b  means summation over both indices a and b; equivalent to Σ a  Σ b      We use lower case subscripts to indicate an element of a tensor e.g. u tkƒ     We denote sub tensors by dropping the appropriate subscript e.g. u tk  denotes the vector formed over ƒ from u tkƒ .
 
ML-ICA
       

     To provide some context we first outline Maximum Likelihood independent component analysis (ML-ICA): If we assume that the demixed output Y is drawn from a complex circular symmetric (CCS) Laplace distribution then we obtain
 
ρ( Y   tkƒ )∝ e   −|Y     tkƒ     | 
 
     Real audio signals tend to be heavy-tailed. The Laplace distribution is the most heavy-tailed distribution that retains the useful convergence property of being log-concave. Assuming independence, the log-likelihood of the observations X f  given the matrix W f  is then given by: 
                     L   ⁡     (       X   f     ;     W   f       )       ⁢     =   Δ     ⁢       2   ⁢   T   ⁢           ⁢   ln   ⁢          det   ⁢           ⁢     W   f              +       ∑     t   ,   k       ⁢     L   ⁡     (     Y   tkf     )                   (   5   )               
where L(Y tkƒ )           −|Y tkƒ |.

     For each frequency ƒ, ML-ICA searches over W ƒ  for a local maximum to eq(5). The result is an estimate for the sources up to a scaling ambiguity (B f ) and an inter-frequency permutation (P f ):
 
 Y   f   ≈S   f   B   f   P   f  
 
     ML-ICA then uses a separate permutation alignment operation to determine the inter-frequency permutations. One permutation alignment mechanism is to maximise the cross correlation of the output across frequencies according to some distance criterion. However one problem with this approach is that the fricatives and voiced parts of speech are normally anti-correlated, and this can lead to them being swapped between output channels. 
     NMF-ICA 
     Non-negative matrix factorisation (NMF) is a technique that can provide a good model for the structure inherent in real audio signals. The techniques we describe here combine NMF and ICA into a single unified approach where they can be jointly optimised. 
     We make the premise that the STFT time-frequency data is drawn from a statistical NMF-ICA model with unknown latent variables (which include the demixing matrices W ƒ ). The NMF-ICA algorithm then has four basic steps.
         Use the STFT to convert the time domain data into a time-frequency representation.   Use statistical inference to calculate either the maximum likelihood or the maximum posterior values for the latent variables. The algorithms work by iteratively improving an estimate for the latent variables.   Given estimates for the latent variables, the procedure can directly calculate the source estimates Y from eq(3).   Depending on the application the procedure can then either:
           use the inverse STFT to convert the estimate of Y back into the time domain; or   use a multi-channel inverse Fourier Transform on W to calculate the demixing time domain filter.
 
Maximum Likelihood NMF-ICA Model
   
               

     In deriving the NMF-ICA model we first express the probability of Y in terms of a generalisation of a complex normal distribution with unknown time-frequency dependent variance σ and: 
     
       
         
           
             
               
                 
                   
                     p 
                     ⁡ 
                     
                       ( 
                       
                         
                           Y 
                           tkf 
                         
                         ; 
                         
                           σ 
                           tkf 
                         
                       
                       ) 
                     
                   
                   ∝ 
                   
                     
                       1 
                       
                         σ 
                         tkf 
                         2 
                       
                     
                     ⁢ 
                     
                       
                         ⅇ 
                         
                           - 
                           
                             
                               
                                  
                                 
                                   Y 
                                   tkf 
                                 
                                  
                               
                               λ 
                             
                             
                               σ 
                               tkf 
                               λ 
                             
                           
                         
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   6 
                   ) 
                 
               
             
           
         
       
     
     The variances σ tkƒ  are then inferred from an NMF model with L non-negative components defined by the latent variables U, V (a set of dictionaries and activations for each source) such that 
     
       
         
           
             
               
                 
                   
                     σ 
                     tkf 
                     λ 
                   
                   = 
                   
                     
                       ∑ 
                       l 
                     
                     ⁢ 
                     
                       
                         V 
                         ltk 
                       
                       ⁢ 
                       
                         
                           U 
                           lfk 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
           
         
       
     
     We factorise σ tkf   λ  as it gives analytically tractable update equations. Assuming independence, we can write the overall log likelihood of the observations given the model parameters as: 
     
       
         
           
             
               
                 
                   
                     L 
                     ⁡ 
                     
                       ( 
                       
                         
                           Y 
                           tkf 
                         
                         ; 
                         
                           σ 
                           tkf 
                         
                       
                       ) 
                     
                   
                   ⁢ 
                   
                     = 
                     Δ 
                   
                   ⁢ 
                   
                     
                       
                         - 
                         
                           2 
                           λ 
                         
                       
                       ⁢ 
                       
                         ln 
                         ⁡ 
                         
                           ( 
                           
                             σ 
                             tkf 
                             λ 
                           
                           ) 
                         
                       
                     
                     - 
                     
                       
                         
                            
                           
                             Y 
                             tkf 
                           
                            
                         
                         λ 
                       
                       
                         σ 
                         tkf 
                         λ 
                       
                     
                   
                 
               
               
                 
                   ( 
                   8 
                   ) 
                 
               
             
           
         
       
     
     
       
         
           
             
               
                 
                   
                     L 
                     ⁡ 
                     
                       ( 
                       
                         
                           X 
                           ; 
                           W 
                         
                         , 
                         U 
                         , 
                         V 
                       
                       ) 
                     
                   
                   ⁢ 
                   
                     = 
                     Δ 
                   
                   ⁢ 
                   
                     
                       
                         ∑ 
                         f 
                       
                       ⁢ 
                       
                         2 
                         ⁢ 
                         T 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         ln 
                         ⁢ 
                         
                            
                           
                             det 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               W 
                               f 
                             
                           
                            
                         
                       
                     
                     + 
                     
                       
                         ∑ 
                         
                           t 
                           , 
                           k 
                           , 
                           f 
                         
                       
                       ⁢ 
                       
                         L 
                         ⁡ 
                         
                           ( 
                           
                             
                               Y 
                               tkf 
                             
                             ; 
                             
                               σ 
                               tkf 
                             
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
     The task is then to search over W, U, V for the maximum likelihood (ML) solution to equation (9). (The factors of 2 in eq. (8) and (9) are due to using complex circular symmetric distributions, although the NMF-ICA algorithm is robust to using a different factor). 
     This NMF-ICA model then has several advantages over ML-ICA:
         It unifies permutation alignment and ICA.   Taking λ&lt;1 will create a more heavy-tailed distribution at the expense of potentially introducing more local maxima.   Having several components allows uncorrelated and anti-correlated behaviour such as fricatives vs voiced behaviour to be modelled.   The latent variables U, V provide a wider solution space which will generally contain better solutions than ML-ICA
 
Related Models: Maximum a Posteriori Model
       

     One can introduce prior information about the latent variables using Bayes rule. Embodiments of this procedure use inverse gamma priors for U and V as they again lead to analytically tractable solutions. 
     
       
         
           
             
               L 
               ⁡ 
               
                 ( 
                 
                   
                     
                       V 
                       ltk 
                     
                     ; 
                     γ 
                   
                   , 
                   α 
                 
                 ) 
               
             
             ⁢ 
             
               = 
               Δ 
             
             ⁢ 
             
               
                 
                   - 
                   
                     ( 
                     
                       α 
                       + 
                       1 
                     
                     ) 
                   
                 
                 ⁢ 
                 ln 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   V 
                   ltk 
                 
               
               - 
               
                 γ 
                 
                   V 
                   ltk 
                 
               
             
           
         
       
       
         
           
             
               L 
               ⁡ 
               
                 ( 
                 
                   
                     
                       U 
                       lfk 
                     
                     ; 
                     
                       γ 
                       ′ 
                     
                   
                   , 
                   
                     α 
                     ′ 
                   
                 
                 ) 
               
             
             ⁢ 
             
               = 
               Δ 
             
             ⁢ 
             
               
                 
                   - 
                   
                     ( 
                     
                       
                         α 
                         ′ 
                       
                       + 
                       1 
                     
                     ) 
                   
                 
                 ⁢ 
                 ln 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   U 
                   lfk 
                 
               
               - 
               
                 
                   γ 
                   ′ 
                 
                 
                   U 
                   lfk 
                 
               
             
           
         
       
       
         
           
             
               L 
               ⁡ 
               
                 ( 
                 
                   
                     V 
                     ; 
                     γ 
                   
                   , 
                   α 
                 
                 ) 
               
             
             = 
             
               
                 ∑ 
                 
                   l 
                   , 
                   t 
                   , 
                   k 
                 
               
               ⁢ 
               
                 L 
                 ⁡ 
                 
                   ( 
                   
                     
                       
                         V 
                         ltk 
                       
                       ; 
                       γ 
                     
                     , 
                     α 
                   
                   ) 
                 
               
             
           
         
       
       
         
           
             
               L 
               ⁡ 
               
                 ( 
                 
                   
                     U 
                     ; 
                     
                       γ 
                       ′ 
                     
                   
                   , 
                   
                     α 
                     ′ 
                   
                 
                 ) 
               
             
             = 
             
               
                 ∑ 
                 
                   l 
                   , 
                   f 
                   , 
                   k 
                 
               
               ⁢ 
               
                 L 
                 ⁡ 
                 
                   ( 
                   
                     
                       
                         U 
                         lfk 
                       
                       ; 
                       
                         γ 
                         ′ 
                       
                     
                     , 
                     
                       α 
                       ′ 
                     
                   
                   ) 
                 
               
             
           
         
       
     
     Note that a side effect of the priors is to resolve the scaling ambiguity in the NMF model between U and V. This ambiguity does not matter from a theoretical point of view, but it can potentially cause numerical instability in practice. 
     Combining the priors with the observation likelihoods eq(9) using Bayes rule gives a posterior likelihood of
 
 L ( W,U,V;X , . . . )             L ( X;W ,σ,λ)+ L ( V ;γ,α)+ L ( U ;γ′,α′).  (10)

     Maximising this equation gives maximum a posteriori (MAP) solution to the problem. 
     Fixed Dictionary 
     Rather than learning U blindly from the data, embodiments of the procedure can use a fixed set of dictionaries, learnt from some representative single source training data. The scaling ambiguities between U, V and W mean that some of the variability that would have been captured by optimising U can be absorbed in the updates of V and W. A fixed dictionary can be a computational saving. 
     Input Channel Noise 
     Embodiments of the procedure can model stationary input channel noise by including an extra noise term in eq(7). This component has activations set to 1 and a fixed spectral dictionary. 
     Optimisation 
     The maximum a posteriori estimation (MAP) is found by maximising eq. (10). Similarly the maximum likelihood estimator (ML) is found by maximising eq. (9). Both these procedures are very similar, so we will derive the MAP estimator first. 
     We iteratively optimise L(W, U, V; X, . . . ) with respect to U, V and W. To optimise with respect to U and V we use a minorisation-maximisation (MM) algorithm. We optimise W using two different algorithms; the first optimises W with respect to permutations and output gain, the second uses a natural gradient method. All of these methods apart from the natural gradient give guaranteed convergence to a local maximum. The natural gradient method is expected to converge for a suitably small step size. 
     Optimisation with Respect to U 
     Looking at the terms in eq(10) that depend upon U one obtains: 
     
       
         
           
             
               
                 
                   
                     L 
                     ⁡ 
                     
                       ( 
                       
                         W 
                         , 
                         U 
                         , 
                         
                           V 
                           ; 
                           X 
                         
                         , 
                         … 
                       
                       ) 
                     
                   
                   ⁢ 
                   
                     = 
                     Δ 
                   
                   ⁢ 
                   
                     
                       - 
                       
                         
                           ∑ 
                           
                             t 
                             , 
                             k 
                             , 
                             f 
                           
                         
                         ⁢ 
                         
                           ( 
                           
                             
                               
                                 2 
                                 λ 
                               
                               ⁢ 
                               
                                 ln 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     σ 
                                     tkf 
                                     λ 
                                   
                                   ) 
                                 
                               
                             
                             + 
                             
                               
                                 
                                    
                                   
                                     Y 
                                     tkf 
                                   
                                    
                                 
                                 λ 
                               
                               
                                 σ 
                                 tkf 
                                 λ 
                               
                             
                           
                           ) 
                         
                       
                     
                     - 
                     
                       
                         ∑ 
                         
                           l 
                           , 
                           f 
                           , 
                           k 
                         
                       
                       ⁢ 
                       
                         ( 
                         
                           
                             
                               ( 
                               
                                 
                                   α 
                                   ′ 
                                 
                                 + 
                                 1 
                               
                               ) 
                             
                             ⁢ 
                             ln 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               U 
                               lfk 
                             
                           
                           + 
                           
                             
                               γ 
                               ′ 
                             
                             
                               U 
                               lfk 
                             
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
           
         
       
     
     If we take a hypothetical function ƒ(x), the first stage of MM is minorisation, which is creating an auxiliary function ƒ + ({circumflex over (x)},x), that satisfies
 
ƒ + ( {circumflex over (x)},x )≦ƒ({circumflex over ( x )})
 
ƒ + ( x,x )=ƒ( x )
 
     The auxiliary function should be one that is easier to maximise than the original. The second stage is then maximising f + ({circumflex over (x)},x) with respect to {circumflex over (x)}. Because of the constraints we know that ƒ({circumflex over (x)})≧ƒ + ({circumflex over (x)},x)≧ƒ(x), which proves that iterating the process will converge on a maximum of ƒ. One important property is that the sum of functions can be minorised by the sum of the individual minorisations. 
     In our case we have four terms. Two of the terms involve −ln x which is convex and can be minorised by simple linearisation about the current point 
     
       
         
           
             
               f 
               ⁡ 
               
                 ( 
                 x 
                 ) 
               
             
             = 
             
               
                 - 
                 ln 
               
               ⁢ 
               
                   
               
               ⁢ 
               x 
             
           
         
       
       
         
           
             
               
                 f 
                 + 
               
               ⁡ 
               
                 ( 
                 
                   
                     x 
                     ^ 
                   
                   , 
                   x 
                 
                 ) 
               
             
             = 
             
               - 
               
                 ( 
                 
                   
                     ln 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     x 
                   
                   + 
                   
                     
                       x 
                       ^ 
                     
                     x 
                   
                   - 
                   1 
                 
                 ) 
               
             
           
         
       
     
     The second term is 
             -       1       ∑   l     ⁢       U   lfk     ⁢     V   ltk           .           
A suitable minorisation for this is:
 
     
       
         
           
             
               σ 
               ftk 
               λ 
             
             = 
             
               
                 ∑ 
                 1 
               
               ⁢ 
               
                 
                   V 
                   ltk 
                 
                 ⁢ 
                 
                   U 
                   lfk 
                 
               
             
           
         
       
       
         
           
             
               f 
               ⁡ 
               
                 ( 
                 U 
                 ) 
               
             
             = 
             
               - 
               
                 1 
                 
                   σ 
                   ftk 
                   λ 
                 
               
             
           
         
       
       
         
           
             
               
                 f 
                 + 
               
               ⁡ 
               
                 ( 
                 
                   
                     U 
                     ^ 
                   
                   , 
                   U 
                 
                 ) 
               
             
             = 
             
               
                 - 
                 
                   1 
                   
                     
                       ( 
                       
                         σ 
                         ftk 
                         λ 
                       
                       ) 
                     
                     2 
                   
                 
               
               ⁢ 
               
                 
                   ∑ 
                   l 
                 
                 ⁢ 
                 
                   
                     
                       U 
                       lfk 
                       2 
                     
                     ⁢ 
                     
                       V 
                       ltk 
                     
                   
                   
                     
                       U 
                       ^ 
                     
                     lfk 
                   
                 
               
             
           
         
       
     
     Lastly we have terms of the form −1/x. These don&#39;t require minorising so we define 
     
       
         
           
             
               f 
               ⁡ 
               
                 ( 
                 x 
                 ) 
               
             
             = 
             
               - 
               
                 1 
                 x 
               
             
           
         
       
       
         
           
             
               
                 f 
                 + 
               
               ⁡ 
               
                 ( 
                 
                   
                     x 
                     ^ 
                   
                   , 
                   x 
                 
                 ) 
               
             
             = 
             
               - 
               
                 1 
                 
                   x 
                   ^ 
                 
               
             
           
         
       
     
     Putting these together gives a minorisation for eq(10) as: 
     
       
         
           
             
               
                 ?? 
                 + 
               
               ⁡ 
               
                 ( 
                 
                   
                     U 
                     ^ 
                   
                   , 
                   U 
                 
                 ) 
               
             
             = 
             
               
                 - 
                 
                   
                     ∑ 
                     
                       t 
                       , 
                       k 
                       , 
                       f 
                     
                   
                   ⁢ 
                   
                     
                       2 
                       λ 
                     
                     ⁢ 
                     
                       ( 
                       
                         
                           ln 
                           ⁡ 
                           
                             ( 
                             
                               σ 
                               tkf 
                               λ 
                             
                             ) 
                           
                         
                         + 
                         
                           
                             
                               ∑ 
                               l 
                             
                             ⁢ 
                             
                               
                                 
                                   U 
                                   ^ 
                                 
                                 lfk 
                               
                               ⁢ 
                               
                                 V 
                                 ltk 
                               
                             
                           
                           
                             σ 
                             tkf 
                             λ 
                           
                         
                         - 
                         1 
                       
                       ) 
                     
                   
                 
               
               + 
               
                 
                   
                     
                        
                       
                         Y 
                         tkf 
                       
                        
                     
                     λ 
                   
                   
                     σ 
                     tkf 
                     λ 
                   
                 
                 ⁢ 
                 
                   
                     ∑ 
                     l 
                   
                   ⁢ 
                   
                     
                       
                         U 
                         lfk 
                         2 
                       
                       ⁢ 
                       
                         V 
                         ltk 
                       
                     
                     
                       
                         U 
                         ^ 
                       
                       lfk 
                     
                   
                 
               
               - 
               
                 
                   ∑ 
                   
                     l 
                     , 
                     f 
                     , 
                     k 
                   
                 
                 ⁢ 
                 
                   
                     ( 
                     
                       
                         α 
                         ′ 
                       
                       + 
                       1 
                     
                     ) 
                   
                   ⁢ 
                   
                     ( 
                     
                       
                         ln 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           U 
                           lfk 
                         
                       
                       + 
                       
                         
                           
                             U 
                             ^ 
                           
                           lfk 
                         
                         
                           U 
                           lfk 
                         
                       
                       - 
                       1 
                     
                     ) 
                   
                 
               
               + 
               
                 
                   γ 
                   ′ 
                 
                 
                   
                     U 
                     ^ 
                   
                   lfk 
                 
               
             
           
         
       
     
     This auxiliary function only has a single maximum, which can be found analytically. Solving for Û gives: 
     
       
         
           
             
               
                 
                   
                     
                       
                         ∂ 
                         
                           ?? 
                           + 
                         
                       
                       
                         ∂ 
                         
                           
                             U 
                             ^ 
                           
                           lfk 
                         
                       
                     
                     = 
                     
                       
                         - 
                         
                           
                             
                               α 
                               ′ 
                             
                             + 
                             1 
                           
                           
                             U 
                             lfk 
                           
                         
                       
                       + 
                       
                         
                           γ 
                           ′ 
                         
                         
                           
                             U 
                             ^ 
                           
                           lfk 
                           2 
                         
                       
                       - 
                       
                         
                           ∑ 
                           t 
                         
                         ⁢ 
                         
                           ( 
                           
                             
                               
                                 2 
                                 λ 
                               
                               ⁢ 
                               
                                 V 
                                 ltk 
                               
                               ⁢ 
                               
                                 
                                   σ 
                                   ^ 
                                 
                                 tkf 
                                 
                                   - 
                                   λ 
                                 
                               
                             
                             - 
                             
                               
                                 
                                   
                                      
                                     
                                       Y 
                                       tkf 
                                     
                                      
                                   
                                   λ 
                                 
                                 
                                   σ 
                                   ftk 
                                   
                                     2 
                                     ⁢ 
                                     λ 
                                   
                                 
                               
                               ⁢ 
                               
                                 
                                   
                                     U 
                                     lfk 
                                     2 
                                   
                                   ⁢ 
                                   
                                     V 
                                     ltk 
                                   
                                 
                                 
                                   
                                     U 
                                     ^ 
                                   
                                   lfk 
                                   2 
                                 
                               
                             
                           
                           ) 
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     
                       
                         U 
                         ^ 
                       
                       lfk 
                     
                     = 
                     
                       
                         
                           
                             γ 
                             ′ 
                           
                           + 
                           
                             
                               U 
                               lfk 
                               2 
                             
                             ⁢ 
                             
                               
                                 ∑ 
                                 t 
                               
                               ⁢ 
                               
                                 
                                   V 
                                   ltk 
                                 
                                 ⁢ 
                                 
                                   
                                      
                                     
                                       Y 
                                       tkf 
                                     
                                      
                                   
                                   λ 
                                 
                                 ⁢ 
                                 
                                   σ 
                                   tkf 
                                   
                                     
                                       - 
                                       2 
                                     
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     λ 
                                   
                                 
                               
                             
                           
                         
                         
                           
                             
                               
                                 α 
                                 ′ 
                               
                               + 
                               1 
                             
                             
                               U 
                               lfk 
                             
                           
                           + 
                           
                             
                               2 
                               λ 
                             
                             ⁢ 
                             
                               
                                 ∑ 
                                 t 
                               
                               ⁢ 
                               
                                 
                                   V 
                                   ltk 
                                 
                                 ⁢ 
                                 
                                   σ 
                                   tkf 
                                   
                                     - 
                                     λ 
                                   
                                 
                               
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
           
         
       
     
     Importantly this update is guaranteed to improve the likelihood eq(10), so it can be interleaved with the other stages of the NMF-ICA algorithm. Having calculated Û lfk  for all l, ƒ, k, we can then update U by assigning U lƒk ←Û lƒk . 
     Optimisation with Respect to V 
     Optimising with respect to V follows the same process as for U. By symmetry we obtain: 
     
       
         
           
             
               
                 
                   
                     
                       V 
                       ^ 
                     
                     ltk 
                   
                   = 
                   
                     
                       
                         γ 
                         + 
                         
                           
                             V 
                             ltk 
                             2 
                           
                           ⁢ 
                           
                             
                               ∑ 
                               f 
                             
                             ⁢ 
                             
                               
                                 U 
                                 lfk 
                               
                               ⁢ 
                               
                                 
                                    
                                   
                                     Y 
                                     tkf 
                                   
                                    
                                 
                                 λ 
                               
                               ⁢ 
                               
                                 σ 
                                 tkf 
                                 
                                   
                                     - 
                                     2 
                                   
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   λ 
                                 
                               
                             
                           
                         
                       
                       
                         
                           
                             
                               α 
                               ′ 
                             
                             + 
                             1 
                           
                           
                             V 
                             ltk 
                           
                         
                         + 
                         
                           
                             2 
                             λ 
                           
                           ⁢ 
                           
                             
                               ∑ 
                               f 
                             
                             ⁢ 
                             
                               
                                 U 
                                 lfk 
                               
                               ⁢ 
                               
                                 σ 
                                 tkf 
                                 
                                   - 
                                   λ 
                                 
                               
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
     Having calculated {circumflex over (V)} ltk  for all l, t, k, we can then update V by assigning V ltk ←{circumflex over (V)} ltk . 
     Rank 1 Solution to U and V 
     When L=1 it is possible to solve 
                 ∂   L       ∂     V   tk         =   0         
directly without minorisation-minimisation. (Note that we can drop the l subscript). The solution is given by
 
     
       
         
           
             
               
                 
                   
                     
                       V 
                       ^ 
                     
                     tk 
                   
                   = 
                   
                     
                       γ 
                       + 
                       
                         
                           ∑ 
                           f 
                         
                         ⁢ 
                         
                           
                             
                                
                               
                                 Y 
                                 tkf 
                               
                                
                             
                             λ 
                           
                           
                             U 
                             fk 
                           
                         
                       
                     
                     
                       α 
                       + 
                       1 
                       + 
                       
                         
                           2 
                           ⁢ 
                           F 
                         
                         λ 
                       
                     
                   
                 
               
               
                 
                   ( 
                   14 
                   ) 
                 
               
             
           
         
       
     
     Similarly the solution for 
                 ∂   L       ∂     U   fk         =   0         
is
 
     
       
         
           
             
               
                 
                   
                     
                       U 
                       ^ 
                     
                     fk 
                   
                   = 
                   
                     
                       
                         γ 
                         ′ 
                       
                       + 
                       
                         
                           E 
                           t 
                         
                         ⁢ 
                         
                           
                             
                                
                               
                                 Y 
                                 tkf 
                               
                                
                             
                             λ 
                           
                           
                             V 
                             tk 
                           
                         
                       
                     
                     
                       
                         α 
                         ′ 
                       
                       + 
                       1 
                       + 
                       
                         
                           2 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           T 
                         
                         λ 
                       
                     
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
       
     
     The rank-1 solution for U ƒk  is redundant as, without any loss of generality, it can be absorbed into the scaling values Λ kkƒ  (described later). 
     Optimisation with Respect to W 
     In optimising W we introduce the following matrix notation:
         σ ƒ  is the appropriate T×K matrix (T-rows; K-columns) formed from σ tkƒ .   Tr A is the trace of matrix A.   I is the identity matrix   Element wise operations are indicated by • as follows:
           A•B is element wise (Hadamard) multiplication,   ln•A take the natural logarithm the elements of A,   A •λ  raises the elements of A to the power λ,   abs •  A takes the absolute values of the elements of A.   
               

     Optimising W is independent of the priors on U and V, so we can rewrite both the MAP and ML likelihood equations as functions of W, σ plus a constant offset as
 
 L ( W   ƒ ;σ •λ , . . . )             T  ln  det ( W   ƒ   H   W   ƒ )− Tr ((σ ƒ   •λ ) T   abs   •   Y   ƒ   •λ )  (16)

     It is also useful to define an intermediate variable 
               Q   f     =       λ     2   ⁢   T       ⁢       (     σ   f     •   -   λ       )     T     ⁢     abs   ·     Y   f   •λ               
Permutation and Scaling
 
     The likelihood equation (16) can be directly optimised with respect to permutations and scaling. Let P ƒ  be a permutation matrix and Λ f  be a diagonal real scaling matrix. The updated value of W ƒ  will be given by
 
 W   ƒ   ←W   ƒ   P   ƒ Λ ƒ .  (17)
 
     Note that for permutation and scaling we have
 
| detP   ƒ |=1
 
 abs   • ( Y   ƒ   P   ƒ Λ ƒ ) •λ =( abs   •   Y   ƒ   •λ ) P   ƒ Λ ƒ   λ .
 
     Treating the likelihood as a function of P ƒ  and Λ ƒ  we get 
     
       
         
           
             
               
                 
                   
                     L 
                     ⁡ 
                     
                       ( 
                       
                         
                           
                             
                               W 
                               f 
                             
                             ⁢ 
                             
                               P 
                               f 
                             
                             ⁢ 
                             
                               Λ 
                               f 
                             
                           
                           ; 
                           
                             σ 
                             f 
                             •λ 
                           
                         
                         , 
                         … 
                       
                       ) 
                     
                   
                   ⁢ 
                   
                     = 
                     Δ 
                   
                   ⁢ 
                   
                     
                       
                         2 
                         ⁢ 
                         T 
                       
                       λ 
                     
                     ⁢ 
                     
                       
                         ( 
                         
                           
                             ln 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               
                                 det 
                                 ⁢ 
                                 Λ 
                               
                               f 
                               λ 
                             
                           
                           - 
                           
                             Tr 
                             ⁡ 
                             
                               ( 
                               
                                 
                                   Q 
                                   f 
                                 
                                 ⁢ 
                                 
                                   P 
                                   f 
                                 
                                 ⁢ 
                                 
                                   Λ 
                                   f 
                                   λ 
                                 
                               
                               ) 
                             
                           
                         
                         ) 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   18 
                   ) 
                 
               
             
           
         
       
     
     For each ƒ we can now maximise L with respect to Λ ƒ  given P ƒ  to show that 
     
       
         
           
             
               
                 
                   
                     Λ 
                     kkf 
                   
                   = 
                   
                     
                       ( 
                       
                         
                           Q 
                           f 
                         
                         ⁢ 
                         
                           P 
                           f 
                         
                       
                       ) 
                     
                     kk 
                     
                       
                         - 
                         1 
                       
                       / 
                       λ 
                     
                   
                 
               
               
                 
                   ( 
                   19 
                   ) 
                 
               
             
           
         
       
     
     If we substitute eq(19) back into the eq(18) we can solve for P ƒ  by
 
 P   ƒ   ←               Tr ( P   ƒ ln •   Q   ƒ )  (20)

     To update W ƒ  we therefore calculate Q f  as defined above, then apply equation (20), equation (19) and finally equation (17). 
     Using this permutation and scaling stage can alleviate the local maxima problem and allows the procedure to use λ&lt;1, as it can jump the solution between local maxima. 
     Natural Gradient 
     Eq. (16) can be differentiated with respect to W using Wirtinger calculus to give
 
∇ L   ƒ =2 T·W   ƒ   −H −λ((σ ƒ   •-λ   •abs   •   Y   ƒ   •λ   •Y   ƒ   •-1 ) T   X   ƒ ) H   (21)
 
     The natural gradient ∂W ƒ  is the direction of steepest ascent of L with respect to distance ds travelled in a Riemannian manifold. The manifold for invertible matrices gives us
 
 ds   2   =∥W   ƒ   −1   ∂W   ƒ ∥ F   2  
 
∂ W   ƒ   =W   ƒ   W   ƒ   H   ∇L   ƒ .  (22)
 
     Using an intermediate variable Ψ ƒ  and a step size μ we can substitute eq(21) into eq(22) to get an update equation 
                       Ψ   f     =       λ     2   ⁢           ⁢   T       ⁢       (       σ   f     •   -   λ       ·   abs   ·     Y   f   •λ     ·     Y   f     •   -   1         )     T     ⁢     Y   f         ⁢     
     ⁢       ∂     W   f       =       W   f     ⁡     (     I   -     Ψ   f   H       )                 (   23   )                 W   ƒ   ←W   ƒ   +μ∂w   ƒ   (24)
 
     In the above update equations the superscript T refers to transpose, and the variable T to the number of frames. The calculation of Ψ ƒ  is deterministic, from Y ƒ  and σ ƒ ; the step size μ can be a fixed value such as 0.1. 
     Diagonal Scaled Natural Gradient 
     The procedure can skip the permutation but still perform the scaling efficiently as part of the natural gradient. Where the procedure skips the permutation we have P ƒ =I. Since Q kkƒ =Ψ kkƒ  (i.e. the diagonal elements of Q ƒ  and Ψ ƒ  are the same), we can calculate the diagonal scaling Λ ƒ  in eq(19) directly from Ψ ƒ . The scaling can then be incorporated into the update by making the following substitutions in eq(23) and eq(24): 
               Λ   kkf     ←     Ψ   kkf     -     1   λ                       Ψ   f     ←       Λ   f     ⁢     Ψ   f   H     ⁢     Λ   f     λ   -   1                       W   f     ←       W   f     ⁢     Λ   f             
Overall Blind Source Separation Algorithm
 
     A preferred embodiment of the overall algorithm recursively uses each of the above optimisation steps to improve the estimates of W, U, V. An example implementation is set out below, noting that the initialisation can change, and that in principle the update steps 2(a,b), 2(c,d), 2(e,f), 2(g,h) may be performed in any order: 
     1. Initialise W, U, V for example as follows:
         a. W ƒ =I for all ƒ.   b. U is randomly initialised with non-negative real values.   c. V is randomly initialised with non-negative real values.   d. Y ƒ ←X ƒ  for all ƒ.   e. σ k   •λ ←V k   T U k  for all k.       

     2. Repeat until convergence:
         a. Update W using the permutation and scaling eq(17) for all ƒ.   b. Y ƒ ←X ƒ  W ƒ  for all ƒ.   c. Update W using the natural gradient eq(24) for all ƒ.   d. Y ƒ ←X ƒ  W ƒ  for all ƒ.   e. Update U using minorisation-maximisation eq(12).   f. σ k   •λ ←V k   T U k  for all k.   g. Update V using minorisation-maximisation eq(13).   h. σ k   •λ ←V k   T U k  for all k.       

     The convergence criterion can be a fixed number of iterations (say 25 to 30), or until there has been no significant change in W, U or V. 
     A preferred embodiment of employs random initialisation of U and V so that each component is initialised with a different profile. Alternatively initialisations from priors or from the data may be employed. 
     In broad terms, embodiments of the procedure aim to maximise eq(9) with respect to W, U and V, or eq(10) if there are priors on U, V. 
     Maximum Likelihood Variant 
     The maximum likelihood criterion is essentially the same as the MAP estimator, but without the priors on U or V. Thus in embodiments the only effect is a minor change to the updates on U and V. The update on U, eq(12), becomes: 
                 U   ^     lfk     =       U   lfk     ⁢           ∑   t     ⁢       V   ltk     ⁢            Y   tkf          λ     ⁢     σ   tkf       -   2     ⁢   λ               2   λ     ⁢       ∑   t     ⁢       V   ltk     ⁢     σ   tkf     -   λ                         
Similarly eq(13) becomes:
 
     
       
         
           
             
               
                 V 
                 ^ 
               
               ltk 
             
             = 
             
               
                 V 
                 ltk 
               
               ⁢ 
               
                 
                   
                     
                       ∑ 
                       t 
                     
                     ⁢ 
                     
                       
                         U 
                         lfk 
                       
                       ⁢ 
                       
                         
                            
                           
                             Y 
                             tkf 
                           
                            
                         
                         λ 
                       
                       ⁢ 
                       
                         σ 
                         tkf 
                         
                           
                             - 
                             2 
                           
                           ⁢ 
                           λ 
                         
                       
                     
                   
                   
                     
                       2 
                       λ 
                     
                     ⁢ 
                     
                       
                         ∑ 
                         f 
                       
                       ⁢ 
                       
                         
                           U 
                           lfk 
                         
                         ⁢ 
                         
                           σ 
                           tkf 
                           
                             - 
                             λ 
                           
                         
                       
                     
                   
                 
               
             
           
         
       
     
     There are also maximum likelihood equivalents to the rank-1 equations. 
     Extensions 
     The above described procedure performs NMF-ICA blind source separation. 
     However extensions are desirable for practical application of the techniques to listening aids and in other fields. 
     Dimension Reduction 
     Embodiments of the above described procedure demix the signals from K audio channels (microphones) into signals from K putative sources. However where the number of sources (eg 2) is less than the number of microphones (eg 8), sources can be fragmented—one real source can be split across two or more presumed sources. For this reason it can be beneficial to reduce the dimensionality of the input data (number of input channels) to match the actual number of sources. 
     Thus where the number of sources is less than the number of microphones the procedure can use dimension reduction to reduce the problem to a smaller number of virtual microphones. The microphone observations X ƒ  are pre-processed by a multichannel linear filter W DR     ƒ    which has fewer columns than rows:
 
 X   ƒ   ′=X   ƒ   W   DR     ƒ   
 
     It is these virtual microphone signals X ƒ ′ which are then passed to the NMF-ICA algorithm. For example, if K R  is the reduced number of channels, then W DR     ƒ    is a K by K R  matrix, X ƒ  is a T by K matrix, and X ƒ ′ is a T by K R  matrix. 
     The simplest form of dimension reduction is to discard microphones, but Principal Component Analysis gives a minimum distortion dimension reduction. It is found by setting W DR     ƒ    to the set of eigenvectors corresponding to the largest eigenvalues of X ƒ   H X ƒ . 
     Scaling Ambiguity 
     Embodiments of the above described procedure extract the source estimates up to an arbitrary diagonal scaling matrix B ƒ . This is an arbitrary filter, since there is a value of B ƒ  at each frequency (this can be appreciated from the consideration that changing the bass or treble would not affect the independence of the sources). There is an unknown filter arising from the transfer function of the room, but the arbitrary filter can be removed by considering what a source would have sounded like at a particular microphone. 
     In one approach the scaling ambiguity can be resolved by taking one source, undoing the effect of the demixing to see what it would have sounded like at one or more of the microphones, and then using the result to adjust the scaling of the demixing matrix to match what was actually received (heard)—that is, applying a minimum distortion principle. This correction can be subsumed into a modified demixing matrix. 
     The procedure can estimate the sources as received at the microphones using a minimum distortion principle as follows:
         Let Ŵ ƒ  be the combined demixing filter including any dimension reduction or other pre processing e.g. Ŵ ƒ =W DR     ƒ    W ƒ .   Let Ŵ ƒ   †  be the pseudo inverse of Ŵ ƒ . This is a minimum distortion projection from the source estimates back to the microphones.   Let D(k) be a selector matrix which is zero everywhere except for one element on the diagonal D(k) kk =1.       

     To project source estimate k back to all the microphones we use
 
 W   ƒ ′( k )= Ŵ   ƒ   D ( k ) Ŵ   ƒ   †   (25)
 
 Ŷ   ƒ ( k )= X   ƒ   W   ƒ ′( k )  (26)
 
     Matrix D(k) selects one source k, and equations (25) and (26) define an estimate for the selected source on all the microphones. In equation (26) Ŷ ƒ  (k) is an estimate of how the selected source would have sounded at microphones, rather than an estimate of the source itself, because the (unknown) room transfer function is still present. 
     Source Selection 
     Oftentimes it is only a subset of the sources that is desired. Because there is still a global permutation, it may be useful to estimate which of the sources are the desired ones—that is, the sources have been separated into independent components but there is still ambiguity as to which source is which (eg in the case of a group of speakers around a microphone, which source k is which speaker). In addition embodiments of the procedure operate on time slices of the audio (successive groups of STFT frames) and it is not guaranteed that the “physical” source labelled as, say, k=1 in one group of frames will be the same “physical” source as the source labelled as k=1 in the next group of frames (this depends upon the initialisation of U, V, and W, which may, for example, be random or based on a previous group of frames). 
     Source selection may be made in various ways, for example on the basis of voice (or other audio) identification, or matching a user selected direction. Other procedures for selecting a source include selecting the loudest source (which may comprise selecting a direction from which there is most power); and selecting based upon a fixed (predetermined) direction for the application. For example the wanted source may be a speaker with a known direction with respect to the microphones. A still further approach is to look for a filter selecting a particular acoustic source which is similar to a filter in an adjacent time-frequency block, assuming that similar filters correspond to the same source. Such approaches enable a consistent global permutation matrix (P) to be determined from one time-frequency block to another. 
     In embodiments to match a user-selected direction knowledge of the expected microphone phase response θ jf  from the indicated direction may be employed. This can either be measured or derived from a simple anechoic model given the microphone geometry relative to an arbitrary origin. A simple model of the response of microphone j may be constructed as follows: 
     Given the known geometry for each microphone we can define
         s is the speed of sound.     x   j  is the position of microphone j relative to an arbitrary origin in real space     d  is a unit vector corresponding to a chosen direction towards the desired source in the same coordinate system as  x   j .   ρ is the sample rate (of digitised samples from the microphone).       

     The far field microphone time delay, τ j , in samples relative to the origin is then given by 
     
       
         
           
             
               τ 
               j 
             
             = 
             
               - 
               
                 
                   ρ 
                   ⁢ 
                   
                     
                       x 
                       _ 
                     
                     j 
                     T 
                   
                   ⁢ 
                   
                     d 
                     _ 
                   
                 
                 s 
               
             
           
         
       
     
     This leads to a phase shift for microphone j of 
     
       
         
           
             
               θ 
               jf 
             
             = 
             
               ⅇ 
               
                 2 
                 ⁢ 
                 
                   πτ 
                   j 
                 
                 ⁢ 
                 
                   f 
                   F 
                 
                 ⁢ 
                 ⅈ 
               
             
           
         
       
     
     However the phase response θ jf  is determined, the chosen source k s  is the source whose corresponding row in Ŵ ƒ   554  maximises the phase correlation: 
               k   s     =     arg   ⁢               ⁢   max     k     ⁢       ∑   f     ⁢              ∑   j     ⁢           W   ^     kjf   †              W   ^     kjf   †            ⁢     θ   jf   *              2               
where the sum j runs over the microphones and θ jf  is the (complex) frequency/phase response of microphone j in the selected direction. In principle this approach could be employed to select multiple source directions.
 
Low Latency Implementation
 
     In embodiments of the above described procedure the output of the procedure may be Y ƒ  or Ŷ ƒ  (k); additionally or alternatively an output may be the demixing filter W ƒ  or W ƒ ′(k). Where the output comprises a demixing filter this may be provided in the time-frequency domain or converted back into the time domain (as used in eq(1) above) and applied to the time domain data x t . Where filtering is performed in the time domain the time domain data may be delayed so that the filtering is applied to the time domain data from which it was derived, or (as the calculations can be relatively time-consuming), the filtering may be applied to the current time domain data (thus using coefficients which are slightly delayed relative to the data) 
     In some real-time applications, such as a listening aid, low latency is desirable. In this case, the filtering may be performed in the time domain using eq(1). The filter coefficients w are updated by using eq(25) to design the filter coefficients asynchronously in the STFT domain. For example, if calculation of the filter coefficients can be performed, say, every second then the coefficients are around 1 second out of date. This presents no problem if the acoustic scene is reasonably static (the speakers do not move around much), so that the filter coefficients are appropriate for later samples. If low latency is not needed, the procedure can use an inverse STFT on eq(26). 
     Stereo Filtering 
     A stereo output signal can be created by selecting an appropriate pair of rows from W ƒ ′ in eq(25). This leads to a more natural sounding output which still retains some of the spatial cues from the source. A listener who has not lost too much of their ability to spatially discriminate sounds can make use of these cues to further aid in discrimination against any residual interference. 
     Example Implementations 
     Referring to  FIG. 2 , this shows the architecture of apparatus  200  to improve the audibility of an audio signal by blind source separation, employing time-domain filtering to provide low latency. The apparatus comprises a microphone array  202  with microphones  202   a - n , coupled to a multi-channel analogue-to-digital converter  204 . This provides a digitised multi-channel audio output  205  to a spatial filter  206  which may be implemented as a multi-channel linear convolutional filter, and to a filter coefficient determiner  208 . The filter coefficient determiner  208  determines coefficients of a demixing filter which are applied by spatial filter  206  to extract audio from one (or more) selected sources for a demixed audio output  210 . The filter determiner  208  accepts optional user input, for example to select a source, and has an output  212  comprising demixing filter coefficients for the selected source. The demixed audio  210  is provided to a digital-to-analogue converter  214  which provides a time domain audio output  216 , for example to headphones or the like, or for storage/further processing (for example speech recognition), communication (for example over a wired or wireless network such as a mobile phone network and/or the Internet), or other uses. In  FIG. 2  the audio signal path is shown in bold. 
     In embodiments it is assumed that the acoustic scene is quasi-static and thus the filter coefficient determiner  208  and spatial filter  206  can operate in parallel. The latency is then determined by the main acoustic path (shown in bold), and depends upon the group delay of the filter coefficients, the latency of the spatial filter implementation, and the input/output transmission delays. Many different types of spatial filter may be used—for example one low latency filter implementation is to use a direct convolution; a more computationally efficient alternative is described in Gardener, W G (1995), “Efficient Convolution without Input-output Delay”,  Journal of the Audio Engineering Society′,  43 (3), 127-136. 
     The skilled person will recognise that the signal processing illustrated in the architecture of  FIG. 2  may be implemented in many different ways. For example the filter designer, in preferred embodiments with a user interface, and/or spatial filter and/or DAC  214  may be implemented on a general purpose computing device such as a mobile phone, tablet, laptop or personal computer. In embodiments the microphone array and ADC  204  may comprise part of such a general purpose computing device. Alternatively some or all of the architecture of  FIG. 2  may be implemented on a dedicated device such as dedicated hardware (for example an ASIC), and/or using a digital signal processor (DSP). A dedicated approach may reduce the latency on the main acoustic path which is otherwise associated with input/output to/from a general purpose computing device, but this may be traded against the convenience of use of a general purpose device. 
     An example spatial filter  206  for the apparatus of  FIG. 2  is shown in  FIG. 3 a   . The illustrated example shows a multi-channel linear discrete convolution filter in which the output is the sum of the audio input channels convolved with their respective filter co-efficients, as described in eq(1) above. In embodiments a multi-channel output such as a stereo output is provided. For a stereo output either the spatial filter output may be copied to all the output channels or more preferably, as shown in  FIG. 3 a   , a separate spatial filter is provided for each output channel. This latter approach is advantageous as it can approximate the source as heard by each ear (since the microphones are spaced apart from one another). This can lead to a more natural sounding output which still retains some spatial cues from the source. Thus a listener who has not lost too much of their ability to spatially discriminate sounds can employ those cues to further aid in discrimination against any residual interference. 
       FIG. 3 b    shows time-frequency and frequency-time domain conversions (not shown in  FIG. 2 ) for the frequency domain filter coefficient determiner  208  of  FIG. 2 . In embodiments each audio channel may be provided with an STFT (Short Time Fourier Transform) module  207   a - n  each configured to perform a succession of overlapping discrete Fourier transforms on an audio channel to generate a time sequence of spectra. Transformation of filter coefficients back into the time domain may be performed by a set of inverse discrete Fourier transforms  209 . 
     The Discrete Fourier Transform (DFT) is a method of transforming a block of data between a time domain representation and a frequency domain representation. The STFT is an invertible method where overlapping time domain frames are transformed using the DFT to a time-frequency domain. The STFT is used to apply the filtering in the time-frequency domain; in embodiments when processing each audio channel, each channel in a frame is transformed independently using a DFT. Optionally the spatial filtering could also be applied in the time-frequency domain, but this incurs a processing latency and thus more preferably the filter coefficients are determined in the time-frequency domain and then inverse transformed back into the time domain. The time domain convolution maps to frequency domain multiplication. 
     Referring now to  FIG. 4 , this shows modules of a preferred implementation of a frequency domain filter coefficient determiner  208  for use in embodiments of the invention. The modules of  FIG. 4  operate according to the procedure as previously described. Thus the filter coefficient determination system receives digitised audio data from the multiple audio channels in a time-frequency representation, from the STFT modules  207   a - n  of  FIG. 3 b   , defining the previously described observation matrix X f . This is provided to an optional dimension reduction module  402  which reduces the effective number of audio channels according to a dimension reduction matrix W DRf . The dimension reduction matrix, which in embodiments has fewer columns than rows, is determined (module  404 ) either in response to user input defining the number of sources to demix or in response to a determination by the system of the number of sources to demix, step  406 . The procedure may determine the number of sources based upon prior knowledge or, for example, on some heuristic measure of the output or, say, based on user feedback on the quality of demixed output. In a simple implementation the dimension reduction matrix may simply discard some of the audio input channels but in other approaches the input channels can be mapped to a reduced number of channels, for example using PCA as previously outlined. The complete or reduced set of audio channels is provided to a blind source separation module  410  which implements a procedure as previously described to perform joint NMF-ICA source separation. 
     The blind source separation module  410  provides a set of demixing matrices as an output, defining frequency domain filter coefficients W f . In embodiments these are provided to module  412  which removes the scaling ambiguity as previously described, providing filter coefficients for a source k at all the microphones (or reduced set of microphones). The user or the procedure then selects one or more of these microphones (by selecting data from one or more rows of W f (k)), which are then output for use by the spatial filter after conversion back into the time domain. 
     In embodiments a source selection module  416  operates on a pseudo inverse of the demixing matrix, using the microphone phase responses to choose a source k s . The source may be selected  418  either by the user, for example the user indicating a direction of the desired source, or by the procedure, for example based on a priori knowledge of the source direction. 
       FIG. 5  shows a flow diagram of a procedure for blind source separation according to an embodiment of the invention; this procedure may be used to implement the blind source separation module  410  and dimension reduction  402  of  FIG. 4 . Thus at step S 100  the procedure inputs audio data and then converts this to the time-frequency domain, optionally reducing the number of audio channels (S 102 ). The procedure also initialises latent variables U, V and the demixing matrices W, for example randomly or as previously outlined, and then calculates initial values for Y and σ. The procedure then repeats a number of update steps until convergence (S 106 ); the convergence criterion may be a fixed number of iterations. 
     Update step S 108  replaces W with a permuted and scaled version by calculating Q f  then Λ (eq19 above), then P f (eq20), using this to update W (eq17). Update step S 110  steps up the slope of W, performing a gradient search according to eq24. In an alternative approach step S 110  may recalculate the NMF model rather than updating the model. Update steps S 112 , S 114 , update the latent variables U, V using, for example, equations 12 and 13 or the maximum likelihood alternatives described above. 
     Once convergence has been achieved preferably the procedure resolves scaling ambiguity (S 114 ; implemented in module  412  of  FIG. 4 ), and optionally converts the filter coefficients back to the time domain (S 114 ). 
       FIG. 6  shows an example of a general purpose computing system  600  programmed to implement a system as described above to improve audibility of an audio signal by blind source separation according to an embodiment of the invention. Thus the computing system comprises a processor  602 , coupled to working memory  604 , program memory  606 , and to storage  608 , such as a hard disk. Program memory  606  comprises code to implement embodiments of the invention, for example operating system code, time to frequency domain conversion code, frequency to time domain conversion code, dimension reduction code, blind source separation code, scaling code, source selection code, and spatial (time domain) filter code. Working memory  604 /storage  608  stores data for the above-described variables W, U, V, X, Y, σ, and θ. Processor  602  is also coupled to a user interface  612 , to a network/communications interface  612 , and to an (analogue or digital) audio data input/output module  614 . The skilled person will recognise that audio module  614  is optional since the audio data may alternatively be obtained, for example, via network/communications interface  612  or from storage  608 . 
     Although in some preferred implementations the above described techniques are applied to audio comprising speech, the techniques are not limited to such applications and can be applied to other acoustic source separation problems, for example processing seismic data. Often a selected source comprises a human speaker to provide a listening aid or to assist teleconferencing, machine hearing or speech recognition, or other in applications such as selectively capturing speech from a driver or passenger in a vehicle for a vehicle phone. In some applications, however, embodiments of the techniques may be employed to identify a noise-like source (for example a source with the most noise-like characteristics may be selected), and this selected source may then be employed for active noise cancellation. 
     In principle the techniques we describe may be employed outside the audio/acoustic domain, for example to mixed-source electrical signal data such as data from sensing apparatus or instrumentation such as laboratory or medical apparatus. Examples include EEG (electroencephalography) data, and mixed source spectral data from a spectrum analyser such as an optical spectrum analyser, mass spectrum analyser or the like. 
     No doubt many other effective alternatives will occur to the skilled person. It will be understood that the invention is not limited to the described embodiments and encompasses modifications apparent to those skilled in the art lying within the spirit and scope of the claims appended hereto.