U.S. patent application number 13/097628 was filed with the patent office on 2012-11-01 for systems and methods for blind localization of correlated sources.
This patent application is currently assigned to Siemens Corporation. Invention is credited to Heiko Claussen, Justinian Rosca, Thomas Wiese.
Application Number | 20120275271 13/097628 |
Document ID | / |
Family ID | 46148952 |
Filed Date | 2012-11-01 |
United States Patent
Application |
20120275271 |
Kind Code |
A1 |
Claussen; Heiko ; et
al. |
November 1, 2012 |
SYSTEMS AND METHODS FOR BLIND LOCALIZATION OF CORRELATED
SOURCES
Abstract
A system and a method for a blind direction of arrival
estimation is provided for a nonlinear 1-dimensional array of M
receivers of J<M unknown signals transmitted through an unknown
channel and subject to Gaussian noise. The method overcomes the
conditions of that signals be disjoint in the time-frequency domain
and/or statistically independent. A likelihood function evaluates a
probability for the directions of arrivals. A weighting and
regularization scheme selects relevant frequency components. A
joint likelihood function of all frequencies is interpreted as a
cost or objective function. The time domain is included into the
estimation procedure by using a particle filter that operates on
the constructed likelihood function and a simple transition kernel
for the directions.
Inventors: |
Claussen; Heiko;
(Plainsboro, NJ) ; Rosca; Justinian; (West
Windsor, NJ) ; Wiese; Thomas; (Munich, DE) |
Assignee: |
Siemens Corporation
Iselin
NJ
|
Family ID: |
46148952 |
Appl. No.: |
13/097628 |
Filed: |
April 29, 2011 |
Current U.S.
Class: |
367/118 |
Current CPC
Class: |
G01S 3/8006 20130101;
G01S 3/74 20130101 |
Class at
Publication: |
367/118 |
International
Class: |
G01S 3/80 20060101
G01S003/80 |
Claims
1. A method for determining a direction of arrival of a signal
generated by a source in a plurality of sources, comprising: a
processor estimating a first direction of arrival for a previous
signal generated by the source with a particle filter operating in
a time-frequency domain; receiving by a plurality of receivers a
plurality of signals from the plurality of sources, the plurality
of signals including the signal, at least two of the received
signals being correlated; the processor determining a
regularization factor and a weighting factor for each of a
plurality of frequency components of the signal; the processor
determining a likelihood for the direction of arrival for each of
the plurality of frequency components of the signal based on the
regularization factor and the weighting factor; the processor
determining a global likelihood for the direction of arrival of the
signal based on the likelihood for the direction of arrival for
each of the plurality of frequency components of the signal; the
processor updating a plurality of weights of the particle filter;
and the processor generating an estimate of the direction of
arrival of the signal by applying the updated weights of the
particle filter.
2. The method of claim 1, wherein the global likelihood L(.theta.)
is determined by an expression: L ( .theta. ) = n = 1 N L n reg (
.theta. ) .lamda. n , ##EQU00061## with .theta. representing a
vector of all directions-of-arrivals for all the sources in the
plurality of sources, n represents a frequency component, N
represents a total number of frequency components,
L.sub.n.sup.reg(.theta.).lamda..sub.n is a regularized likelihood
of all directions-of-arrivals for all the sources adjusted by a
normalized weighting factor .lamda..sub.n for a frequency component
n.
3. The method of claim 1, wherein the weighting factor for a
frequency component n is a normalized weighting factor
.lamda..sub.n is determined by an expression: .lamda. n = .tau. n n
= 1 N .tau. n , ##EQU00062## wherein .tau. n = .PHI. ( P ^ n (
.theta. ) X n X n , ##EQU00063## with .tau..sub.n representing an
un-normalized weighting factor, {circumflex over
(P)}.sub.n(.theta.) representing a regularized projector for a
frequency component n on a signal subspace, .theta. representing a
vector of all directions-of-arrivals for all the sources in the
plurality of sources, X.sub.n representing a vector of frequency
components n of the plurality of signals received by the plurality
of receivers and .phi. representing a non-negative non-decreasing
weighting function.
4. The method of claim 1, further comprising the processor
determining a Maximum A Posteriori likelihood (MAP) estimate of a
frequency component of the signal.
5. The method of claim 4, wherein the MAP estimate represented by
S.sub.n(.theta.) is determined by an expression:
S.sub.n(.theta.)=(A.sub.n.sup.HA+.lamda..sub.nI).sup.-1A.sub.n.sup.HX.sub-
.n, with .theta. representing a vector of all
directions-of-arrivals for all the sources in the plurality of
sources, A representing a steering matrix, A.sub.n.sup.H
representing a Hermitian of a steering matrix for a frequency
component n, .lamda..sub.n representing a regularization factor for
a frequency component n, I representing an identity matrix and
X.sub.n representing a vector of frequency components n of the
plurality of signals received by the plurality of receivers.
6. The method of claim 1, wherein the updated weights w.sub.i.sup.k
of the particle filter are determined by an expression: w i k = w i
k - 1 p ( X 1 : N | .theta. i k ) i = 1 P w i k - 1 p ( X 1 : N |
.theta. i k ) , ##EQU00064## with k representing an update time, i
representing a particle index, p representing a probability
distribution function, P representing a total number of particles,
X.sub.1:N representing a vector of frequency components at moment
k, and .theta..sub.i.sup.k representing a direction of arrival
estimates of particle i at moment k.
7. The method of claim 6, wherein the global likelihood is adjusted
with a scaling parameter.
8. The method of claim 1, wherein at least one source in the
plurality of sources is a moving source.
9. The method of claim 1, wherein a signal is an acoustic
signal.
10. The method of claim 1, wherein at least one source in the
plurality of sources is associated with a power generator.
11. A system to determine a direction of arrival of a signal
generated by a source in a plurality of sources, comprising: a
memory enabled to store data and instructions; a processor enabled
to retrieve and execute instructions from the memory to perform the
steps: estimating a first direction of arrival for a previous
signal generated by the source with a particle filter operating in
a time-frequency domain; processing a plurality of signals from the
plurality of sources received by a plurality of receivers, the
plurality of signals including the signal, at least two of the
received signals being correlated; determining a regularization
factor and a weighting factor for each of a plurality of frequency
components of the signal; determining a likelihood for the
direction of arrival for each of the plurality of frequency
components of the signal based on the regularization factor and the
weighting factor; determining a global likelihood for the direction
of arrival of the signal based on the likelihood for the direction
of arrival for each of the plurality of frequency components of the
signal; updating a plurality of weights of the particle filter; and
generating an estimate of the direction of arrival of the signal by
applying the updated weights of the particle filter.
12. The system of claim 11, wherein the global likelihood
L(.theta.) is determined by an expression: L ( .theta. ) = n = 1 N
L n reg ( .theta. ) .lamda. n , ##EQU00065## with .theta.
representing a vector of all directions-of-arrivals for all the
sources in the plurality of sources, n represents a frequency
component, N represents a total number of frequency components,
L.sub.n.sup.reg(.theta.).lamda..sub.n is a regularized likelihood
of all directions-of-arrivals for all the sources adjusted by a
normalized weighting factor .lamda..sub.n for a frequency component
n.
13. The system of claim 11, wherein the weighting factor for a
frequency component n is a normalized weighting factor
.lamda..sub.n that is determined by an expression: .lamda. n =
.tau. n n = 1 N .tau. n , ##EQU00066## wherein .tau. n = .PHI. ( P
^ n ( .theta. ) X n X n , ##EQU00067## with .tau..sub.n
representing an un-normalized weighting factor, {circumflex over
(P)}.sub.n(.theta.) representing a regularized projector for a
frequency component n on a signal subspace, .theta. representing a
vector of all directions-of-arrivals for all the sources in the
plurality of sources, X.sub.n representing a vector of frequency
components n of the plurality of signals received by the plurality
of receivers and .phi. representing a non-negative non-decreasing
weighting function.
14. The system of claim 11, further comprising the processor
determining a Maximum A Posteriori likelihood (MAP) estimate of a
frequency component of the signal.
15. The system of claim 14, wherein the MAP estimate represented by
S.sub.n(.theta.) is determined by an expression:
S.sub.n(.theta.)=(A.sub.n.sup.HA+.lamda..sub.nI).sup.-1A.sub.n.sup.HX.sub-
.n, with .theta. representing a vector of all
directions-of-arrivals for all the sources in the plurality of
sources, A representing a steering matrix. A.sub.n.sup.H
representing a Hermitian of a steering matrix for a frequency
component n, .lamda..sub.n representing a regularization factor for
a frequency component n, I representing an identity matrix and
X.sub.n a vector of frequency components n of the plurality of
signals received by the plurality of receivers.
16. The system of claim 11, wherein the updated weights
w.sub.i.sup.k of the particle filter are determined by an
expression: w i k = w i k - 1 p ( X 1 : N | .theta. i k ) i = 1 P w
i k - 1 p ( X 1 : N | .theta. i k ) , ##EQU00068## with k
representing an update time, i representing a particle index, p
representing a probability distribution function, P representing a
total number of particles, X.sub.1:N a vector of frequency
components at moment k, and .theta..sub.i.sup.k representing a
direction of arrival estimates of particle i at moment k.
17. The system of claim 16, wherein the global likelihood is
adjusted with a scaling parameter.
18. The system of claim 11, wherein at least one source in the
plurality of sources is a moving source.
19. The system of claim 11, wherein a signal is an acoustic
signal.
20. The system of claim 11, wherein at least one source in the
plurality of sources is associated with a power generator.
Description
BACKGROUND
[0001] The present invention relates to blind separation of
signals. More specifically, this invention relates to blind
separation of signals also called direction of arrival of signals
generated by sources that are correlated.
[0002] Power generation main components (e.g., turbines,
generators, boilers, transformers as well as related auxiliary
equipment such as boiler feed pumps, cooling water pumps, fans,
valves, exhaust gas cleaning systems) operational issues are
normally characterized by changes in their acoustic signatures (in
sonic and most importantly, ultrasonic ranges), e.g., arcing,
bearing or lubrication issues, loose parts, leaks, etc. Contact
based acoustic emission monitoring of all subcomponents of e.g.,
auxiliary equipment that is non-critical for the continuation of
the plant operation is generally uneconomic. Non-contact,
microphone-based systems could monitor a large area with a
multitude of different types of equipments simultaneously. The
concept of this area monitoring is illustrated in FIG. 1. Here, a
fault is identified to come from the device under test (DUT) number
2 using a microphone array containing sensors or microphones
numbers 1, 2 and 3 and source localization instructions implemented
on a processor.
[0003] Challenges of this scenario are as follows. State-of-the art
source separation and localization approaches that use
time-frequency masking will have problems due to the broad band
characteristics of acoustic emissions from power generation
equipment and the ensuing signal overlap in the time-frequency
domain. In particular, the commonly used sparseness assumption that
signals be disjoint in the time-frequency domain is violated and
algorithms that build on this assumption, like e.g. DUET, will
fail. Another assumption required by many algorithms like ESPRIT,
MUSIC or Weighted Subspace Fitting such as disclosed in "M.
Jansson, A. L. Swindlehurst, and B. Ottersten. Weighted subspace
fitting for general array error models. IEEE Transactions on Signal
Processing, 46: 2484-2498, 1997" and "T. Melia and S. Rickard.
Underdetermined blind source separation in echoic environments
using desprit, 2007" is that all signals be independent. This
assumption may be approximately true for speech signals but not
necessarily for machine vibrations.
[0004] Accordingly, improved and novel methods and systems for
blind separation of signals that are not disjoint and/or
independent and/or uncorrelated are required.
SUMMARY
[0005] In a context of acoustic plant surveillance, a blind
direction of arrival estimation scheme is provided for a nonlinear
1-dimensional array of M receivers of J<M unknown signals
transmitted through an unknown channel and subject to Gaussian
noise. The methods and systems provided in accordance with at least
one aspect of the present invention overcome the conditions of
current state-of-the-art methods that signals be disjoint in the
time-frequency domain (W-orthogonality) and/or statistically
independent. Based on a model that relates signals and
measurements, a likelihood function for the directions of arrivals
is constructed. To select relevant frequencies, a new weighting and
regularization scheme is provided and the joint likelihood function
of all frequencies is interpreted as a cost or objective function.
The time domain is included into the estimation procedure by using
a particle filter that operates on the constructed likelihood
function and a simple transition kernel for the directions. An
approach to evaluate the results is also provided.
[0006] In accordance with an aspect of the present invention a
method is provided for determining a direction of arrival of a
signal generated by a source in a plurality of sources, comprising:
a processor estimating a first direction of arrival for a previous
signal generated by the source with a particle filter operating in
a time-frequency domain, receiving by a plurality of receivers a
plurality of signals from the plurality of sources, the plurality
of signals including the signal, at least two of the received
signals being correlated, the processor determining a
regularization factor and a weighting factor for each of a
plurality of frequency components of the signal, the processor
determining a likelihood for the direction of arrival for each of
the plurality of frequency components of the signal based on the
regularization factor and the weighting factor, the processor
determining a global likelihood for the direction of arrival of the
signal based on the likelihood for the direction of arrival for
each of the plurality of frequency components of the signal, the
processor updating a plurality of weights of the particle filter,
and the processor generating an estimate of the direction of
arrival of the signal by applying the updated weights of the
particle filter.
[0007] In accordance with a further aspect of the present invention
the method is provided, wherein the global likelihood L(.theta.) is
determined by an expression
L ( .theta. ) = n = 1 N L n reg ( .theta. ) .lamda. n ,
##EQU00001##
with .theta. representing a vector of all directions-of-arrivals
for all the sources in the plurality of sources, n represents a
frequency component, N represents a total number of frequency
components, L.sub.n.sup.ref(.theta.).lamda..sub.n is a regularized
likelihood of all directions-of-arrivals for all the sources
adjusted by a normalized weighting factor .lamda..sub.n for a
frequency component n.
[0008] In accordance with yet a further aspect of the present
invention the method is provided, wherein the weighting factor for
a frequency component n is a normalized weighting factor
.lamda..sub.n that is determined by an expression
.lamda. n = .tau. n n = 1 N .tau. n , ##EQU00002##
wherein
.tau. n = .PHI. ( P ^ n ( .theta. ) X n X n , ##EQU00003##
with .tau..sub.n representing an un-normalized weighting factor,
{circumflex over (P)}.sub.n(.theta.) representing a regularized
projector for a frequency component n on a signal subspace, .theta.
representing a vector of all directions-of-arrivals for all the
sources in the plurality of sources, X.sub.n representing a vector
of frequency components n of the plurality of signals received by
the plurality of receivers and .phi. representing a non-negative
non-decreasing weighting function.
[0009] In accordance with a further aspect of the present invention
the method is provided, further comprising the processor
determining a Maximum A Posteriori likelihood (MAP) estimate of a
frequency component of the signal.
[0010] In accordance with a further aspect of the present invention
the method is provided, wherein the MAP estimate represented by
S.sub.n(.theta.) is determined by an expression
S.sub.n(.theta.)=(A.sub.n.sup.HA+.lamda..sub.nI).sup.-1A.sub.n.sup.HX.sub-
.n, with .theta. representing a vector of all
directions-of-arrivals for all the sources in the plurality of
sources, A representing a steering matrix. A.sub.n.sup.H
representing a Hermitian of a steering matrix for a frequency
component n, .lamda..sub.n representing a regularization factor for
a frequency component n, I representing an identity matrix and
X.sub.n representing a vector of frequency components n of the
plurality of signals received by the plurality of receivers.
[0011] In accordance with a further aspect of the present invention
the method is provided, wherein the updated weights w.sub.i.sup.k
of the particle filter are determined by an expression:
w i k = w i k - 1 p ( X 1 : N | .theta. i k ) i = 1 P w i k - 1 p (
X 1 : N | .theta. i k ) , ##EQU00004##
with k representing an update time, i representing a particle
index, p representing a probability distribution function. P
representing a total number of particles, X.sub.1:N representing a
vector of frequency components at moment k, and .theta..sub.i.sup.k
representing a direction of arrival estimate of particle i at
moment k.
[0012] In accordance with a further aspect of the present invention
the method is provided, wherein the global likelihood is adjusted
with a scaling parameter.
[0013] In accordance with a further aspect of the present invention
the method is provided, wherein at least one source in the
plurality of sources is a moving source.
[0014] In accordance with a further aspect of the present invention
the method is provided, wherein a signal is an acoustic signal.
[0015] In accordance with a further aspect of the present invention
the method is provided, wherein at least one source in the
plurality of sources is associated with a power generator.
[0016] In accordance with an aspect of the present invention a
system is provided to determine a direction of arrival of a signal
generated by a source in a plurality of sources, comprising a
memory enabled to store data and instructions, a processor enabled
to retrieve and execute instructions from the memory to perform the
steps: estimating a first direction of arrival for a previous
signal generated by the source with a particle filter operating in
a time-frequency domain, processing a plurality of signals from the
plurality of sources received by a plurality of receivers, the
plurality of signals including the signal, at least two of the
received signals being correlated, determining a regularization
factor and a weighting factor for each of a plurality of frequency
components of the signal, determining a likelihood for the
direction of arrival for each of the plurality of frequency
components of the signal based on the regularization factor and the
weighting factor, determining a global likelihood for the direction
of arrival of the signal based on the likelihood for the direction
of arrival for each of the plurality of frequency components of the
signal, updating a plurality of weights of the particle filter, and
generating estimates of the direction of arrival of the signal by
applying the updated weights of the particle filter.
[0017] In accordance with another aspect of the present invention
the system is provided, wherein the global likelihood L(.theta.) is
determined by an expression
L ( .theta. ) = n = 1 N L n reg ( .theta. ) .lamda. n ,
##EQU00005##
with .theta. representing a vector of all directions-of-arrivals
for all the sources in the plurality of sources, n represents a
frequency component, N represents a total number of frequency
components, L.sub.n.sup.ref(.theta.).lamda..sub.n is a regularized
likelihood of all directions-of-arrivals for all the sources
adjusted by a normalized weighting factor .lamda..sub.n for a
frequency component n.
[0018] In accordance with yet another aspect of the present
invention the system is provided, wherein the weighting factor for
a frequency component n is a normalized weighting factor
.lamda..sub.n that is determined by an expression:
.lamda. n = .tau. n n = 1 N .tau. n , ##EQU00006##
wherein
.tau. n = .PHI. ( P ^ n ( .theta. ) X n X n , ##EQU00007##
with .tau..sub.n representing an un-normalized weighting factor,
{circumflex over (P)}.sub.n(.theta.) representing a regularized
projector for a frequency component n on a signal subspace, .theta.
representing a vector of all directions-of-arrivals for all the
sources in the plurality of sources, X.sub.n representing a vector
of frequency components n of the plurality of signals received by
the plurality of receivers and .phi. representing a non-negative
non-decreasing weighting function.
[0019] In accordance with yet another aspect of the present
invention the system is provided, further comprising the processor
determining a Maximum A Posteriori likelihood (MAP) estimate of a
frequency component of the signal.
[0020] In accordance with yet another aspect of the present
invention the system is provided, wherein the MAP estimate
represented by S.sub.n(.theta.) is determined by an expression
S.sub.n(.theta.)=(A.sub.n.sup.HA+.lamda..sub.nI).sup.-1A.sub.n.sup.HX.sub-
.n, with .theta. representing a vector of all
directions-of-arrivals for all the sources in the plurality of
sources. A representing a steering matrix, A.sub.n.sup.H,
representing a Hermitian of a steering matrix for a frequency
component n, .lamda..sub.n representing a regularization factor for
a frequency component n, I representing an identity matrix and
X.sub.n representing a vector of frequency components n of the
plurality of signals received by the plurality of receivers.
[0021] In accordance with et another aspect of the present
invention the system is provided, wherein the updated weights
w.sub.i.sup.k of the particle filter are determined by an
expression:
w i k = w i k - 1 p ( X 1 : N | .theta. i k ) i = 1 P w i k - 1 p (
X 1 : N | .theta. i k ) , ##EQU00008##
with k representing an update time, i representing a particle
index, p representing a probability distribution function. P
representing a total number of particles, X.sub.1:N representing a
vector of frequency components at moment k, and .theta..sub.i.sup.k
representing a direction of arrival estimate of particle i at
moment k.
[0022] In accordance with yet another aspect of the present
invention the system is provided, wherein the global likelihood is
adjusted with a scaling parameter.
[0023] In accordance with yet another aspect of the present
invention the system is provided, wherein at least one source in
the plurality of sources is a moving source.
[0024] In accordance with yet another aspect of the present
invention the system is provided, wherein a signal is an acoustic
signal.
[0025] In accordance with yet another aspect of the present
invention the system is provided, wherein at least one source in
the plurality of sources is associated with a power generator.
DESCRIPTION OF THE DRAWINGS
[0026] FIGS. 1 and 3 illustrate a system set-up in accordance with
an aspect of the present invention;
[0027] FIG. 3 illustrate devices applied in creating signals in
accordance with an aspect of the present invention;
[0028] FIG. 4 illustrates signal properties in accordance with an
aspect of the present invention;
[0029] FIGS. 4-20 illustrate results in accordance with an aspect
of the present invention;
[0030] FIG. 21 illustrates a system with a processor in accordance
with an aspect of the present invention;
[0031] FIG. 22 illustrates an application in accordance with an
aspect of the present invention; and
[0032] FIGS. 23-25 illustrate a performance of methods and systems
provided in accordance with one or more aspects of the present
invention.
DESCRIPTION OF AN EMBODIMENT
[0033] Herein a direction-of-arrival (DOA) estimation method and
system is provided for wideband sources that are based on the
Bayesian paradigm. The method provided in accordance with an aspect
of the present invention approximates the posterior probability
density function of the DOAs in time-frequency domain with a
particle filter. In contrast to other methods, no time-averaging is
necessary, therefore moving sources can be tracked. The herein
provided method uses a new low complexity weighting and
regularization scheme to fuse information from different
frequencies and to overcome the problem of overfitting when few
sensors are available.
[0034] The direction of arrival (DOA) problem is an important
problem in array signal processing and closely intertwined with the
discipline of blind source separation (BSS). Both fields have given
rise to rich research and abundant number of methods, undoubtedly
for the many applications in areas as different as radar, wireless
communications or speech recognition as described in "H. Krim and
M. Viberg. Two decades of array signal processing research: the
parametric approach. Signal Processing Magazine, IEEE, 13(4):
67-94, July 1996." DOA estimation requires an antenna array and
exploit time differences of arrival between sensors. Narrowband DOA
algorithms use phase shifts to approximate time differences of
arrival between noisy measurements of different microphones. This
establishes an algebraic statistical relationship between received
and emitted signals. Methods exist that attempt to solve this
problem. Most of these are variants of ESPRIT as described in "R.
Roy and T. Kailath. Esprit-estimation of signal parameters via
rotational invariance techniques. Acoustics, Speech and Signal
Processing, IEEE Transactions on, 37(7):984, July 1989" or MUSIC as
described in "R. Schmidt. Multiple emitter location and signal
parameter estimation. Antennas and Propagation, IEEE Transactions
on, 34(3):276, March 1986" that use subspace fitting techniques as
described in "M. Viberg and B. Ottersten. Sensor array processing
based on subspace fitting. Signal Processing, IEEE Transactions on,
39(5):1110-1121, May 1991" and are fast to compute a solution.
Statistically optimal methods such as Maximum Likelihood (ML) as
described in "P. Stoica and K. C. Sharman. Acoustics, speech and
signal processing, Acoustics, Speech and Signal Processing, IEEE
Transactions on, 38(7):1132, July 1990" or Bayesian methods as
described in "J. Lasenby and W. J. Fitzgerald. A bayesian approach
to high-resolution beamforming. Radar and Signal Processing, IEE
Proceedings F, 138(6):539-544, December 1991" were long considered
intractable as described in "James A. Cadzow. Multiple source
localization--the signal subspace approach. IEEE Transactions on
Acoustics, Speech, and Signal Processing, 38(7):1110-1125. July
1990", but have been receiving more attention recently in for
instance "C. Andrieu and A. Doucet. Joint Bayesian model selection
and estimation of noisy sinusoids via reversible jump mcmc. Signal
Processing, IEEE Transactions on, 47(10):2667-2676, October 1999"
and "J. Huang, P. Xu, Y. Lu, and Y. Sun. A novel bayesian
high-resolution direction-of-arrival estimator. OCEANS, 2001.
MTS/IEEE Conference and Exhibition, 3:1697-1702, 2001."
[0035] Many methods for wideband DOA estimation are formulated in
the time-frequency domain. The narrowband assumption is then valid
for each subband or frequency bin. Incoherent signal subspace
methods (ISSM) compute DOA estimates that fulfill the signal and
noise subspace orthogonality condition in all subbands
simultaneously. On the other hand, coherent signal subspace methods
(CSSM) as described in "H. Wang and M. Kaveh. Coherent
signal-subspace processing for the detection and estimation of
angles of arrival of multiple wide-band sources. Acoustics, Speech
and Signal Processing, IEEE Transactions on, 33(4):823, August
1985" compute a universal spatial covariance matrix (SCM) from all
data. Any narrowband signal subspace method can then be used to
analyze the universal SCM and estimate the DOAs.
[0036] However, good initial DOA estimates are necessary to
correctly cohere the subband SCMs into the universal SCM as
described in "D. N. Swingler and J. Krolik. Source location bias in
the coherently focused high-resolution broad-band beamformer.
Acoustics, Speech and Signal Processing, IEEE Transactions on,
37(1):143-145, January 1989." Methods like BI-CSSM as described in
"Ta-Sung Lee. Efficient wideband source localization using
beamforming invariance technique. Signal Processing. IEEE
Transactions on, 42(6):1376-1387, June 1994" or TOPS as described
in "Yeo-Sun Yoon, L. M. Kaplan, and J. H. McClellan. Tops: new doa
estimator for wideband signals. Signal Processing, IEEE
Transactions on, 54(6):1977, June 2006" were developed to alleviate
this problem. Subspace algorithms use orthogonality of signal and
noise subspaces as the criterion of optimality. Yet, a
mathematically more appealing approach is to ground the estimation
process on a decision theoretic framework, enabling cross-layer
optimization. A prerequisite is the computation of the posterior
probability density function (pdf) of the DOAs, which can be
achieved by employing particle filters. Such an approach is taken
in "W. Ng, J. P. Reilly, and T. Kirubarajan. A bayesian approach to
tracking wideband targets using sensor arrays and particle filters.
Statistical Signal Processing, 2003 IEEE Workshop on, pages
510-513, 2003" where a Bayesian maximum a posteriori (MAP)
estimator is formulated in the time domain.
[0037] As an aspect of the present invention, a Bayesian MAP
estimator is provided using the time-frequency representation of
the signals. The advantage of time-frequency analysis is shown by
techniques used in the field of blind source separation (BSS) such
as DUET described in "Scott Rickard, Radu Balan, and Justinian
Rosca. Real-time time-frequency based blind source separation. In
Proc. of International Conference on Independent Component Analysis
and Signal Separation (ICA2001), pages 651-656, 2001" and DESPRIT
described in "Thomas Melia and Scott Rickard. Underdetermined blind
source separation in echoic environments using desprit, 2007."
These methods exploit dissimilar signal fingerprints to separate
signals and, in a second step, estimate their DOAs. For example,
speech signals are generally disjoint in time-frequency domain and
single DOAs can be computed for each frequency bin separately. A
power weighted histogram combines frequency information and
clustering is used to find the true DOAs.
[0038] A method provided as an aspect of the present invention uses
a novel, heuristical weighting scheme to combine information across
frequencies. A particle filter approximates the posterior density
of the DOAs, which is then used to extract a maximum a posteriori
estimate.
[0039] First a method will be explained in terms of blind
localization. An explanation in terms of direction-of-arrival is
also provided.
[0040] FIG. 1 illustrates a floor plan of a room including a device
under test (DDT) number 1 to 3 as top view with microphone
installation, including microphones 1, 2 and 3. The location of a
faulty subsystem in DUT#2 is found by triangulation using its
airborne acoustic emissions. A method provided herein in accordance
with an aspect of the present invention overcomes the condition of
time-frequency orthogonality and does not impose any conditions on
signal correlations.
[0041] In accordance with at least one aspect of the present
invention a signal model is provided that is used in a localization
method. A source localization approach and a particle filter
extension to this approach will be provided below. Subsequently,
results of simulations as well as field experiments are also
provided.
[0042] Problem Formulation
[0043] Time-discrete mixtures x.sub.1(t), . . . , x.sub.M(t) from
an array of M closely-spaced microphones recording J source signals
s.sub.1(t), . . . , s.sub.J(t), so that
x M ( t ) = j = 1 J s j ( t - .tau. mj ) ( 1 ) ##EQU00009##
[0044] Here, .tau..sub.mj is the delay of signal j at microphone m
with regard to microphone 1. The microphones are linearly aligned
but not necessarily equally spaced with known distances d.sub.m
between microphones 1 and m. The distances are assumed to be small
with regard to source distances. In particular, attenuation of
source signals between microphones is considered to be negligible
and the wave fronts are plane. The delays .tau..sub.mj then follow
the relation .tau..sub.mj=.delta..sub.jd.sub.m, where .delta..sub.j
is the delay per meter of signal j with regard to the array axis.
Note that .delta..sub.j holds the direction of arrival information
of signal j and therefore .delta..sub.j is referred to as
direction. A discrete windowed Fourier transform is used and it is
assumed that F{x(t-.tau.)}=F(x(t)}e.sup.-i.omega..sup.n.sup..tau.
to hold for each discrete frequency .omega..sub.n. In this case,
the relationship between signals and recordings is linear in the
time-frequency (t-f) domain and can be described on a per frequency
basis.
[0045] Using uppercase letters for Fourier representations the
relationship reads
X m ( .omega. n , .tau. k ) = j = 1 J S j ( .omega. n , .tau. k ) -
.omega. n .delta. j d m + V m ( .omega. n , .tau. k ) ( 2 )
##EQU00010##
where .tau..sub.k is the window center and independently and
identically distributed (iid.) zero mean Circularly Symmetric
Complex Gaussian (CSCG) noise V.sub.n:=(V.sub.1(.omega..sub.n,
.tau..sub.k), . . . , V.sub.M(.omega..sub.n, .tau..sub.k)).sup.T
was added to model measurement noise, windowing errors, and model
uncertainties. One may compare (2) with Equations (30)-(31) in "S.
Rickard, R. Balan, and J. Rosca. Real-time time-frequency based
blind source separation. In Proc. of International Conference on
Independent Component Analysis and Signal Separation (ICA2001),
pages 651-656, 2001."
[0046] One goal is to estimate the directions
.delta.:=(.delta..sub.1, . . . , .delta..sub.J).sup.T, where the
number of sources J is known, based on all available measurements
up to time .tau..sub.k. No assumptions are made on source
characteristics S.sub.n:=(S.sub.1(.omega..sub.n,.tau..sub.k), . . .
, S.sub.j(.omega..sub.n,.tau..sub.k)).sup.T.epsilon.C.sup.J and
noise variances.
[0047] Method Outline
[0048] First addressed is an estimation problem for each time-point
.tau..sub.k independently and to construct a log-likelihood
function for the directions .delta..sub.j as a weighted sum of
log-likelihoods of all frequency bins .omega..sub.n. A score factor
will be used as weights to select frequency bins with pronounced
directional characteristics, see further below In a second step,
the problem of overfitting at frequencies where not all signals are
active will be reduced by introducing a regularization scheme. A
time dimension will be integrated and a (nonlinear) filtering
problem will be posed, which will be solved using a Sequential
Importance Resampling (SIR) particle filter.
[0049] Log-Likelihood
[0050] Under a model as defined by equation (2), the observations
X.sub.n:=(X.sub.1(.omega..sub.n,.tau..sub.k), . . . ,
X.sub.M(.omega..sub.n,.tau..sub.k)).sup.T are iid. CSCG random
variables if conditioned on S.sub.n and .delta., and the joint
probability density function (pdf) of all M observations factorizes
into the marginal ones. (See "R. G. Gallager. Circularly-symmetric
gaussian random vectors. Class Notes, January 2008" for a review of
CSCG random variables).
[0051] For the log-likelihood one then has
J - log p ( X n | S n , .delta. ) .varies. m = 1 M X M - j = 1 J S
j - .omega. n .delta. j d m 2 = : L n ( .delta. , S n ) ( 3 )
##EQU00011##
where the t-f index (.omega..sub.n,.tau..sub.k) has been dropped
from X.sub.m(.omega..sub.n,.tau..sub.k) and
S.sub.j(.omega..sub.n,.tau..sub.4). The negative log-likelihood
function is used according to the cost function framework described
below. If only a single frequency is considered, an ML-estimate for
.delta. and S.sub.n can be obtained by minimizing the sum (cost
function L.sub.n). A high number J of degrees of freedom can lead
to overfitting if the sample size M is small. The parameters
.delta..sub.j are shared throughout all frequencies and their ML
estimates can be found by solving a global optimization problem,
which combines all frequencies--or equivalently by maximizing the
joint likelihood of all frequencies. A global log-likelihood
function is obtained by summation over the log-likelihoods of each
frequency, because all noise terms are assumed to be iid.:
L ( .delta. , S 1 , , S N ) = n = 1 N L n ( .delta. , S n ) ( 4 )
##EQU00012##
[0052] Uncovering the Signals
[0053] The approach discussed above faces a problem that comes
twofold: not all sources are active on all frequencies. First,
there may be frequencies on which no source is active but a high
amount of energy is received nonetheless. This can be due to
microphone characteristics or, considering the small sample size M,
be caused by outliers. One intent is to discard frequency bins on
which no signals are present. Second, if signal energies differ
significantly or if not all sources are active on a given
frequency, the small sampling size is prone to overfitting and the
algorithm may deliver misleading results. On such frequencies, the
use of all J available degrees of freedom in estimating the source
signals is restrained by introducing a regularization scheme.
[0054] Both modifications of the likelihood, rejection of frequency
bins and regularization, build on a factor c.sup.opt which is
defined as the maximal amount of received energy that can be
explained by a source at direction .delta.. Now defined is an
energy reduction that is achieved by a fixed angle/signal pair
(.delta.,S),
c ( .omega. , .delta. , S ) = 1 - m = 1 M X m - S - .omega..delta.
d m 2 m = 1 M X m 2 ( 5 ) ##EQU00013##
[0055] This factor will be used later for the evaluation of the
estimates. By maximizing over S, one obtains the maximal achievable
reduction for a given direction
c opt ( .omega. , .delta. ) = max S .di-elect cons. C c ( .omega. ,
.delta. , S ) ( 6 ) ##EQU00014##
[0056] This factor is related to the Bayes Factor used in model
choice, weighing here the assumptions of a one-signal model with
direction .delta. and a no-signal model. (See "C. P. Robert. The
Bayesian Choice. Springer, 2nd edition, 2001" Chapter 7 for a
review of model choice and Bayes Factors.) If
c.sup.opt(.omega.,.delta.) is small, only little cost function
reduction can be achieved for a direction .delta.. Hence, this
frequency bin is of limited help in estimating a source from that
direction. On the other hand, if c.sup.opt(.omega.,.delta.) is
close to 1 the measurements can be matched to a signal coming from
direction .delta.. In other words, it is likely that a signal is
present and in the direction .delta. and that this frequency can be
of use for the estimation. However, care must be taken since some
cost function reduction is possible even if no signals are present.
In fact, the average percentage achievable cost function reduction
under the assumption that no signal is present: i.e.,
X.sub.m=V.sub.m (H.sub.0 hypothesis) can be calculated to be
E H 0 [ c opt ( .omega. , .delta. ) ] = E [ max S .di-elect cons. C
( 1 - 1 m = 1 M V m 2 m = 1 M V m - S - .omega..delta. d m 2 ) ] =
1 M ( 8 ) ##EQU00015##
[0057] It is quickly verified that
S = 1 M m = 1 M W m ##EQU00016##
with W.sub.m=V.sub.me.sup.i.omega..delta.d.sup.m, solves this MMSE
problem. Hence, using the independence of W.sub.m:
E [ max S .di-elect cons. C ( 1 - 1 m = 1 M V m 2 m = 1 M V m - S -
.omega..delta. d m 2 ) ] = E [ ( 1 - 1 m = 1 M W m 2 m = 1 M W m -
1 M m ' = 1 M W m ' 2 ) ] = E [ ( 1 - 1 m = 1 M W m 2 m = 1 ( W m *
W m - 1 M W m * m ' = 1 M W m ' - 1 M W m m ' = 1 M W m ' * + 1 M 2
m ' = 1 M m '' = 1 M W m ' * W m '' * ) ) ] = E [ ( 1 - 1 m = 1 M W
m 2 m = 1 ( W m * W m - 1 M W m * W m ) ) ] = E [ 1 - ( 1 - 1 M ) ]
= 1 M ##EQU00017##
[0058] Regularized Log-Likelihood
[0059] The percentage achievable cost function reduction factors
are now used to create a regularized log-likelihood function. One
intent is to address the problem of large signal energy
disparities: high-energy frequency bins have a large influence on
the global likelihood function and will affect direction estimates
of all signals, even if only a single signal is present at that
frequency. To mitigate these undesired cross-influences, a penalty
.lamda.(.omega..sub.n,.delta..sub.j) for the signal powers
|S.sub.j(.omega..sub.n,.tau..sub.k)|.sup.2
.lamda. ( .omega. n , .delta. j ) = .alpha. ( max j ' c opt (
.omega. n , .delta. j ' ) - c opt ( .omega. n , .delta. j ) ) ( 9 )
##EQU00018##
where .alpha.>0 is a design parameter that can be adapted to the
problem. In doing so one favors directions that have a high
percentage achievable cost function reduction towards directions
that do not and this has proven to work well in cases where signals
have dissimilar fingerprints in the frequency domain.
[0060] Adding this penalty to the log-likelihood (3) yields the
regularized (negative) log-likelihood
L n reg ( .delta. , S n ) = L n ( .delta. , S n ) + j = 1 J .lamda.
( .omega. n , .delta. j ) S j ( .omega. n , .tau. k ) 2 ( 10 )
##EQU00019##
[0061] From a statistical viewpoint, the regularization can also be
understood as using a non-flat exponential prior for the signal
powers; i.e.,
p(|S.sub.j(.omega..sub.n,.tau..sub.k)|.sup.2)=.lamda.(.omega..sub.n,.del-
ta..sub.j)e.sup.-.lamda.(.omega..sup.n.sup.,.delta..sup.j.sup.)|S.sup.j.su-
p.(.omega..sup.n.sup.,.tau..sup.k.sup.)|.sup.2l.sub.[0,.infin.)(|S.sub.j(.-
omega..sub.n,.tau..sub.k)|.sup.2) (11)
and the regularized log-likelihood as posterior density for S.sub.n
conditioned on .delta.. Here, l.sub.[0,.infin.) is the indicator
function of [0,.infin.).
[0062] Global Log-Likelihood
[0063] Next, a global log-likelihood or cost function for .delta.
and S.sub.1, . . . , S.sub.N will be provided, which is the main
building block for the estimator. Equation (4) will be modified to
include the c.sup.opt factors from above as weighting factors
w(.omega..sub.n,.tau..sub.k),
w ^ ( .omega. n , .tau. k , .delta. ) = j = 1 J c opt ( .omega. n ,
.delta. j ) ( 12 ) w ( .omega. n , .tau. k , .delta. ) = w ^ (
.omega. n , .tau. k , .delta. ) n = 1 N w ^ ( .omega. n , .tau. k ,
.delta. ) ( 13 ) ##EQU00020##
Combining these weights with the regularized log-likelihoods for
each frequency yields the global log-likelihood
L ( .delta. , S 1 , , S N ) = n = 1 N w ( .omega. n , .tau. k ,
.delta. ) L n reg ( .delta. , S n ) ( 14 ) ##EQU00021##
[0064] Direct minimization of this J+2NJ dimensional cost function
is computationally prohibitive if real-time applications are in
mind. However, for fixed .delta. each addend is a convex function
of S.sub.n and a minimizer can be found explicitly (see below). The
remainder is a superposition of sine functions of .delta. that are
costly to evaluate, non-convex and possibly having multiple
minima.
[0065] Particle Filter Implementation
[0066] In order to include the time dimension into the estimation
procedure, a Markovian stochastic dynamic system is constructed;
i.e., state equations p(.delta..sub.k, S.sub.1,k, . . . ,
S.sub.N,k|S.sub.1,k-1, . . . , S.sub.N,k-1) for the estimation
variables .delta. and S.sub.n where k indexes time. It is first
assumed that .delta..sub.k be independent of S.sub.n,l for all n,
k, l and then immediately drop S.sub.n from the system, because not
any assumptions are made on the time structure of the source
signals. Dropping S.sub.n from the system can be achieved formally
by using p(S.sub.1,k, . . . , S.sub.N,k|S.sub.1,k-1, . . . ,
S.sub.N,k-1).varies.1. This prior does not cause problems although
it is not a real pdf because the likelihood is well-behaved.
[0067] A particle filter is used to obtain a discrete approximation
of p(.delta..sub.k, S.sub.1,k, . . . , S.sub.N,k|X.sub.1,0:k, . . .
, X.sub.N,0:k), which uses all available data from time .tau..sub.0
to .tau..sub.k. The estimation of this pdf can then be used to
obtain estimates of .delta..sub.k and S.sub.1,k, . . . , S.sub.N,k.
A particle filter iteratively updates p(.delta..sub.k, S.sub.1,k, .
. . , S.sub.N,k|X.sub.1,0:k, . . . , X.sub.N,0:k) as new
measurements become available and a Sampling Importance Resampling
(SIR) variant with resampling will be used only if the
approximation becomes poor. An estimate of the effective number of
particles is used as criterion. A tutorial overview on particle
filters can be found in [1] S. Arulampalam, S. Maskell, N. Gordon,
and T. Clapp. A tutorial on particle filters for on-line
non-linear/non-gaussian bayesian tracking. IEEE Transactions on
Signal Processing, 50: 174-188, 2001.
[0068] Transitions Kernels and Importance Sampling
[0069] The transition probabilities for the source locations serve
two purposes: to describe the source dislocation and to search the
state space for possible sources. It is assumed that source
movements are independent and, therefore, use 1-dimensional random
walks without drift for each state location; i.e.,
p ( .delta. k | .delta. k - 1 ) = j = 1 J p ( .delta. j , k |
.delta. j , k - 1 ) = j = 1 J ( .delta. j , k - 1 , .sigma. 2 ) (
15 ) ##EQU00022##
where N(.XI.,.sigma..sup.2) is used to represent the pdf of a
Gaussian random variable with mean .mu. and variance .sigma..sup.2.
A particle filter ideally samples from
p(.delta..sub.k|.delta..sub.k-1,
X.sub.k).varies.p(X.sub.k|.delta..sub.k)p(X.sub.k|.delta..sub.k-1)
to update its pdf estimate. However, this distribution is too
difficult for sampling. Thus, importance sampling is employed with
p(.delta..sub.k|.delta..sub.k-1) as importance density and
p(X.sub.k|.delta..sub.k) used to update the particle weights
w.sub.k.sup.i=w.sub.k-1.sup.iexp(-.gamma.L.sup.ref(.delta..sub.k.sup.i,
S.sub.1,k.sup.i, . . . , S.sub.N,k.sup.i)) (16)
.gamma.>0 is a design parameter, which is connected to the
measurement noise variance. Resampling takes place if the effective
number of particles
N eff = 1 i = 1 N particles ( w i ) 2 ##EQU00023##
drops below a threshold N.sub.th.
[0070] Signal Estimates
[0071] As the signal transition kernels have been defined to be
flat, no information is conveyed between time steps. Consequently,
an ML estimator is used for the signals (or rather a MAP estimator
if the regularization scheme is viewed as a prior). Even though no
signal estimates are needed for localization, explicit estimates
are intentionally used instead of marginalization to avoid costly
J-dimensional complex integration. ML/MAP estimates can be obtained
analytically thanks to the quadratic exponent of the Gaussian. One
has
L n reg ( .delta. , S n ) = m = 1 M X m - j = 1 J S j -
.omega..delta. j d m 2 + j = 1 J .lamda. ( .omega. n , .delta. j )
S j 2 = m = 1 M ( X m 2 + j = 1 J j ' = 1 J S j * S j ' - .omega. (
.delta. j ' - .delta. j ) d m - j = 1 J ( X m * S j -
.omega..delta. j d m + X m S j * .omega..delta. j d m ) ) + j = 1 J
.lamda. ( .omega. n , .delta. j ) S j 2 ( 17 ) ##EQU00024##
Defining
[0072] ( H ) ij = m = 1 M - .omega. ( .delta. j - .delta. i ) d m (
18 ) ( .LAMBDA. ) ij { .lamda. ( .omega. n , .delta. i ) i = j 0 i
.noteq. j ( 19 ) g i = - m = 1 M X m .omega..delta. i d m ( 20 )
##EQU00025##
(17) can be rewritten in matrix-vector notation as
L.sub.n.sup.ref(.delta.,S.sub.n)=X.sup.HX+S.sup.H(H+.LAMBDA.)S+g.sup.HS+-
S.sup.Hg (21)
and a minimizer can be found to be
S-(H+.LAMBDA.).sup.-1g (22)
[0073] L.sub.n.sup.reg(.delta.,S.sub.n) is clearly convex in S,
hence the minimizer can be found by solving
.gradient.L.sub.n.sup.ref(.delta.,S)=0 where
.gradient.L.sub.n.sup.ref(.delta.,S) is the Wirthinger derivative
of L.sub.n.sup.reg(.delta.,S) w.r.t. S.sup.H; i.e.,
.gradient.L.sub.n.sup.reg(.delta.,S)=(H+.LAMBDA.)S+g=0.
[0074] The signal estimates can be used to update the particle
weights using (16).
[0075] Estimate Evaluation
[0076] To improve interpretability of the particle filter result
and for fine-tuning of the filter, the following evaluation scheme
is provided in accordance with an aspect of the present invention.
Each of the source locations
.delta. j = 1 N particles i = 1 N particles .omega. i .delta. j i
##EQU00026##
is assigned a score c(.delta..sub.j) that quantifies the belief
that a signal is tracked at the correct location. The calculation
of c(.delta..sub.j) follows that of c.sup.opt(.omega.,.delta.) in
(6) and instead of finding the optimal signals for each direction
.delta..sub.j the ML-estimates S from (22) are used. To concentrate
on those frequencies that are relevant for the respective
direction, the factors c.sup.opt(.omega.,.delta..sub.j) and the
received energy E(.omega..sub.n,.tau..sub.k) are used as weights
.kappa. in accordance with an aspect of the present invention. The
c(.delta..sub.j) factors have to be calculated for each
particle.
E ( .omega. n , .tau. k ) = m = 1 M X m ( .omega. n , .tau. k ) (
23 ) .kappa. n , j i = E ( .omega. n , .tau. k ) c opt ( .omega. n
, .delta. j i ) n = 1 N E ( .omega. n , .tau. k ) c opt ( .omega. n
, .delta. j i ) ( 24 ) c ( .delta. j i ) = n = 1 N .kappa. n , j i
c ( .omega. n , .delta. j i , S j i ( .omega. n , .tau. n ) ) ( 25
) c ( .delta. j ) = 1 N particles i = 1 N particles .omega. i c (
.delta. j i ) ( 26 ) ##EQU00027##
[0077] Modified Transition Kernels
[0078] To harmonize both purposes of the transition kernels
p(.delta..sub.j,k|.delta..sub.j,k-1)--state space search and source
dislocation--, a transition variance .sigma..sup.2 is used that
depends on the evaluation measure described above. A small variance
.sigma..sup.2 is preferred when a source is being tracked and a
large variance if this is not the case. Hence, the following
expression is used to replace .sigma..sup.2 in (15) and
.sigma..sub.max.sup.2 as design parameter:
.sigma..sup.2(.delta..sub.k-1, S.sub.1,k-1, . . . ,
S.sub.N,k-1)=(1-c(.delta..sub.j,k-1)).sigma..sub.max.sup.2 (27)
[0079] This will formally break the independence assumption was
made on S.sub.n and .delta..
[0080] Method Summary
[0081] One method of localization provided in accordance with an
aspect of the present invention can be described by the following
steps:
Method 1:
TABLE-US-00001 [0082] Initialize particle states
.delta..sub.j,0.sup.i i = l,..., N.sub.particles, j = l,..., J
.cndot. for k = l, ., . do: - for i = l, ... ,N.sub.particles do: *
Generate proposals .delta..sub.j,k.sup.i according to (15). *
Calculate c.sup.opt (.omega..sub.n, .delta..sub.j,k.sup.i)
according to (6). * Calculate regularization factors
.lamda.(.omega..sub.n,.delta..sub.j,k.sup.i) according to (9). *
Calculate signal ML/MAP estimates S.sub.n,k.sup.i according to
(22). * Calculate likelihoods L.sub.n.sup.reg
(.delta..sup.i,S.sub.n.sup.i) according to (3) and (10). *
Calculate global likelihood L.sup.reg
(.delta..sub.k.sup.i,S.sub.l,k.sup.i,...,S.sub.N,k.sup.i) according
to (14). * Calculate evaluation factors c(.delta..sub.j,k.sup.i)
according to (23). * Update particle weights w.sup.i according to
(16). - end for - if N.sub.eff < N.sub.th * Resample particles.
- end if .cndot. end for
[0083] Simulation Results
[0084] The development of a method as provided herein in accordance
with an aspect of the present invention was influenced by its
performance on a real data set. However, before going into the real
data the performance of the method in a controlled environment
using artificially generated data is examined first. FIG. 2a
illustrates in diagram a floor plan of a boiler feed pump (BFP)
house in the field. The locations of the additive test sources and
the microphones are shown as well as the activity of power
equipment e.g., BFP#12 and BFP#13. FIG. 2b shows a recording setup
#1 with a loudspeaker source. Equipment is pointed out by arrows.
The two BFPs in the back of the room are active and generate high
broadband noise.
[0085] Artificial Data
[0086] To probe the method under ideal conditions; i.e., in
accordance with a system model that is provided in accordance with
an aspect of the present invention, we two test scenarios have been
generated. For the first scenario two different mixtures of
deterministic sine signals are used as sources. The signals
partially overlap in the frequency domain. The source locations are
varying over time and intersect in the second half of the
simulation. Three different signal to noise ratios are used. The
second scenario is composed of three sources that emit the same
signal, which is again a deterministic mixture of sines. Here, the
sources are slowly moving but do not cross.
[0087] Discussion of Results
[0088] The method succeeds in separating signals of similar powers
even if noise was added to the microphones. However, in case of a
signal to signal ratio of 2 the estimates of the weak signal are
affected. If the signal to signal ratio is further increased no
valid results can be obtained.
[0089] Real World Measurements
[0090] Test recordings were performed and this data was used to
evaluate the above method on application relevant data from the
field. In the following the recorded data and the performance of
the current localization method will be discussed. The recordings
were taken in a boiler feed pump house with two different
microphone arrays. The first array consists of three nonlinearly
spaced ultrasonic microphones (4 Hz-80 kHz). The spacing between
the centers of the microphones is 7 mm and 8 mm, respectively. The
second array uses 4 linearly spaced microphones (40 Hz-16 kHz) with
7 mm spacing between the centers of neighboring units. Different
sound signals are generated at seven locations close to the power
equipment to simulate out-of-normal emissions. The setup of the
experiment is illustrated in FIG. 2a. Note that BFP#12, #13 and
some of the auxiliary equipment is running during all experiments
and is generating intense background noise. The recording scenario
is further clarified by FIG. 2b. Here, the recording setup of test
location 1 is shown with a loudspeaker source.
[0091] FIG. 3 provides a view of the various recording equipment
used in the experiments (top row). Moreover, FIG. 3 also
illustrates the different sources that are employed in the
recordings. A number of repeatable artificial signals with defined
properties are emitted using the loudspeakers. These test signals
are a chirp, a set of sinus signals at different frequencies, broad
band white noise, high frequency noise and recorded noise from a
damaged engine. The time-frequency representations of the first
four test signals are illustrated in FIG. 4. The engine-defect
signal was not used in the analysis due to the limited frequency
range. A circular saw and a hammer impact on a metal tube were used
to generate high energy signals additional to the loudspeakers. The
reason for adding sources despite the repeatable loudspeaker
signals is that the high background noise intensity resulted in
inaudibility of some of the test signals. The different test
recordings are listed in Table 1. Note that the circular saw was
only used with the PSC microphone array as its emissions were only
barely audible in this loud environment.
[0092] Overall, 95 experiments were performed resulting in 1.65 GB
of acoustic data. One focus is on a subset of these recordings due
to their large number and the limited gain in information by an
exhaustive evaluation. Results are presented below for the sine
signal and different test locations. Results for used signal types
on test location 2 are also provided.
[0093] Source Localization Results
[0094] In the following, results for the experiments with
sinusoidal signals and different test locations are presented. The
sine test signal was selected, because it has the highest signal to
noise ratio (SNR). Experiment results are documented in FIGS. 5-20
and will be summarized below. Each of these figures shows in a
curve the particle filter result. The marker curve indicates the
score function values c(.delta..sub.j) described earlier. The
spectrogram of the received mixture signal is displayed also in the
result figure. Also the estimated spectrograms of the tracked
signals will be shown. Note that the different test locations are
illustrated in FIG. 2a. FIG. 4 shows time-frequency plots of test
signals: (a) chirp test signal, (b) sine signals at different
frequencies, (c) broadband noise, and (d) narrow-band
high-frequency noise.
[0095] Discussion of Results
[0096] The first four sources were located in front of the machines
at angles of approximately -60, -50, -5, and 40 degrees as
indicated by the dashed lines in the result figures. Locations 2
and 3 are closest to the microphone array and achieve the highest
SNR; therefore leading to the most accurate results. The
loudspeakers produce the clearest signal at 10 kHz between seconds
2 and 3. During this time interval, the signal can be located
correctly at locations 1-3. Some of the sine signals at 2, 5, 15,
and 20 kHz are detected at locations 2 and 3. Being closest to the
microphone array, location 3 is detected correctly during the whole
signal duration. In addition to the small distance, one of the
active machines (BFP#12) is placed behind location 3, which
facilitates signal detection. Throughout all recordings, BFP#13 is
detected consistently at approximately 60 degrees. Locations 5-7
are behind the machines and the SNR at the microphone array is too
low for the method provided herein to detect the appropriate
signals.
[0097] A summary of experiment with two microphone arrays is
provided in the following Table 1.
TABLE-US-00002 TABLE 1 Number of Position Signal Types Recordings
Ultrasonic 1, 2, 3, 4, Chirp, sine tones, white 35 Microphone 5, 6,
7 noise, high freq. noise, Array defect engine noise Hammer on
metal 7 -- Normal operation 1 (BFP#11 & BFP#13) PSC 1, 2, 3, 4,
Chirp, sine tones, white 35 Microphone 5, 6, 7 noise, high freq.
noise, Array defect engine noise Circular saw, hammer on 14 metal
-- Normal operation 1 (BFP#12 & BFP#13) Normal operation 1
(BFP#11 & BFP#13) Switch from BFP#12 to 1 BFP#11
[0098] The following provides an overview of results using the
provided method based on different signals.
[0099] Simulation Results with Artificial Data
[0100] FIG. 5: Simulation with two sources using different
combinations of sine signals (bottom left). Particle filter result
curve with dashed lines indicating the true source positions
(dashed line) and curves for the estimation result (top) as well as
estimated signals (bottom right).
[0101] FIG. 6: Simulation with two sources using different
combinations of sine signals (bottom left) and additive noise
with
P noise P signal 1 = P noise P signal 2 = 10 ##EQU00028##
at all microphones. Particle filter result indicating the true
source positions (dashed line) and the estimation result (top) as
well as estimated signals (bottom right).
[0102] FIG. 7: Simulation with two sources using different
combinations of sine signals (bottom left) and additive noise
with
P noise P signal 1 = 10 and P noise P signal 2 = 5 ##EQU00029##
at all microphones. Particle filter result with the true source
positions and the estimation result (top) as well as estimated
signals (bottom right).
[0103] FIG. 8: Simulation with three sources using the same sine
signal (top left). Particle filter result indicating the true
source positions (dashed line) and the estimation result (top
right) as well as estimated signals (bottom).
[0104] Results for Recordings with Sine Signals
[0105] FIG. 9: Test location 1, sine signal. Particle filter
results indicating the true source (dashed line) and the estimation
result (top). Mean spectrogram of measurements (bottom left) and
estimated spectrograms of first and second source on linear scale
(bottom center and right).
[0106] FIG. 10: Test location 2, sine signal. Particle filter
results indicating the true source position (dashed line) and the
estimation result (top). Mean spectrogram of measurements (bottom
left) and estimated spectrograms of first and second source on
linear scale (bottom center and right).
[0107] FIG. 11: Test location 3, sine signal. Particle filter
results indicating the true source position (dashed line) and the
estimation result (top). Mean spectrogram of measurements (bottom
left) and estimated spectrograms of first and second source on
linear scale (bottom center and right).
[0108] FIG. 12: Test location 4, sine signal. Particle filter
results indicating the true source position (dashed line) and
circles for the estimation result (top). Mean spectrogram of
measurements (bottom left) and estimated spectrograms of first and
second source on linear scale (bottom center and right).
[0109] FIG. 13: Test location 5, sine signal. Particle filter
results indicating the true source position (dashed line) and
circles for the estimation result (top). Mean spectrogram of
measurements (bottom left) and estimated spectrograms of first and
second source on linear scale (bottom center and right).
[0110] FIG. 14: Test location 6, sine signal. Particle filter
results indicating the true source position (dashed line) and the
estimation result (top). Mean spectrogram of measurements (bottom
left) and estimated spectrograms of first and second source on
linear scale (bottom center and right).
[0111] FIG. 15: Test location 7, sine signal. Particle filter
results indicating the true source position (dashed line) and the
estimation result (top). Mean spectrogram of measurements (bottom
left) and estimated spectrograms of first and second source on
linear scale (bottom center and right).
[0112] Results for Recordings with Different Signal Types
[0113] FIG. 16: Test location 2, chirp signal. Particle filter
results indicating the true source position (dashed line) and the
estimation result (top). Mean spectrogram of measurements (bottom
left) and estimated spectrograms of first and second source on
linear scale (bottom center and right).
[0114] FIG. 17: Test location 2, sine signal. Particle filter
results indicating the true source position (dashed line) and the
estimation result (top). Mean spectrogram of measurements (bottom
left) and estimated spectrograms of first and second source on
linear scale (bottom center and right).
[0115] FIG. 18: Test location 2, noise signal. Particle filter
results indicating the true source position (dashed line) and the
estimation result (top). Mean spectrogram of measurements (bottom
left) and estimated spectrograms of first and second source on
linear scale (bottom center and right).
[0116] FIG. 19: Test location 2, high-frequency noise signal.
Particle filter results indicating the true source position (dashed
line) and the estimation result (top). Mean spectrogram of
measurements (bottom left) and estimated spectrograms of first and
second source on linear scale (bottom center and right).
[0117] FIG. 20: Test location 2, hammer on metal signal. Particle
filter results indicating the true source position (dashed line)
and the estimation result (top). Mean spectrogram of measurements
(bottom left) and estimated spectrograms of first and second source
on linear scale (bottom center and right).
[0118] Comparison of Results from Different Test Signals
[0119] To better understand and improve a method provided herein in
accordance with an aspect of the present invention, a broad range
of different source signals was employed. Sines at piecewise
constant frequencies were used to examine whether narrowband high
SNR signals can be detected. In addition to that, a chirp signal
was used, which has the same characteristics as the sine above, but
is centered at a slowly time-varying frequency. The case of
balanced signal powers and low SNRs is evaluated using a broadband,
almost white noise signal with contributions up to 24 kHz and a
high-frequency noise signal with contributions between 20 and 25
kHz. The high-frequency noise signal had a low SNR due to
limitations of the loudspeakers. Finally, in order to generate high
SNR signals, a steel hammer was hit repeatedly on a metal tube.
[0120] Discussion of Results
[0121] The chirp and sine signals can be correctly detected for
most of the signal duration due to their high SNR. The broadband
noise signal provokes a visible influence on the location
estimates. The signal from the hammer impacts on the metal tube has
a high SNR during times when the hammer strokes are audible. The
hammer can be correctly located at these time instances.
[0122] In summary, a method for blind separation of correlated
sources provided herein in accordance with an aspect of the present
invention addresses the shortcomings of previous and known methods
for blind separation of sources. The provided method can utilize
nonlinear and/or non-linearly spaced microphone arrays to cancel
disambiguates in the estimated source location. Furthermore, it
enables the localization of multiple signals simultaneously and
overcomes the assumptions of known methods that signals are
disjoint in the time frequency domain or statistically independent.
Finally, the implemented method automatically separates the
predicted source signals. These and their estimated source
locations can be further analyzed e.g., to detect out-of-normal
sounds and reason about their origin. The source localization
results of the herein provided method show improvements on both
synthetic and real world data.
[0123] In a further embodiment of the present invention the
robustness of source localization and separation is increased by
assuming signal coherence over time. In yet a further embodiment of
the present invention, the channel estimate for each
source-receiver pair is used to address problems in echoic
environments.
[0124] The methods as provided herein are in one embodiment of the
present invention implemented on a system or a computer device. A
system illustrated in FIG. 21 and as provided herein is enabled for
receiving, processing and generating data. The system is provided
with data that can be stored on a memory 1801. Data may be obtained
from a sensor such as a microphone or a microphone array or may be
provided from a data source. Data may be provided on an input 1806.
The processor is also provided or programmed with an instruction
set or program executing the methods of the present invention is
stored on a memory 1802 and is provided to the processor 1803,
which executes the instructions of 1802 to process the data from
1801. Data, such as a source localization data or estimated source
signals or any other signal resulting from the processor can be
outputted on an output device 1804, which may be a display to
display data or a loudspeaker to provide an acoustic signal.
However, the data generated by the processor may also be used to
enable further processing of signals, including to generate a
classification of a signal estimated from a source. The processor
can also have a communication channel 1807 to receive external data
from a communication device and to transmit data to an external
device. The system in one embodiment of the present invention has
an input device 1805, which may be a sensor, a microphone, a
keyboard, a mouse or any other device that can generate data to be
provided to processor 1803. The processor can be dedicated
hardware. However, the processor can also be a CPU or any other
computing device that can execute the instructions of 1802.
Accordingly, the system as illustrated in FIG. 21 provides a system
for data processing resulting from a sensor or any other data
source and is enabled to execute the steps of the methods as
provided herein as an aspect of the present invention.
[0125] Thus, a system and methods have been described herein for
blind localization of correlated sources, for instance in a system
to detect out-of-normal operations in a machine that generates
signals. Such a system is illustrated in FIG. 22. The system 2200
includes a receiver 2202, and a processor 2204. The receiver 2202
has multiple microphones and receives audio signals generated by
different audio sources related to a machine. The audio signals are
sent to the processor 2204. The processor 2204 can be in the form
shown in the preceding figure. The processor 2204 processes the
received signals in accordance with the steps described herein.
Depending on the results of the processing, the processor 2204
generates a signal on output 2206 that indicates an out-of-normal
operating mode.
[0126] Formulation in Terms of Direction of Arrival
[0127] Previously, aspects of the invention are described with
respect to blind localization. Various aspects in the invention can
also be expressed in terms of a direction of arrival problem. The
rest of the specification does so.
[0128] Consider a linear array of M sensors with element spacing
d.sub.m-d.sub.m-1. Impinging on this array are J unknown wavefronts
from different directions .theta..sub.j. The propagation speed of
the wavefronts is c. The number J of sources is assumed to be known
and J.ltoreq.M. The microphones are assumed to be in the far-field
of the sources. In DFT domain, the received signal at the m.sup.th
sensor in the n.sup.th subband can be modeled as:
X m ( .omega. n ) = j = 0 J - 1 S j ( .omega. n ) - .omega. n v m
si n ( .theta. j ) + N m ( .omega. n ) ( 28 ) ##EQU00030##
where S.sub.j(.omega..sub.n) is the jth source signal,
N.sub.m(.omega..sub.n) is noise, and
v m = d m c . ##EQU00031##
The noise is assumed to be circularly symmetric complex Gaussian
(CSCG) and iid within each frequency, that is, variances may vary
with .omega..
[0129] The following definitions are used:
X.sub.n=[X.sub.0(.omega..sub.n) . . .
X.sub.M-1(.omega..sub.n)].sup.T (29)
N.sub.n=[N.sub.0(.omega..sub.n) . . .
N.sub.M-1(.omega..sub.n)].sup.T (30)
S.sub.n=[S.sub.0(.omega..sub.n) . . .
S.sub.J-1(.omega..sub.n)].sup.T (31)
Eq. (29) can be rewritten in matrix-vector notation as:
X.sub.n=A.sub.n(.theta.)S.sub.n+N.sub.n (32)
with the M.times.J steering matrix A.sub.n(.theta.):
A.sub.n(.theta.)=[a(.omega..sub.n,.theta..sub.0), . . . ,
a(.omega..sub.n,.theta..sub.J-1)] (33)
whose columns are the M.times.1 array manifolds:
a(.omega..sub.n,.theta..sub.j)=[1e.sup.-i.omega..sup.n.sup.v.sup.1.sup.s-
in(.theta..sup.j.sup.) . . .
e.sup.-i.omega..sup.n.sup.v.sup.M-1.sup.sin(.theta..sup.j.sup.)].sup.T
(34)
and
.theta.=[.theta..sub.0, . . . , .theta..sub.J-1] (35)
[0130] Subspace Methods
[0131] The most commonly used methods to solve the DOA problem
compute signal and noise subspaces from the sample covariance
matrix of the received data and choose those .theta..sub.j whose
corresponding array manifolds a(.theta..sub.j) are closest to the
signal subspace, i.e.,
.theta. ^ = argmax .theta. P n ( .theta. ) X n ( 37 )
##EQU00032##
where the columns of E.sub.N form an orthonormal basis of the noise
subspace. Reference to .omega..sub.n has been omitted
intentionally. Incoherent methods compute signal and noise
subspaces E.sub.N(.omega..sub.n) for each subband and the
.theta..sub.j are chosen to satisfy (36) on average. Coherent
methods, conversely, compute reference signal and noise subspaces
by transforming data from all frequencies to a reference frequency.
The orthogonality condition (36) is then verified for the reference
array manifold only. These methods, of which CSSM and WAVES are two
representatives, show significantly better performance than
incoherent methods, especially for highly correlated and low SNR
signals. But they require good initial DOA estimates for the
focusing procedure to be effective and it is not obvious how these
can be obtained.
[0132] Maximum Likelihood Methods
[0133] In contrast to subspace algorithms, ML methods compute the
signal subspace from the A matrix and choose .theta..sub.0, . . . ,
.theta..sub.J-1 that best fit the observed data in terms of
maximizing its projection on that subspace:
.theta. ^ j = argmin .theta. a ( .theta. ) H E N E N H a ( .theta.
) ( 36 ) ##EQU00033##
where P.sub.n=A(A.sup.HA).sup.-1A.sup.H is a projection matrix on
the signal subspace spanned by the columns of A.sub.n(.theta.).
[0134] If noise variances are equal across frequencies, an overall
log-likelihood function for the wideband problem can be obtained by
summing (37) across frequencies. The problem of varying noise
variances across frequencies has not been addressed to date.
[0135] In "C. E. Chen. F. Lorenzelli, R. E. Hudson, and K. Yao.
Maximum likelihood doa estimation of multiple wideband sources in
the presence of nonuniform sensor noise. EURASIP Journal on
Advances in Signal Processing, 2008: Article ID 835079, 12 pages,
2008. doi:10.1155/2008/835079, 2008" the case of non-uniform noise
across sensors, but constant across frequencies, has been
investigated in great detail. ML methods offer higher flexibility
regarding array layouts and signal correlations than subspace
methods and generally show better performance for small sample
sizes, but the nonlinear multidimensional optimization in (37) is
computationally complex. Recently, importance sampling methods have
been proposed for the narrowband case to solve the optimization
problem efficiently as described in "Huigang Wang, S. Kay, and S.
Saha. An importance sampling maximum likelihood direction of
arrival estimator. Signal Processing, IEEE Transactions on, 56(10):
5082-5092, 2008." The particle filter employed in a method provided
herein in accordance with an aspect of the present invention
approaches the nonlinear optimization in a similar way.
[0136] DOA Estimation
[0137] Under the model (28), the observations X.sub.0, . . . ,
X.sub.M-1 are iid CSCG random variables, if conditioned on S.sub.n
and .theta. the joint pdf, factorize into the marginal ones. For
the negative log-likelihood one has:
- log p ( X n | S n , .theta. ) .varies. .varies. m = 0 M - 1 X m -
j = 0 J - 1 S j - .omega. n v m si n ( .theta. j ) 2 = .DELTA. Ln (
.theta. , S N ) ( 38 ) ##EQU00034##
It is common to compute the solution for S.sub.n:
S.sub.n(.theta.)=A.sub.n(.theta.).sup..dagger.X.sub.n (39)
given fixed .theta. and solve for .theta. by maximizing the
remaining concentrated log-likelihood (A.sup..dagger. denotes the
Moore-Penrose inverse of A). Eq. (39) is a simple least squares
regression and great care must be taken with the problem of
overfitting the data, especially if at some frequency there are
signals with significantly more power than other signals.
[0138] If the noise variances .sigma..sub.n.sup.2 were known, a
global (negative) log-likelihood could be computed by summing the
frequency likelihoods:
L ( .theta. , S 0 : N - 1 ) = n = 0 N - 1 L n ( .theta. , S n )
.sigma. n 2 ( 40 ) ##EQU00035##
This criterion function has been stated previously and was
considered intractable (in 1990) as described in "James A. Cadzow.
Multiple source localization--the signal subspace approach. IEEE
Transactions on Acoustics, Speech, and Signal Processing,
38(7):1110-1125, July 1990." Below a particle filter method is
provided as an aspect of the present invention to solve the
filtering problem for multiple snapshots that naturally solves the
optimization problem as a byproduct. In practical applications, a
good weighting scheme is necessary, as will be provided below as an
aspect of the present invention, as well as a regularization scheme
that improves performance.
[0139] Weighting
[0140] The noise variance .sigma..sup.2 in equation (40) cannot be
estimated from a single snapshot. Instead, the noise variances are
re-interpreted as weighting factors .tau..sub.n, a viewpoint that
is taken by BSS algorithms like DUET. In practice, the signal
bandwidths may not be known exactly and in some frequency bins the
assumption of J signals breaks down. The problem of overfitting
becomes severe in these bins and including them in the estimation
procedure can distort results. The following weights are provided
to account for inaccurate modeling, high-noise bins, and outlier
bins:
.tau. ^ n = .PHI. ( P n ( .theta. ) X n X n ( 41 ) .tau. n = .tau.
^ n n = 0 N - 1 .tau. ^ n ( 42 ) ##EQU00036##
where .phi. is a non-negative non-decreasing weighting function.
Its argument measures the portion of the received signal that can
be explained given the DOA vector .theta.. .tau..sub.n are the
normalized weights.
[0141] Overfitting
[0142] As indicated above, computing the ML estimate for S.sub.n in
(37) can be subject to overfitting, especially if the number of
sensors is small or the J signal assumption breaks down. In
ridge-regression, penalty terms are introduced for the estimation
variables and in Bayesian analysis these translate to prior
distribution for the S.sub.n. Accordingly, the ML estimates of
S.sub.n as an aspect of the present invention are substituted with
MAP approximates under exponential priors for the signal powers. If
one chooses:
p(|S.sub.j(.omega..sub.n)|.sup.2)=.lamda.e.sup.-.lamda.|S.sup.j.sup.(.om-
ega..sup.n.sup.)|.sup.2II.sub.[0,.infin.)(|S.sub.j(.omega..sub.n)|.sup.2)
(43)
with .lamda.>0 and II being the indicator function, the negative
log-posterior for S.sub.n given .theta. can be written:
L n reg ( .theta. , S n ) = L n ( .theta. , S n ) + j = 0 J - 1
.lamda. S j ( .omega. n ) 2 . ( 44 ) ##EQU00037##
[0143] Particle Filter
[0144] The likelihood function for S.sub.n and .theta. can be
computed from equation (40):
p(X.sub.0:N-1|.theta.,S.sub.0:N-1).varies.e.sup.-.gamma.L(.theta.,S.sup.-
0:N-1.sup.) (45)
This is the true likelihood function only if the true noise
variance at frequency n is
.sigma..sub.n.sup.2=(.gamma..tau..sub.n).sup.-1. For what follows
this is assumed to be the case.
[0145] The time dimension is included into the estimation procedure
by means of a particle filter. Starting with the old posterior
pdf:
p.sup.k-1=p(.theta..sup.k-1,S.sub.0:N-1.sup.k-1|I.sup.k-1) (46)
where I.sup.k-1 denotes all information (measurements) up to time
instant k-1, Bayes' rule and elementary calculations are used to
update the posterior as new measurements become available:
p.sup.k.varies.l.sup.kz.sup.k,k-1p.sup.k-1 (47)
with the likelihood:
l.sup.k=p(X.sub.0:N-1.sup.k|.theta..sup.k,S.sub.0:N-1.sup.k)
(48)
and a first order Markov transition kernel for the estimation
variables .theta. and S.sub.n:
z.sup.k,k-1=p(.theta..sup.k,S.sub.0:N-1.sup.k|.theta..sup.k-1,S.sub.0:N--
1.sup.k-1) (49)
[0146] It is assumed that source movements and signals are
independent and that signals are independent between time steps. It
is furthermore assumed that sources move independently. Hence, the
transition kernel can be separated:
z.sup.k,k-1=p(S.sub.0:N-1.sup.k).PI..sub.j=0.sup.J-1p(.theta..sub.j.sup.-
k|.theta..sub.j.sup.k-1) (50)
[0147] The problem of particle depletion is addressed through
resampling if the effective number of particles
N.sub.eff=(.SIGMA..sub.i=0.sup.P-1(w.sub.(i).sup.k).sup.2).sup.-1
falls below a predetermined threshold, as is for instance described
in "Sanjeev Arulampalam, Simon Maskell, Neil Gordon, and Tim Clapp.
A tutorial on particle filters for on-line non-linear/non-gaussian
Bayesian tracking. IEEE Transactions on Signal Processing,
50:174-188, 2001."
[0148] For .theta..sub.j the following modification of a random
walk is used as transition density:
p ( .theta. j k | .theta. j k - 1 ) = .alpha. U [ - .pi. 2 , .pi. 2
] + ( 1 - .alpha. ) P [ - .pi. 2 , .pi. 2 ] N ( .theta. j k - 1 ,
.sigma. .theta. 2 ) ( 51 ) ##EQU00038##
where 0.ltoreq..alpha.<<1,
U [ - .pi. 2 , .pi. 2 ] ##EQU00039##
denotes the pdf of a uniform distribution on
[ - .pi. 2 , .pi. 2 ] , ##EQU00040##
and
P [ - .pi. 2 , .pi. 2 ] N ( .theta. j k - 1 , .sigma. .theta. 2 )
##EQU00041##
denotes the pdf of a normal distribution with mean
.theta..sub.j.sup.k-1 and variance .sigma..sub..theta..sup.2
projected on
[ - .pi. 2 , .pi. 2 ] .theta. ' = a sin ( sin ( .theta. ) ) .
##EQU00042##
A small world proposal density according to "Yongtao Guan, Roland
Fleissner, Paul Joyce, and Stephen M. Krone. Markov chain monte
carlo in small worlds. Statistics and Computing, 16:193-202, June
2006" is used, because this is likely to speed up convergence,
especially in the present case with multimodal likelihood
functions.
[0149] Computer Simulations
[0150] The method herein provided in accordance with an aspect of
the present invention was tested with three different computer
simulated scenarios. In all scenarios, equipowered sources that
emit temporally white Gaussian noise signals are in the far-field
of a linear array of M sensors. Any two different source signals
are correlated with .rho..epsilon.[-1,1]; sensors are subject to
iid white Gaussian noise. An L-point FFT with rectangular window
was performed on the sensor output and further processing was done
on N frequency bins within the sensor passband f.sub.0.+-..DELTA.f.
SNR values are per sensor and for the passband. In the first two
scenarios, intersensor spacing is
d = .lamda. 0 2 ##EQU00043##
between all elements (ULA) where
.lamda. 0 = c f 0 . ##EQU00044##
[0151] For the particle filter, a uniform initial distribution was
used for the .theta..sub.(i),j of all P particles and
.phi.(x)=x.sup.4 as weighting function. The bin size of the
histogram used for extraction of the MAP estimate is .DELTA..sub.h.
For .gamma. the following heuristic was used:
.gamma. k = 10 i = 0 P - 1 L ( .theta. ( i ) k , S ( i ) , 0 : N -
1 k ( 52 ) ##EQU00045##
[0152] WAVES, CSSM, IMUSIC, and TOPS compute DOA estimates based on
the current and the Q preceding blocks. This allows for on-line
dynamic computations and is substantially different from processing
a complete time series offline. WAVES and CSSM use RSS focusing
matrices as described in "H. Hung and M. Kaveh. Focusing matrices
for coherent signal-subspace processing. Acoustics, Speech and
Signal Processing, IEEE Transactions on, 36(8):1272-1281, August
1988" to cohere the sample SCMs with the true angles as focusing
angles. This is an unrealistic assumption but provides an upper
bound on performance for coherent methods. The WAVES algorithm is
implemented as described in "E. D. di Claudio and R. Parisi. Waves:
weighted average of signal subspaces for robust wideband direction
finding. Signal Processing, IEEE Transactions on, 49(10):2179,
October 2001" and Root-MUSIC was used herein as narrowband
technique for both CSSM and WAVES.
[0153] The parameter values of the three scenarios are listed in
the following Table 2.
TABLE-US-00003 TABLE 2 M Source Positions f.sub.x f.sub.0 .DELTA.f
B L N Q P .sigma..sub.0.sup.2 .alpha. .DELTA..sub.h Scenario 1 10
8.degree., 13.degree., 33.degree., 37.degree. 400 Hz 100 Hz 40 Hz
200 256 52 25 2000 (0.5.degree.).sup.2 0.03 0.5.degree. Scenario 2
7 8.degree., 13.degree., 33.degree. 44 kHz 10 kHz 9.9 kHz 220 1024
462 88 200 (1.degree.).sup.2 0.03 .sup. 3.degree. Scenario 3 4
moving 400 Hz 100 Hz 100 Hz 400 256 52 -- 1000 (3.degree.).sup.2
0.05 0.5.degree.
[0154] As an evaluation measure p.sub.all(.DELTA..theta.) was used,
which reflects the percentage of blocks where all sources are
detected within .+-..DELTA..theta. of the true DOA. Each simulation
consists of (B+Q)L samples and evaluation p.sub.all(.DELTA..theta.)
is averaged over the last B blocks. All results are based on 200
Monte Carlo runs for each combination of parameters.
[0155] The first scenario was used in "H. Wang and M. Kaveh.
Coherent signal-subspace processing for the detection and
estimation of angles of arrival of multiple wide-band sources.
Acoustics, Speech and Signal Processing, IEEE Transactions on,
33(4):823, August 1985" and "E. D. di Claudio and R. Parisi. Waves:
weighted average of signal subspaces for robust wideband direction
finding. Signal Processing, IEEE Transactions on, 49(10):2179,
October 2001."
[0156] FIG. 23 illustrates a scenario 1: Percentage of blocks where
all sources are detected within 2.degree. versus SNR for different
values of the source correlation factor .rho.. The particle filter
can discriminate correlated signals at various SNRs in contrast to
CSSM. The results in FIG. 23 show that the particle filter based
method can resolve closely spaced signals at low SNR values and for
arbitrary correlations. In contrast, the resolution threshold for
CSSM increases with correlation (no spatial averaging was applied
as this restricts the possible array geometries). The WAVES results
were nearly identical with the CSSM results and are not shown for
clarity. None of the incoherent methods succeeded in resolving all
four sources in this low SNR scenario. The graphs show that the
proposed algorithm outperforms CSSM in dynamic scenarios averaging
length Q corresponds to approximately 4 seconds--even when perfect
focusing is used.
[0157] Scenario 2 is illustrated in FIG. 24. p.sub.all(2.5.degree.)
versus SNR for .rho.=0 (straight lines) and .rho.=0.1 (dashed
lines). For the second scenario (see FIG. 24), parameters were used
relevant for audio signals. To achieve realtime computations on a
dual-core laptop, the particle filter method uses only the top 10%
high energy frequency bins (N=47). The bin size of the histogram
used for MAP extraction was 3.degree.. The results show that a
suboptimal particle filter with few particles using only 10% of the
available information is a valid alternative to CSSM (again, the
WAVES results were nearly identical and are not shown) and
outperforms (MUSIC.
[0158] A Scenario 3 is illustrated in FIG. 25, showing in a graph
source positions against time. In this third scenario, as
illustrated in FIG. 25, the potential of the herein provided method
to track moving sources is shown. A non-uniform, linear array of
M=4 sensors with coordinates d.sub.0, . . . , d.sub.M-1 and
distances d.sub.m-d.sub.m-1=d+.DELTA.d where
d = .lamda. 0 c ##EQU00046##
and .DELTA.d.about.U.sub.[-0.2d,0.2d]. The signals are concentrated
in the signal passband [f.sub.0-.DELTA.f.sub.SRC,
f.sub.0+.DELTA.f.sub.SRC].OR right.[f.sub.0-.DELTA.f,
f.sub.0+.DELTA.f] with .DELTA.f.sub.SRC=20 Hz and an SNR of 0 dB
total signal power to total noise power was used. The current
method succeeds in estimating the correct source locations of
moving sources before they cross and after they cross.
[0159] In one embodiment of the present invention the method as
provided herein includes an estimation of the number of sources in
the wake of "C. Andrieu and A. Doucet. Joint Bayesian model
selection and estimation of noisy sinusoids via reversible jump
mcmc. Signal Processing, IEEE Transactions on, 47(10):2667-2676,
October 1999" using reversible jump Markov Models. This will
prevent the "ghost" estimates when sources cross.
[0160] Aspects of the present invention are also explained in yet
the following description.
[0161] The following provides a further description of a direction
of arrival estimation scheme for wideband sources that is based on
the Bayesian paradigm in accordance with a aspect of the present
invention. The provided method of which an implementation is
identified under the acronym "MUST" approximates the posterior
probability density function of the DOAs in timefrequency domain
with a particle filter. In contrast to other prior methods, no
time-averaging is necessary, therefore moving sources can be
tracked. The herein provided method uses a new low complexity
weighting and regularization scheme to fuse information from
different frequencies and to overcome the problem of overfitting
when few sensors are available. Decades of research have given rise
to many methods that solve the direction of arrival (DOA)
estimation problem and these methods find application in fields
like radar, wireless communications or speech recognition as
described in "H. Krim and M. Viberg. Two decades of array signal
processing research: the parametric approach. Signal Processing
Magazine, IEEE, 13(4): 67-94, July 1996." DOA estimation requires a
sensor array and exploits time differences of arrival between
sensors. Narrowband methods approximate these differences with
phase shifts. Most of the existing methods for this problem are
variants of ESPRIT as described in "R. Roy and T. Kailath.
Esprit-estimation of signal parameters via rotational invariance
techniques. Acoustics, Speech and Signal Processing, IEEE
Transactions on, 37(7):984, July 1989" or MUSIC as described in "R.
Schmidt. Multiple emitter location and signal parameter estimation.
Antennas and Propagation, IEEE Transactions on, 34(3):276, March
1986" that use subspace fitting techniques as described in "M.
Viberg and B. Ottersten. Sensor array processing based on subspace
fitting. Signal Processing, IEEE Transactions on, 39(5):1110-1121,
May 1991" and are fast to compute a solution. In general, the
performance of subspace based methods degrades with signal
correlation. Statistically optimal methods such as Maximum
Likelihood (ML) as described in "P. Stoica and K. C. Sharman.
Acoustics, speech and signal processing, Acoustics, Speech and
Signal Processing, IEEE Transactions on, 38(7):1132. July 1990" or
Bayesian methods as described in "J. Lasenby and W. J. Fitzgerald.
A bayesian approach to high-resolution beamforming. Radar and
Signal Processing, IEE Proceedings F, 138(6):539-544, December
1991" were long considered intractable as described in "James A.
Cadzow. Multiple source localization--the signal subspace approach.
IEEE Transactions on Acoustics, Speech, and Signal Processing,
38(7):1110-1125, July 1990", but have been receiving more attention
recently such as in "C. Andrieu and A. Doucet. Joint Bayesian model
selection and estimation of noisy sinusoids via reversible jump
mcmc. Signal Processing, IEEE Transactions on, 47(10):2667-2676,
October 1999" and "J. Huang. P. Xu, Y. Lu, and Y. Sun. A novel
bayesian high-resolution direction-of-arrival estimator. OCEANS,
2001. MTS/IEEE Conference and Exhibition, 3:1697-1702, 2001."
[0162] Methods for wideband DOA are mostly formulated in the
time-frequency (t-f) domain. The narrowband assumption is then
valid for each subband or frequency bin. Incoherent signal subspace
methods (ISSM) compute DOA estimates that fulfill the signal and
noise subspace orthogonality condition in all subbands
simultaneously. On the other hand, coherent signal subspace methods
(CSSM) such as described in "H. Wang and M. Kaveh. Coherent
signal-subspace processing for the detection and estimation of
angles of arrival of multiple wide-band sources. Acoustics. Speech
and Signal Processing, IEEE Transactions on, 33(4):823, August
1985" compute a universal spatial covariance matrix (SCM) from all
data. Any narrowband signal subspace method can then be used to
analyze the universal SCM. However, good initial estimates are
necessary to correctly cohere the subband SCMs into the universal
SCM such as described in "D. N. Swingler and J. Krolik. Source
location bias in the coherently focused high-resolution broad-band
beamformer. Acoustics. Speech and Signal Processing, IEEE
Transactions on, 37(1):143-145, January 1989." Methods like BI-CSSM
as described in "Ta-Sung Lee. Efficient wideband source
localization using beamforming invariance technique. Signal
Processing, IEEE Transactions on, 42(6):1376-1387, June 1994" or
TOPS as described in "Yeo-Sun Yoon, L. M. Kaplan, and J. H.
McClellan. Tops: new doa estimator for wideband signals. Signal
Processing, IEEE Transactions on, 54(6):1977, June 2006" were
developed to alleviate this problem.
[0163] Subspace methods use orthogonality of signal and noise
subspaces as the criterion of optimality. Yet, a mathematically
more appealing approach is to ground the estimation on a decision
theoretic framework. A prerequisite is the computation of the
posterior probability density function (pdf) of the DOAs, which can
be achieved with particle filters. Such an approach is taken in "W.
Ng, J. P. Reilly, and T. Kirubarajan. A bayesian approach to
tracking wideband targets using sensor arrays and particle filters.
Statistical Signal Processing, 2003 IEEE Workshop on, pages
510-513, 2003", where a Bayesian maximum a posteriori (MAP)
estimator is formulated in the time domain.
[0164] Herein, a Bayesian MAP estimator is presented using the
time-frequency representation of the signals. The advantage of
time-frequency analysis is shown by techniques used in BSS such as
DUET as described in "Scott Rickard, Radu Balan, and Justinian
Rosca. Real-time time-frequency based blind source separation. In
Proc. of International Conference on Independent Component Analysis
and Signal Separation (ICA2001), pages 651-656, 2001" and DESPRIT
as described in "Thomas Melia and Scott Rickard. Underdetermined
blind source separation in echoic environments using desprit,
2007." These methods exploit dissimilar signal fingerprints to
separate signals and work well for speech signals.
[0165] The presented multiple source tracking of which an
implementation or embodiment is called MUST method uses a novel,
heuristical weighting scheme to combine information across
frequencies. A particle filter approximates the posterior density
of the DOAs and a MAP estimate is extracted. Some widely used other
methods will be discussed and a detailed description of MUST will
be provided. Simulation results of using the herein provided MAP
estimator in accordance with an aspect of the present invention and
in comparison with WAVES as described in "E. D. di Claudio and R.
Parisi. Waves: weighted average of signal subspaces for robust
wideband direction finding. Signal Processing, IEEE Transactions
on, 49(10):2179, October 2001", CSSM, and IMUSIC will be
provided.
[0166] A Further Problem Formulation and Related Work
[0167] Consider a linear array of M sensors with distances between
sensor 1 and m denoted as d.sub.m. Impinging on this array are J
unknown wavefronts from different directions .theta..sub.j. The
propagation speed of the wavefronts is c. The number J of sources
is assumed to be known and J.ltoreq.M. The microphones are assumed
to be in the far-field of the sources. In DFT domain, the received
signal at the m.sup.th sensor in the n.sup.th subband can be
modeled as:
X m ( .omega. n ) = j = 0 J - 1 S j ( .omega. n ) - .omega. n v m s
i n ( .theta. j ) + N m ( .omega. n ) ( 53 ) ##EQU00047##
where S.sub.j(.omega..sub.n) is the jth source signal,
N.sub.m(.omega..sub.n) is noise and
v m = d m c . ##EQU00048##
The noise is assumed to be circularly symmetric complex Gaussian
(CSCG) and iid within each frequency, that is, variances may vary
with .omega.. If one defines X.sub.n=[X.sub.0(.omega..sub.n) . . .
X.sub.M-1(.omega..sub.n)].sup.T and N.sub.n and S.sub.n similarly,
eq. (53) can be rewritten in matrix-vector notation as:
X=A.sub.n(.theta.)S.sub.n+N.sub.n (54)
with the M.times.J steering matrix:
A.sub.n(.theta.)=[a(.omega..sub.n,.theta..sub.1), . . . ,
a(.omega..sub.n,.theta..sub.j)] (55)
whose columns are the M.times.1 array manifolds
a(.omega..sub.n,.theta..sub.j)=[1e.sup.-i.omega..sup.n.sup.v.sup.1.sup.s-
in(.theta..sup.j.sup.) . . .
e.sup.-i.omega..sup.n.sup.v.sup.M-1.sup.sin(.theta..sup.j.sup.)].sup.T
(56)
[0168] Subspace Methods
[0169] The most commonly used methods to solve the DOA problem
compute signal and noise subspaces from the sample covariance
matrix of the received data and choose those .theta..sub.j whose
corresponding array manifolds a(.theta..sub.j) are closest to the
signal subspace, i.e.,
.theta. ^ j = argmin .theta. a ( .theta. ) H E N E N H a ( .theta.
) ( 57 ) ##EQU00049##
where the columns of E.sub.N form an orthonormal basis of the noise
subspace. Reference to .omega..sub.n has been omitted
intentionally. Incoherent methods compute signal and noise
subspaces E.sub.N(.omega..sub.n) for each subband and the
.theta..sub.j are chosen to satisfy (57) on average. Coherent
methods, conversely, compute reference signal and noise subspaces
by transforming data from all frequencies to a reference frequency.
The orthogonality condition (57) is then verified for the reference
array manifold only. These methods, of which CSSM and WAVES are two
representatives, show significantly better performance than
incoherent methods, especially for highly correlated and low SNR
signals. But they require good initial DOA estimates for the
focusing procedure to be effective and it is not obvious how these
can be obtained.
[0170] Maximum Likelihood Methods
[0171] In contrast to subspace methods, ML methods compute the
signal subspace from the A.sub.n matrix and choose {circumflex over
(.theta.)} that best fits the observed data in terms of maximizing
its projection on that subspace:
.theta. ^ = argmax .theta. P n ( .theta. ) X n ( 58 )
##EQU00050##
where P=A(A.sup.HA).sup.-1 A.sup.H is a projector on the signal
subspace spanned by the columns of A.sub.n(.theta.). This
deterministic ML estimates presumes no knowledge of the signals. If
signals statistics were known, a stochastic ML estimates could be
computed.
[0172] If noise variances are equal for all frequencies, an overall
log-likelihood function for the wideband problem can be obtained by
summing (58) across frequencies. The problem of varying noise
variances has not been addressed to date. In "C. E. Chen, F.
Lorenzelli, R. E. Hudson, and K. Yao. Maximum likelihood doa
estimation of multiple wideband sources in the presence of
nonuniform sensor noise. EURASIP Journal on Advances in Signal
Processing, 2008: Article ID 835079, 12 pages, 2008.
doi:10.1155/2008/835079, 2008" the case of non-uniform noise w.r.t.
sensors, but constant across frequencies, has been investigated. ML
methods offer higher flexibility regarding array layouts and signal
correlations than subspace methods and generally show better
performance for small sample sizes, but the nonlinear
multidimensional optimization in (58) is computationally complex.
Recently, importance sampling methods have been proposed for the
narrowband case to solve the optimization problem efficiently as
described in "Huigang Wang, S. Kay, and S. Saha. An importance
sampling maximum likelihood direction of arrival estimator. Signal
Processing, IEEE Transactions on, 56(10):5082-5092, 2008." The
particle filter employed as an aspect of the present invention and
in the method provided as an aspect of the present invention, for
instance in the MUST implementation, tackles the optimization along
these lines.
[0173] Multiple Source Tracking (MUST)
[0174] Under the model (53), the observations X.sub.0, . . . ,
X.sub.M-1 are iid CSCG random variables if conditioned on S.sub.n
and .theta.. Therefore, the joint pdf factorizes into the
marginals. For .omega..sub.n, the negative log-likelihood is given
by
-log
P(X.sub.n|S.sub.n,.theta.).varies..parallel.X.sub.n-A.sub.n(.theta.-
)S.sub.2).parallel..sup.2 (59)
It is common to compute the ML solution for S.sub.n as
S.sub.n(.theta.)=A.sub.n(.theta.).sup. X.sub.n (60)
with A.sup..dagger. denoting the Moore-Penrose inverse of A.sub.n.
A solution to .theta. is found by solving (60) for a fixed .theta.
and iteratively improving on this estimate by minimizing the
remaining concentrated negative log-likelihood
L.sub.n(.theta.).parallel.X.sub.n-A.sub.n(.theta.)S.sub.n).parallel..sup-
.2 (61)
[0175] If the noise variances .sigma..sub.n.sup.2 were known, a
global (negative) log-likelihood could be computed by summing the
frequency likelihoods:
L ( .theta. ) = n = 1 N L n ( .theta. ) .sigma. n 2 ( 62 )
##EQU00051##
[0176] This criterion function has been stated previously and was
considered intractable (in 1990) as described in "James A. Cadzow.
Multiple source localization--the signal subspace approach. IEEE
Transactions on Acoustics, Speech, and Signal Processing,
38(7):1110-1125, July 1990." Below a particle filter method is
provided as an aspect of the present invention to solve the
filtering problem for multiple snapshots that naturally solves the
optimization problem as a byproduct. In practical applications, a
regularization scheme as provided in accordance with an aspect of
the present invention will improve performance. Furthermore,
weighting of the frequency bins is necessary. The low-complexity
approach provided in accordance with an aspect of the present
invention is provided below.
[0177] Regularization
[0178] Eq. (60) is a simple least squares regression and great care
must be taken with the problem of overfitting the data. This
problem is accentuated if the number of microphones is small or if
the assumption of J signals breaks down in some frequency bins.
[0179] In ridge-regression, penalty terms are introduced for the
estimation variables and in Bayesian analysis these translate to
prior distributions for the S.sub.n.
[0180] Similarly to Eq. (60), a MAP estimate of S.sub.n is
S.sub.n(.theta.)=(A.sub.n.sup.HA+.lamda..sub.nI).sup.-1A.sub.n.sup.HX.su-
b.n (63)
where .lamda..sub.n is a regularization parameter for frequency
component .omega..sub.n. We can now eliminate S.sub.n and work
exclusively with the concentrated log-likelihoods that can be
written
L.sub.n.sup.reg(.theta.)=.parallel.(I-P.sub.n(.theta.))X.sub.n.parallel.-
.sup.2 (64)
with
{circumflex over
(P)}.sub.n(.theta.)=A.sub.n(A.sub.n.sup.HA.sub.n+.lamda..sub.nI).sup.-1A.-
sub.n.sup.H (65)
[0181] Weighting
[0182] The noise variance .sigma..sup.2 in equation (62) cannot be
estimated from a single snapshot. Instead, the noise variances are
re-interpreted as weighting factors .lamda..sub.n=.sigma..sup.-2, a
viewpoint that is taken by BSS algorithms like DUET. Equation (62)
can thus be rewritten as
L ( .theta. ) = n = 1 N L n reg ( .theta. ) .lamda. n ,
##EQU00052##
reflecting both the regularization and the weighting. In practice,
the signal bandwidths may not be known exactly and in some
frequency bins the assumption of J signals breaks down. The problem
of overfitting becomes severe in these bins and including them in
the estimation procedure can distort results. The following weights
are provided to account for inaccurate modeling, high-noise bins,
and outlier bins:
.tau. n = .PHI. ( P ^ n ( .theta. ) X n X n ( 66 ) .lamda. n =
.tau. n n = 1 N .tau. n ( 67 ) ##EQU00053##
where .phi. is a non-negative non-decreasing weighting function.
Its argument measures the portion of the received signal that can
be explained given the DOA vector .theta.. .lamda..sub.n are the
normalized weights.
[0183] Particle Filter
[0184] Based on the weighting and regularization schemes, the
concentrated likelihood function reads:
p(X.sub.1:N|.theta.).varies.e.sup.-.gamma.L(.theta.) (68)
where a scaling parameter .gamma. is introduced that determines the
sharpness of the peaks of the likelihood function. A heuristic for
.gamma. is provided below. However, this is the true likelihood
function only if the true noise variance at frequency n is
.sigma..sub.n.sup.2=(.gamma..lamda..sub.n).sup.-1. In what follows
this is assumed to be the case. Now, the time dimension into the
estimation procedure is included.
[0185] First, a Markov transition kernel is defined for the DOAs to
relate information between snapshots:
p ( .theta. j k | .theta. j k - 1 ) = .alpha. U [ - .pi. 2 , .pi. 2
] + ( 1 - .alpha. ) N ( .theta. j k - 1 , .sigma. .theta. 2 ) ( 69
) ##EQU00054##
where 0.ltoreq..alpha.<<1,
U [ - .pi. 2 , .pi. 2 ] ##EQU00055##
denotes the pdf of a uniform distribution on
[ - .pi. 2 , .pi. 2 ] , ##EQU00056##
and N(.theta..sub.j.sup.k-1, .sigma..sub..theta..sub.2) denotes the
pdf of a normal distribution with mean .theta..sub.j.sup.k-1 and
variance .sigma..sub..theta..sup.2. A small world proposal density
according to "Yongtao Guan. Roland Fleissner, Paul Joyce, and
Stephen M. Krone. Markov chain monte carlo in small worlds.
Statistics and Computing, 16:193-202, June 2006" is used, because
this is likely to speed up convergence, especially in the present
case with multimodal likelihood functions.
[0186] Let I.sup.k denote all measurements (information) until
snapshot k. Assume that a discrete approximation of the old
posterior pdf is available:
p ( .theta. k - 1 | I k - 1 ) = i = 1 P w i k - 1 .delta. .theta. i
k - 1 ( 70 ) ##EQU00057##
where the
.delta. .theta. i k - 1 ##EQU00058##
are Dirac masses at .theta..sub.i.sup.k-1. The
.theta..sub.i.sup.k-1 are called particles with associated weights
w.sub.i.sup.k-1. According to Bayes' rule:
p(.theta..sup.k|I.sup.k).varies.p(X.sub.1:N.sup.k|.theta..sup.k)p(.theta-
..sup.k|.theta..sup.k-1) (71)
[0187] An approximation of the new posterior can be obtained as
described in "Sanjeev Arulampalam, Simon Maskell, Neil Gordon, and
Tim Clapp. A tutorial on particle filters for on-line
non-linear/non-gaussian Bayesian tracking. IEEE Transactions on
Signal Processing, 50:174-188, 2001" by resampling each particle
from the transition kernel:
.theta..sub.i.sup.k.about.p(.theta..sub.i.sup.k|.theta..sub.i.sup.k-1),
(72)
updating the weights with the likelihood, and renormalizing
w i k = w i k - 1 p ( X 1 : N | .theta. i k ) i = 1 P w i k - 1 p (
X 1 : N | .theta. i k ) ( 73 ) ##EQU00059##
[0188] The parameter .gamma. influences the reactivity of the
particle filter. A small value puts small confidence into new
measurements while a big value rapidly leads to particle depletion,
i.e., all weight is accumulated by few particles. Through
experimentation it was found that a good heuristic for .gamma. that
reduces the necessity for resampling of the particles while
maintaining the method's speed of adaptation is:
.gamma. = 10 i = 1 P L ( .theta. i ) ( 74 ) ##EQU00060##
[0189] A MAP estimate of .theta. can be obtained from the particles
through use of histogram based methods. However, the particles are
not spared from the permutation invariance problem as described in
"H. Sawada, R. Mukai, S. Araki, and S. Makino. A robust and precise
method for solving the permutation problem of frequency-domain
blind source separation. Speech and Audio Processing, IEEE
Transactions on, 12(5):530-538, 2004." The likelihood function does
not change its value if for some particle .theta..sub.i,j and
.theta..sub.i,j' are interchanged. A simple EM-like clustering
technique is used to account for this problem.
[0190] Complexity
[0191] The main load of the MUST implementation of a method
provided herein in accordance with an aspect of the present
invention is the computation of
(A.sub.n.sup.HA+.lamda..sub.nI).sup.-1 in Eq. (63), which has to be
done for P particles and N frequency bins. Inversion of a real,
positive definite J.times.J matrix requires O(J.sup.2) operations
and can be carried out efficiently using BLAS routines.
Accordingly, the complexity of updating the MAP estimates of
.theta. is O(NPJ.sup.2). Note that the number J of sources also
determines the number P of particles necessary for a good
approximation.
[0192] Results of computer simulations are illustrated in FIGS.
23-25.
[0193] A practical technique to solve the wideband DOA estimation
problem with MAP estimates has been provided above. In contrast to
widely used subband techniques. MAP methods are robust to
correlations and work with few samples, which is important to track
moving sources. A reduced complexity weighting and regularization
scheme has been provided herein to address for model inaccuracies
and small number of microphones, two problems often encountered in
practice. Recent advances in computing power make it feasible to
use this method in the real world as shown by the simulations of
scenario 2 as illustrated in FIG. 24. Further improvement provided
herein over prior methods are: a tuning of the weighting and of the
regularization heuristics and by including an estimation of the
number of sources. Finally, the herein provided method in
accordance with an aspect of the present invention provides the
extension to J>M signals if no more than M signals are active at
a given frequency.
[0194] A summary of the process in accordance with one aspect of
the present invention is provided below. (step 1) Initialization of
equation 69, with .alpha.=1. (step 2) Computation of
S.sub.n(.theta.), is provided in accordance with equation 60. (step
3) Computation of L.sub.n(.theta.) is provided in accordance with
equation 61 (and 60). (step 4) Computation of weights which may be
regularization weights .tau..sub.n to account for noise
.sigma..sub.n is performed in accordance with equations 66 and 67.
(step 5) Computation of likelihood L(.theta.) is performed in
accordance with equation 62. (step 6) Computation of parameter
.gamma. gamma is performed in accordance with equation 74. (step 7)
Computation of a concentrated likelihood is then performed in
accordance with equation 68. (step 8) Updating a prior distribution
is done in two steps. One calculates equation 71 or
p(.theta..sup.k|I.sup.k).varies.p(X.sub.1:N.sup.k|.theta..sub.k)p(.theta.-
.sup.k|.theta..sub.k-1)p(.theta..sup.k-1|I.sup.k-1), given equation
70 and transition kernel equation 69 with .alpha.<1. Weight
updates are provided using equation 73.
[0195] The vector .theta. is updated for each particle i.
[0196] The following references provide background information
generally related to the present invention and are hereby
incorporated by reference: [1] S. Arulampalam, S. Maskell, N.
Gordon, and T. Clapp. A tutorial on particle filters for on-line
non-linear/non-gaussian bayesian tracking. IEEE Transactions on
Signal Processing, 50: 174-188, 2001. [2] R. G. Gallager.
Circularly-symmetric gaussian random vectors. Class Notes, January
2008. [3] M. Jansson, A. L. Swindlehurst, and B. Ottersten.
Weighted subspace fitting for general array error models. IEEE
Transactions on Signal Processing, 46: 2484-2498, 1997. [4] T.
Melia and S. Rickard. Underdetermined blind source separation in
echoic environments using desprit, 2007. [5] S. Rickard, R. Balan,
and J. Rosca. Real-time time-frequency based blind source
separation. In Proc. of International Conference on Independent
Component Analysis and Signal Separation (ICA2001), pages 651-656,
2001. [6] C. P. Robert. The Bayesian Choice. Springer, 2nd edition,
2001. [7] H. Krim and M. Viberg. Two decades of array signal
processing research: the parametric approach. Signal Processing
Magazine, IEEE, 13(4): 67-94, July 1996. [8] R. Roy and T. Kailath.
Esprit-estimation of signal parameters via rotational invariance
techniques. Acoustics, Speech and Signal Processing, IEEE
Transactions on, 37(7):984, July 1989. [9] R. Schmidt. Multiple
emitter location and signal parameter estimation. Antennas and
Propagation, IEEE Transactions on, 34(3):276, March 1986. [10] M.
Viberg and B. Ottersten. Sensor array processing based on subspace
fitting. Signal Processing, IEEE Transactions on, 39(5):1110-1121,
May 1991. [11] P. Stoica and K. C. Sharman. Acoustics, speech and
signal processing, Acoustics, Speech and Signal Processing, IEEE
Transactions on, 38(7):1132, July 1990. [12]J. Lasenby and W. J.
Fitzgerald. A bayesian approach to high-resolution beamforming.
Radar and Signal Processing, IEE Proceedings F, 138(6):539-544.
December 1991. [13] James A. Cadzow. Multiple source
localization--the signal subspace approach. IEEE Transactions on
Acoustics, Speech, and Signal Processing, 38(7):1110-1125, July
1990. [14] C. Andrieu and A. Doucet. Joint Bayesian model selection
and estimation of noisy sinusoids via reversible jump mcmc. Signal
Processing, IEEE Transactions on, 47(10):2667-2676, October 1999.
[15] J. Huang, P. Xu, Y. Lu, and Y. Sun. A novel bayesian
high-resolution direction-of-arrival estimator. OCEANS, 2001.
MTS/IEEE Conference and Exhibition, 3:1697-702, 2001. [16] H. Wang
and M. Kaveh. Coherent signal-subspace processing for the detection
and estimation of angles of arrival of multiple wide-band sources.
Acoustics. Speech and Signal Processing, IEEE Transactions on,
33(4):823, August 1985. [17] D. N. Swingler and J. Krolik. Source
location bias in the coherently focused high-resolution broad-band
beamformer. Acoustics. Speech and Signal Processing, IEEE
Transactions on, 37(1):143-145. January 1989. [18] Ta-Sung Lee.
Efficient wideband source localization using beamforming invariance
technique. Signal Processing, IEEE Transactions on,
42(6):1376-1387, June 1994. [19] Yeo-Sun Yoon, L. M. Kaplan, and J.
H. McClellan. Tops: new doa estimator for wideband signals. Signal
Processing, IEEE Transactions on, 54(6):1977. June 2006. [20] W.
Ng, J. P. Reilly, and T. Kirubarajan. A bayesian approach to
tracking wideband targets using sensor arrays and particle filters.
Statistical Signal Processing, 2003 IEEE Workshop on, pages
510-513, 2003. [21] E. D. di Claudio and R. Parisi. Waves: weighted
average of signal subspaces for robust wideband direction finding.
Signal Processing, IEEE Transactions on, 49(10):2179, October 2001.
[22] C. E. Chen, F. Lorenzelli. R. E. Hudson, and K. Yao. Maximum
likelihood doa estimation of multiple wideband sources in the
presence of nonuniform sensor noise. EURASIP Journal on Advances in
Signal Processing, 2008: Article ID 835079, 12 pages, 2008.
doi:10.1155/2008/835079, 2008. [23] Huigang Wang, S. Kay, and S.
Saha. An importance sampling maximum likelihood direction of
arrival estimator. Signal Processing, IEEE Transactions on,
56(145082-5092, 2008. [24] Yongtao Guan, Roland Fleissner, Paul
Joyce, and Stephen M. Krone. Markov chain monte carlo in small
worlds. Statistics and Computing, 16:193-202, June 2006. [25] H.
Hung and M. Kaveh. Focusing matrices for coherent signal-subspace
processing. Acoustics, Speech and Signal Processing, IEEE
Transactions on, 36(8):1272-1281, August 1988. [26] H. Sawada, R.
Mukai, S. Araki, and S. Makino. A robust and precise method for
solving the permutation problem of frequency-domain blind source
separation. Speech and Audio Processing, IEEE Transactions on,
12(5):530-538, 2004.
[0197] A five page appendix providing further explanation of
various aspect of the present invention is included.
[0198] While there have been shown, described and pointed out
fundamental novel features of the invention as applied to preferred
embodiments thereof, it will be understood that various omissions
and substitutions and changes in the form and details of the
methods and systems illustrated and in its operation may be made by
those skilled in the art without departing from the spirit of the
invention. It is the intention, therefore, to be limited only as
indicated by the scope of the claims.
* * * * *