Methods and systems for tracking of amplitudes, phases and frequencies of a multi-component sinusoidal signal

Abstract

System and methods are provided based on optimization of the weighted log-likelihood. These systems are able to efficiently track dominant sinusoidal components of a real signal in Gaussian noise, provided that the number of the components is known. The algorithm is implemented using simple parallel building blocks involving narrow-band filters that are adaptively self-tuned around the frequencies of the signal components. The algorithm has low computational complexity and provides high estimation accuracy and is also able to track chirp signals. The algorithm is flexible enough to be adjusted to operate in different environments such as for speech signals, by selecting a proper window function. Simulation results confirm that the proposed algorithm is reliable in tracking the frequencies as well as in estimation of the amplitudes of the components. In a chirp environment, the algorithm is able to recognize some frequency cross-overs as long as the amplitudes are different enough around the cross-over moment. Simulations show that the LIR of the algorithm is not affected by the SNR and is inversely proportional to the window length. The effects of the length and type of the window on the frequency resolution and the LIR of the algorithm are discussed. The algorithm is efficiently capable of decomposing speech voiced signals.

Description

FIELD OF THE INVENTION

[0002] The present invention relates generally to the decomposition of an observed signal into its constituent components and, in particular, to the estimation and tracking of characteristic parameters of the constituent components of the observed signal.

BACKGROUND OF THE INVENTION

[0003] The estimation of the characteristic parameters of a dominant sinusoidal component of a received signal is a common feature within radio receivers and many other systems. The characteristic parameters of the dominant sinusoidal component are typically required to facilitate the demodulation and extraction of information carried by the received signal.

[0004] In radio communications Amplitude Modulation (AM) and Frequency Modulation (FM) are commonly used, either in conjunction with one another or independently. If either FM or AM is used exclusively, it is frequency or complex amplitude information, respectively, that is to be extracted from the received signal that has been corrupted after traversing a transmission medium (i.e. a channel). Thus, a single independent demodulation method for either FM or AM would suffice.

[0005] However, a more challenging problem is to devise a method to estimate both the complex amplitude (amplitude and phase) and the frequency of multiple dominant sinusoidal components of a signal, or to devise a radio receiver that can, for a number of prominent sinusoids contained in a received signal, estimate the complex amplitude and frequency of each of the prominent sinusoids contained in the received signal. Such methods have applications in the field of speech-processing, for example speech recognition and speech decomposition. By decomposing a speech signal into a set of amplitudes, frequencies and phases of the dominant sinusoidal components, speech information can be efficiently stored or transmitted and then reconstructed after being read from a memory, or after being transmitted and received.

[0006] Given the already complex nature of this problem it is usually assumed that the environment (i.e. channel) and the parameters to be estimated are reasonably stationary. In other words, it is assumed that the statistical characteristics of the channel and the received signal do not change significantly over a short duration of time.

[0007] Different approaches have been used to address the problem of estimating both the complex amplitude and frequency of each of a number of prominent frequency components contained in a received signal. The many approaches employed include the use of the Discrete Fourier Transform (DFT), Adaptive Line Enhancement (ALE), Extended Kalman Filter (EKF), Maximum Likelihood Estimation (MLE), and Joint Time Frequency Analysis (JTFA).

[0008] The DFT was one of the first methods employed because it readily enabled the separation of the prominent sinusoidal components in the frequency domain. This method is particularly suitable in applications where the parameters are constant, as the DFT can be used in concert with the MLE approach. The combination of the DFT and MLE approaches results in a method that is the equivalent of maximizing the periodgram spectrum. Using DFT and MLE approaches together, the periodgram can be calculated and maximized at discrete frequency points. The prominent components of the received signal can be operated upon independently, avoiding cross-interference between them.

[0009] The combined DFT and MLE method outlined above is not very attractive for real time applications because of the high computational cost and complexity associated with the method. Also, because the method has no memory the tracking is not efficient. Several methods have been suggested to enhance the performance of this approach, however they have been employed with diminishing returns and do not completely resolve the problems. For instance, the hidden Markov model was proposed to improve the efficiency of tracking the parameters and the computational cost of this approach can be reduced by use of the Fast Fourier Transform (FFT), in place of the DFT.

[0010] The MLE method on its own is a powerful approach to parameter estimation and is widely used in signal processing. An iterative method for the MLE extraction of parameters of a harmonic series using the Expectation Maximization (EM) algorithm, where the number of harmonics is assumed to be known, has been suggested for AWGN (Additive White Gaussian Noise) channels with known noise variance. Such methods suffer from poor performance because the effect of the cross-interference of harmonics is ignored. These cross-interferences become significant contributors to the degradation of a signal within a low SNR (signal-to-noise ratio) environment or for a short data length.

[0011] If the characteristic parameters of the signal are slightly non-stationary or if the additive noise has a time varying variance, employing an approach that relies on finite length window, such as the DFT, will inherently lead to the loss of some information in regard to the dynamics of the signal. In previous works using the harmonic model for a speech signal, the amplitudes, phases and frequencies that make-up the speech signal have been estimated. Furthermore, it has also been shown that information about the nature of the slowly time varying parameters can be obtained from the distortions caused by windowing.

[0012] In yet another approach called Adaptive Line Enhancement (ALE), the frequency of the dominant component is estimated by minimizing the output energy of a notch filter. An Adaptive Comb Filter (ACF) that is an extension of ALE has been used to estimate the frequency of harmonic signals.

[0013] Kalman Filtering has also been employed in prior work for different scenarios. Specifically, EKF has been used to track the frequencies, the amplitudes and the phases of harmonic components within a periodic signal corrupted by AWGN. Using similar principles, several non-linear filters have been proposed for the decomposition of signals that are modeled as a sum of jointly modulated amplitude and frequency cosines in an additive noise environment, where the centre frequencies are very slowly time varying. Furthermore, assuming that the noise statistics and the number of superimposed signals are known, an EKF can be designed to track frequency formats of speech.

[0014] In another existing approach the signal to be analyzed is assumed to be a Polynomial Phase Signal (PPS) and unknown parameters are estimated. Several techniques could be employed to resolve this problem, such as FFT. High-resolution frequency estimation methods such as Kumaresan-Tufts, MUSIC and Matrix Pencil are alternatives to estimate the polynomial phase coefficients.

[0015] Exponentially-damped Polynomial Phase Signals (EPPS) have been treated as a special case of PPS's. In such a case, it is typically assumed that a PPS can be modeled as having constant amplitude, thus allowing the use of JTFA tools such as Wigner-Ville distribution and an associated Ambiguity Function. Using the Wigner-Ville distribution an estimation method that selects an optimal time domain window length to resolve the trade-off between the estimation bias and the variance of the unknown frequency can be used.

[0016] In another existing approach researchers have considered non-stationary signals as time-dependent ARMA (auto-regressive moving average) processes, and suggested a general estimation procedure for the ARMA parameters using a set of basis functions. In one such work, estimates for the signal parameters for any non-stationary signal can be obtained using the Non-linear Instantaneous Least Squares (NILS) method. NILS has relation with ML as well as with signal-subspace fitting and linear prediction based estimation approaches.

SUMMARY OF THE INVENTION

[0017] According to one broad aspect, the invention provides a method of tracking amplitude, phase and frequency of a plurality of sinusoidal components in a signal, the method comprising: a) processing the signal to produce a new set of amplitude and phase estimates using a weighted likelihood method; and b) processing the signal to produce a new set of frequency estimates using a weighted likelihood method.

[0018] In some embodiments, the method further comprises sampling the signal to produce a sequence of real-valued samples, wherein steps a) and b) are performed in the digital domain.

[0019] In some embodiments, the method further comprises sampling the signal to produce a sequence of complex-valued samples, wherein steps a) and b) are performed in the digital domain.

[0020] In some embodiments, steps a) and b) are performed in the continuous time domain.

[0021] According to one broad aspect, the invention provides a method of tracking amplitude, phase and frequency of a plurality of sinusoidal components in a signal, the method comprising: for a current update period: i) processing the signal to produce a new set of complex amplitude estimates by: a) for a first input set of estimated complex sinusoidal components, separating components to produce component estimates; ii) processing the signal to produce a new set of estimated complex sinusoidal components by: b) for each component of a second input set of estimated complex sinusoidal components, estimating a frequency deviation estimate; c) adapting a previous set of frequency estimates taking into account an input set of component estimates and the frequency deviation estimates to produce a new set of frequency estimates; and d) converting the new set of frequency estimates to a new set of estimated complex sinusoidal components.

[0022] In some embodiments, the signal is a sequence of samples and processing is done in the digital domain.

[0023] In some embodiments, the processing is done in the continuous time domain.

[0024] In some embodiments, the further comprises: performing cross-interference cancellation on the component estimates to produce a new set of cross-interference cancelled component estimates, and using the new set of cross-interference cancelled estimates as the input set of component estimates in an execution of step c).

[0025] In some embodiments, the method further comprises: performing complex envelope extraction on the component estimates to produce a new set of complex amplitude estimates.

[0026] In some embodiments, the method further comprises: performing complex envelope extraction on the cross-interference cancelled component estimates to produce a new set of complex amplitude estimates.

[0027] In some embodiments, for the first input set of estimated complex sinusoidal components, separating components to produce component estimates is done using a weighted log-likelihood function with a first weighting sequence; for each of the second input set of estimated complex sinusoidal components, estimating the frequency deviation estimate is done using a weighted log-likelihood function with a second weighting sequence.

[0028] In some embodiments, the first and second weighting sequences are the same.

[0029] In some embodiments, step i) comprises a first half-iteration, and step ii) comprises a second half iteration, one first half-iteration and one second half-iteration comprising a complete iteration and wherein for each update period, a plurality of complete iterations are performed to produce the new set of complex amplitude estimates and the new set of estimated complex sinusoidal components.

[0030] In some embodiments, the first input set of estimated complex sinusoidal components and the second set of estimated complex sinusoidal components are initially set to initial values, and thereafter are set to estimated complex sinusoidal components produced by a previous iteration of the method.

[0031] In some embodiments, for each update of the complex amplitude and frequency: the step of processing samples of the sequence of samples to produce a new set of complex amplitude estimates is performed before the step of processing the sequence of samples to produce a new set of estimated complex sinusoidal components; the first input set and the second input set of estimated complex sinusoidal components comprise the new set of estimated complex sinusoidal components determined during a previous update period; wherein the input set of cross-interference cancelled estimates comprises the new set of cross-interference cancelled estimates determined during the current update period.

[0032] In some embodiments, for each update of the amplitude, phase and frequency: the step of processing the signal to produce a new set of estimated complex sinusoidal components is performed before the step of processing the sequence of samples to produce a new set of complex amplitude estimates; the input set of component estimates comprises the set of cross-interference cancelled estimates determined during a previous update period; the first input set of estimated complex sinusoidal components comprises the new set of estimated complex sinusoidal components determined during the current update period and the second input set of estimated complex sinusoidal components comprises the new set of estimated complex sinusoidal components determined during a previous update period.

[0033] In some embodiments, for the first input set of estimated complex sinusoidal components, performing component extraction using a weighted log-likelihood function with the first weighting sequence comprises filtering the samples with a respective component extraction filter tuned to a respective one of the first input set of estimated complex sinusoidal components.

[0034] In some embodiments, performing cross-interference cancellation on the component estimates to produce a new set of cross-interference cancelled component estimates comprises multiplying the component estimates by a cross-interference cancellation matrix.

[0035] In some embodiments, performing complex envelope extraction on the cross-interference cancelled component estimates to produce the new set of complex amplitude estimates comprises multiplying each cross-interference cancelled component estimate by the respective estimated complex sinusoidal component with negative exponent.

[0036] In some embodiments, for each of the second input set of estimated complex sinusoidal components, estimating a frequency deviation estimate using the weighted log-likelihood function with the second weighting sequence comprises filtering the sampled sequence with a respective frequency deviation filter tuned to the estimated complex sinusoidal component.

[0037] In some embodiments, adapting the previous set of frequency estimates taking into account an input set of component estimates and the frequency deviation estimates to produce a new set of frequency estimates comprises applying an adaptation value to each previous frequency estimate, the adaptation value being a function of both the input set of component estimates and the frequency deviation estimates.

[0038] In some embodiments, applying an adaptation value to each previous frequency estimate, the adaptation value being a function of both the input set of component estimates and the frequency deviation estimates comprises: determining a partial derivative with respect to each estimated complex sinusoidal component of a function based on the weighted log-likelihood function; for each frequency estimate, determining the adaptation value from the respective partial derivative.

[0039] In some embodiment, adapting the previous set of frequency estimates taking into account the input set of component estimates and the frequency deviation estimates to produce the new set of frequency estimates comprises: applying an adaptation value to each frequency estimate in the previous set of frequency estimates, the adaptation value being a function of both the component estimates and the frequency deviation estimates to produce an intermediate set of frequency estimates; using the frequency deviation estimates and previous frequency deviation estimates to produce an estimate of chirp for each sinusoidal component; for each sinusoidal component, combining the frequency deviation estimate and the estimate of chirp to produce a new frequency estimate.

[0040] In some embodiments, converting the new set of frequency estimates to new estimated complex sinusoidal components comprises combining previous estimated complex sinusoidal component estimates with the new frequency estimates.

[0041] In some embodiments, combining the previous estimated complex sinusoidal component estimates with the new frequency estimates comprises: multiplying each previous estimated complex sinusoidal component estimate by e{circumflex over ( )}(j×new frequency estimate).

[0042] In some embodiments, one or more ASICs (application specific integrated circuit) adapts to implement a method.

[0043] In some embodiments, one or more DSPs (digital signal processors) adapts to implement a method.

[0044] In some embodiments, one or more FPGAs (field programmable gate arrays) adapts to implement a method.

[0045] In some embodiments, one or more general purpose processors adapts to implement a method.

[0046] In some embodiments, a combination of at least two circuits selected from a group consisting of ASIC, FPGA, DSP, and general purpose processor adapts to implement a method.

[0047] In some embodiments, a computer readable medium having executable code embodied therein for causing a processing platform to execute a method.

[0048] According to one broad aspect, the invention provides an apparatus for tracking amplitude, phase and frequency of a plurality of sinusoidal components in a signal, the apparatus comprising: a first processing path adapted to process the signal to produce a new set of amplitude and phase estimates using a weighted likelihood method; and a second processing path adapted to process the signal to produce a new set of frequency estimates using a weighted likelihood method.

[0049] In some embodiments, the apparatus further comprises: a sampler adapted to sample the signal to produce a sequence of real-valued samples, wherein the first and second processing paths perform signal processing in the digital domain.

[0050] In some embodiments, an apparatus further comprises: a sampler adapted to sample the signal to produce a sequence of complex-valued samples, wherein the first and second processing paths perform signal processing in the digital domain.

[0051] In some embodiments, the first and second processing paths perform signal processing in the continuous time domain.

[0052] According to one broad aspect, the invention provides an apparatus for tracking amplitude, phase and frequency of a plurality of sinusoidal components in a signal, the apparatus comprising: at least one component extraction filter adapted process the signal to produce component estimates for each of a first input set of estimated complex sinusoidal components, each component extraction filter being tuned to a respective one of the first input set of estimated complex sinusoidal components; at least one frequency deviation filter adapted to process the signal to produce a frequency deviation estimate for each of a second input set of estimated complex sinusoidal components, each frequency deviation filter being tuned to a respective one of the second input set of estimated complex sinusoidal components; at least one adaptive frequency tracker adapted to produce a new set of frequency estimates by adapting a previous set of frequency estimates taking into account an input set of component estimates and the frequency deviation estimates; and at least one component generator adapted convert the new set of frequency estimates to a new set of estimated complex sinusoidal components.

[0053] In some embodiments, the signal is a sequence of samples and processing is done in the digital domain, and wherein the at least one component generator comprises at least one digital controlled oscillator.

[0054] In some embodiments, the apparatus further comprises: a cross-interference canceller adapted to perform cross-interference cancellation on the component estimates to produce a new set of cross-interference cancelled component estimates; wherein the new set of cross-interference cancelled estimates are used as the input set of component estimates to the adaptive frequency tracker.

[0055] In some embodiments, the apparatus further comprises: at least one complex envelope estimator adapted to perform complex envelope extraction on the component estimates to produce a new set of complex amplitude estimates.

[0056] In some embodiments, the apparatus further comprises: at least one complex envelope estimator adapted to perform complex envelope extraction on the cross-interference cancelled component estimates to produce a new set of complex amplitude estimates.

[0057] In some embodiments, each component extraction filter implements a weighted log-likelihood function with a first weighting sequence; each frequency deviation filter implements a weighted log-likelihood function with a second weighting sequence.

[0058] In some embodiments, the first and second weighting sequences are the same.

[0059] In some embodiments, the first input set of estimated complex sinusoidal components and the second set of estimated complex sinusoidal components are initially set to initial values, and thereafter are set to previously determined estimated complex sinusoidal components.

[0060] In some embodiments, for each time a new set of complex amplitude estimates is produced by the apparatus: the component extraction filter(s) operate to produce the new set of complex amplitude estimates before the frequency deviation filter(s) operate to produce the new set of estimated complex sinusoidal components; the first input set and the second input set of estimated complex sinusoidal components comprise the new set of estimated complex sinusoidal components determined during a previous update period; wherein the input set of cross-interference cancelled estimates comprises the new set of cross-interference cancelled estimates determined during the current update period.

[0061] In some embodiments, for each time a new set of complex amplitude estimates is produced by the apparatus: the component extraction filter(s) operate to produce the new set of estimated complex sinusoidal components before the frequency deviation filters operate to produce the new set of complex amplitude estimates; the input set of component estimates comprises the set of cross-interference cancelled estimates determined during a previous update period; the first input set of estimated complex sinusoidal components comprises the new set of estimated complex sinusoidal components determined during the current update period and the second input set of estimated complex sinusoidal components comprises the new set of estimated complex sinusoidal components determined during a previous update period.

[0062] In some embodiments, the cross-interference canceller produces the new set of cross-interference cancelled component estimates by multiplying the component estimates by a cross-interference cancellation matrix.

[0063] In some embodiments, the complex envelope estimator(s) produce the new set of complex amplitude estimates by multiplying each cross-interference cancelled component estimate by the respective estimated complex sinusoidal component with negative exponent.

[0064] In some embodiments, the adaptive frequency tracker(s) apply an adaptation value to each previous frequency estimate, the adaptation value being a function of both the component estimates and the frequency deviation estimates.

[0065] In some embodiments, the adaptive frequency tracker(s) determine a partial derivative with respect to each estimated complex sinusoidal component of a function based on a weighted log-likelihood function and for each frequency estimate, determine the adaptation value from the respective partial derivative.

[0066] In some embodiments, the adaptive frequency tracker(s) produce a new set of frequency estimates by applying an adaptation value to each frequency estimate in a previous set of frequency estimates, the adaptation value being a function of both the component estimates and the frequency deviation estimates to produce an intermediate set of frequency estimates, and using the frequency deviation estimates and previous frequency deviation estimates to produce an estimate of chirp for each sinusoidal component, and for each sinusoidal component combine the frequency deviation estimate and the estimate of chirp to produce a new frequency estimate.

[0067] In some embodiments, the component generator(s) convert the new set of frequency estimates to new estimated complex sinusoidal components by combining previous estimated complex sinusoidal component estimates with the new frequency estimates.

[0068] In some embodiments, a computer in combination with a computer readable medium compatible with the computer are provided which are cooperatively adapted to implement any of the above methods.

BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Preferred embodiments of the invention will now be described in greater detail with reference to the accompanying diagrams, in which:

[0070] FIG. 1 is a block diagram of an apparatus for the adaptive estimation and tracking of a multi-component sinusoidal real-valued observed signal provided by an embodiment of the invention;

[0071] FIG. 2 is a flow chart of a method provided by an embodiment of the invention for the adaptive estimation and tracking of a multi-component sinusoidal real-valued observed signal;

[0072] FIG. 3A is a block diagram of a first Component Extraction Filter (CEF) usable in the apparatus of FIG. 1;

[0073] FIG. 3B is a block diagram of a second CEF usable in the apparatus of FIG. 1;

[0074] FIG. 4A is a block diagram of a first Frequency Deviation Filter (FDF) usable in the apparatus of FIG. 1;

[0075] FIG. 4B is a block diagram of a second FDF usable in the apparatus of FIG. 1;

[0076] FIG. 5 is a block diagram showing details of the CIC and CEEs of FIG. 1;

[0077] FIG. 6 is a flow chart of a method for dynamically detecting and updating the number of prominent sinusoidal components to be tracked in the real-valued observed signal;

[0078] FIG. 7 is a block diagram of a speech coder employing the adaptive estimation and tracking method of FIG. 2;

[0079] FIG. 8 is an example implementation of one of the DCOs of FIG. 1;

[0080] FIGS. 9A and 9B contain plots of frequency tracking for different SNRs for a signal with two components. Estimated frequencies (solid) and true frequencies (dotted) are superimposed on background of the spectrum of the main signal, 9A: Tracked frequencies of the first scenario for SNR=0 dB, 9B: SNR=20 dB;

[0081] FIG. 10 contains plots of estimation errors of the proposed algorithm; solid: mean squared error |xn−{circumflex over (x)}n|2 (averaged over 50 runs), dotted: sum of MSE of components |x1,n−{circumflex over (x)}1,n|2+|x2,n−x2,n−{circumflex over (x)}2,n|2;

[0082] FIG. 11 contains plots of estimated amplitudes in 20 dB, using a Hamming window with a length of 129, dotted: true values, solid and dashed: estimated values;

[0083] FIG. 12 contains plots of average of squared frequency estimation error of the proposed algorithm: 1$\frac{{\&LeftBracketingBar;f_1-{\hat{f}}_1\&RightBracketingBar;}^2+{\&LeftBracketingBar;f_2-{\hat{f}}_2\&RightBracketingBar;}^2}{2}$

[0084] (averaged over 50 runs). Solid: one iteration per time instance with μ=10−4, Dashed: two iterations per time instance μ1=10−4 and μ2=5×10−5; and

[0085] FIGS. 13A and 13B contain plots of results of decomposing a segment of speech to four components using a μ sample length Hamming window, two iteration and different p for the components, 13A: Tracked frequencies on background of the spectrogram of the main speech, 13B: Spectrogram of the constructed signal.

DETAILED DESCRIPTION OF THE INVENTION

[0086] The problem to be solved can be conceptualized in a discrete time mathematical model. In such a model, the signal to be examined can be considered a real-valued observed signal xn having L sinusoidal AM-FM components and corrupted by additive white noise. The real-valued observed signal xn consists of a sequence of real-valued samples of a multi-component real-valued signal. The samples are obtained by sampling with a sampling period of T and sampling frequency fs=1/T Hz. The real-valued observed signal can be represented mathematically as follows: 2$\begin{array}{cc}{x}_n=\sum _{l=1}^L\mathrm{Re}\left(a_l{e}^{j{\omega }_{l}n}\right)+N_n;n\in Z& \left(1\right)\end{array}$

[0087] where Nnε is a real additive white noise sample, the complex signal a1=|a16|ejφsi 1ε represents the amplitude and phase of the lth component to be estimated, and ωlε[−π,π] is the frequency (radian/sample) of the lth component to be estimated. Re(.) denotes the real part of a complex number. In equation (1), it is noted that when changing the pair (ζ1, a1) to (−ω1, al*), where (.)* is the complex conjugate operator, the observed signal is unchanged. Thus, the sign of the frequency is not identifiable. Because of this, it can be assumed that ωlε[0,π], and if the method used for estimation results in a negative value for the radian frequency ωl the pair (ωl, al) can be changed to (−ωl, al*) without a loss of any information. The relationship between ωl and the real continuous frequency fs is given by 3$f_l=\frac{{\omega }_l}{2\pi }{f}_s.$

[0088] The problem that is addressed here is to estimate the complex amplitude (i.e. amplitude and phase) and the frequency of all prominent sinusoidal components of such a signal.

[0089] It is noted that an embodiment of the invention also provides a similar but simpler solution for the case of complex observed signals, where 4$x_n=\sum _{l=1}^L{a}_l{e}^{j{\omega }_{l}n}+N_n\in C.$

[0090] The solution for the complex case is a simplified version of the real case solution, since in the complex case the quadrature components of the signal are also observed.

[0091] It is assumed that the amplitudes, the phases and the frequencies of the prominent sinusoidal components are very slowly time-varying or equivalently they are assumed to be band limited and smooth signals. Furthermore, it is assumed that each of the prominent sinusoidal components may disappear or appear, but does so in such a manner that the number of prominent sinusoidal components rarely changes. It is also assumed that the number of prominent sinusoidal components is known initially. The method may still work should one or more of these assumptions fail.

[0092] A likelihood function can be evaluated which represents the amount of information about the received (observed) signal that is available to the receiver. Evaluation of the likelihood function can provide an estimation of the unknown parameters. If it is assumed that Nn is a zero-mean white Gaussian random process with variance σN2(n), the log-likelihood function at time n of the observed xn can be expressed as follows: 5$\begin{array}{cc}L\left(x_n|{\left\{a_n\left(n\right),{\omega }_l\left(n\right),{\sigma }_N^2\left(n\right)\right\}}_{l=1}^L\right)=-\frac{1}{2}\log \left(2\pi {\sigma }_N^2\right)-\frac{{\&LeftBracketingBar;x_n-\sum _{l=1}^L\mathrm{Re}\left(a_l{e}^{j{\omega }_{l}n}\right)\&RightBracketingBar;}^2}{2{\sigma }_N^2}& \left(2\right)\end{array}$

[0093] It is noted that within a stationary environment the maximization of the above likelihood function over a rectangular window (e.g. by expectation maximization method) could provide an appropriate estimate for the characteristic parameters of each of the prominent sinusoids contained in the observed real-valued signal. However, this approach requires a large amount of computation and also suffers from the problem of local optima which can result in inaccurate results.

[0094] Embodiments of the invention provide a method and system for evaluating the likelihood function that is computationally feasible and provides accurate estimates for the characteristic parameters of each of the prominent sinusoids contained in the observed real-valued signal.

[0095] Referring now to FIG. 1, shown is a block diagram of an apparatus provided by an embodiment of the invention for the adaptive estimation and tracking of a multi-component real-valued observed signal xn. It is noted that in what follows, two different indices for time are used. The index “n” is used to refer to the processing performed by the method at current time n. The index “i” is a dummy variable used to refer to the observed signal at times other than the current time n. Typically, the processing at time n uses multiple different observed signals xi.

[0096] A first block of the apparatus, shown in FIG. 1, is a set of Component Extraction Filters (CEFs) 110. There is one CEF for each prominent sinusoidal component being tracked. It is assumed there are L such components where L≧1. The CEFs 110 have a first input 111 and a second input 112. The first input 111 accepts the real-valued observed signal xn, while the second input 112 accepts a set of estimated complex sinusoidal components {{circumflex over (Ω)}in}l=1, . . . , L corresponding to the prominent frequency components in the real-valued observed signal xi being tracked. The CEFs 110 also have an output 113 from which they deliver a set of component estimates Yn made up of estimates of the prominent signal components of the real-valued observed signal xn. Contained within CEFs 110 are a number of filters, preferably band pass filters, each of which is tuned to a respective frequency corresponding to one of the prominent signal components of the real-valued observed signal xn. More specifically, in one embodiment at time n these are L band pass filters having impulse responses {hl,n}l=1,2, . . . L. These band pass filters are described in further detail below.

[0097] A second block of the apparatus is a Cross-Interference Canceller (CIC) 130. The CIC 130 has a first input 131 and a second input 132. The first input 131 accepts the set of component estimates Yn, while the second input 132 accepts the set of estimated complex sinusoidal components {{circumflex over (Ω)}ln}. The CIC 130 also has an output 133 from which it delivers a set of cross-interference cancelled component estimates {circumflex over (X)}n of the prominent signal components of the real-valued observed signal xn.

[0098] A third block of the apparatus is a set of Complex Envelope Estimators (CEEs) 150. The CEEs 150 have a first input 151 and a second input 152. The first input 151 accepts the set of cross-interference cancelled signal component estimates {circumflex over (X)}n, while the second input 152 accepts the set of estimated complex sinusoidal components {{circumflex over (Ω)}ln}. The CEEs 150 also have an output 153 from which they deliver a set of complex amplitude estimates {âl,n}. The CEEs 150 are made up of a number of signal multipliers, each of which is usable for a respective frequency corresponding to one of the prominent signal components of the real-valued observed signal xn.

[0099] A fourth block of the apparatus is a set of Frequency Deviation Filters (FDFs) 120. The FDFs 120 have a first input 121 and a second input 122. The first input 121 accepts the real-valued observed signal xn, while the second input 122 accepts the set of estimated complex sinusoidal components The FDFs 120 also have an output 123 from which they deliver a set of frequency deviation estimates {overscore (Y)}n each a measure of a frequency deviation of a prominent signal component of the real-valued observed signal xn. The FDFs 120 consist of a set of L filters, each of which is tuned to a unique frequency corresponding to one of the prominent signal components of the real-valued observed signal xn. The filters have impulse responses {{overscore (h)}l,n}l=1,2, . . . , L.

[0100] A fifth block of the apparatus is a set of Digital Controlled Oscillators (DCOs) 140. The DCOs 140 have an input 141 and an output 142. The input 141 accepts a set of frequency estimates {{circumflex over (ω)}l,n}l=1L corresponding to the frequencies of prominent frequency components contained in up the real-valued observed signal xn. The output of the DCOs 140 available at the output 142 is a new set of estimated complex sinusoidal components {{circumflex over (Ω)}ln} produced by combining the previous estimated complex sinusoidal components {{circumflex over (Ω)}ln−1} with the frequency estimates {{circumflex over (ω)}l,n}.

[0101] A sixth and last block of the apparatus that is shown in FIG. 1 is a set of Adaptive Frequency Trackers (AFTs) 160. The AFTs 160 have a first input 161 and a second input 162. The first input 161 accepts the set of frequency deviation estimates {overscore (Y)}n, while the second input 162 accepts the set of cross-interference cancelled component estimates {circumflex over (X)}n. The AFTs 160 also have an output 163 from which they deliver the set of frequency estimates {{circumflex over (ω)}l,n}.

[0102] It is noted for the single component case, only one CEF, CEE, FDF, DCO ad AFT are required and there is a single CIC so long as there are two or more components, no CIC being required for the single component case.

[0103] In operation, the real-valued observed signal xn is simultaneously fed to the CEFs 110 and the FDFs 120 via inputs 111 and 121 respectively. The adaptive joint estimation and tracking method provided by the invention is recursive, and the method may equivalently start at either the CEFs 110 or the FDFs 120. In alternate half-iterations of the method, either the complex amplitude estimates {âl,n} are updated as a function of the observed signal xn and previous estimated complex sinusoidal components {{circumflex over (Ω)}ln} or the estimated complex sinusoidal components {{circumflex over (Ω)}ln} and the frequency estimates {{circumflex over (ω)}l,n} are updated as a function of the observed signal xn and previous cross-interference cancelled estimates {circumflex over (X)}n. The CEFs 110 and the FDFs 120 generate values for the component estimates Yn and frequency deviation estimates Ŷn respectively. Both the CEFs 110 and the FDFs 120 are initialized with estimates of the estimated complex sinusoidal components, corresponding to each of the prominent sinusoidal components contained in the real-valued observed signal xn. Furthermore, both the CEFs 110 and the FDFs 120 indirectly supply the DCOs 140 a feedback signal that allows the DCOs 140 to update the frequency estimates.

[0104] The CEFs 110, CIC 130 and CEEs 150 collectively process the observed signal xn to produce complex amplitude estimates âl,n of the prominent sinusoidal components. Each of the filters in CEFs 110 filters the observed signal xn to produce a respective initial estimate Yn of each frequency component. The filters are tuned to look at frequencies specified by the estimated complex sinusoidal components {{circumflex over (Ω)}ln}. The CIC 130 accepts the component estimates Yn generated by the CEFs 110 and the estimated complex sinusoidal components generated by the DCOs 140. The CIC 130 can be basically described as a matrix processor that combines the estimated complex sinusoidal components with the component estimates Yn to produce the cross-interference cancelled component estimates {circumflex over (X)}n. The mathematical details of this block will be given in detail in what follows.

[0105] The CEEs 150 operate by multiplying the estimated complex sinusoidal components and the corresponding cross-interference cancelled component estimates {circumflex over (X)}n to produce the set of complex amplitude estimates, each complex amplitude estimate corresponding to a respective prominent sinusoidal component contained in the observed signal xn.

[0106] The FDFs 120, DCOs 140 and AFTs 160 collectively process the observed signal xn to produce the estimated complex sinusoidal component {{circumflex over (Ω)}ln} and frequency estimates {{circumflex over (ω)}l,n}. The FDFs 120 generate estimates of the prominent sinusoidal components frequency deviations from the real-valued observed signal and previous estimates of the frequencies present in the real-valued observed signal xn. The filters in the FDFs 120 filter the observed signal to produce the set of estimates Ŷn of deviations in the modulated frequency of the prominent sinusoidal components from the previous estimates. The FDFs 120 pass the frequency deviation estimates Ŷn they have generated to the AFTs 160 as shown in FIG. 1. The AFTs 160 use the frequency deviation estimates Ŷn from the FDFs 120 and the cross-interference cancelled component estimates {circumflex over (X)}n from the CIC 130 to generate the set of frequencies {{circumflex over (ω)}l,n} corresponding to the respective components contained the real-valued observed signal.

[0107] The frequencies {{circumflex over (ω)}l,n} generated by the AFTs 160 are output by the method as the new set of frequency estimates. They are also passed to the DCOs 140, which contain recursive complex sinusoidal signal generators that modulate the frequency estimates by combining them with previous estimates to produce estimated complex sinusoidal components and in so doing reduces the effect of noise in the frequency estimates. Filtering the frequency estimates in this way minimizes the amount of computational error propagated by erroneous frequency estimation. As previously indicated the DCOs 140 then send the new estimated complex sinusoidal components {{circumflex over (Ω)}in} to the CEFs 110, the FDFs 120, the CIC 130 and the CEEs 150, so that the next iteration of processing can proceed.

[0108] The mathematics behind the approach embodied in FIG. 1 will now be presented. The adaptive estimation and tracking method provided by an embodiment of the invention is obtained by maximizing a weighted average of equation (2), the likelihood equation given above. The weighted-average of equation (2) is taken over a finite amount of time, or rather within a window in time. A recursive adaptive filter with a low order of computational complexity is employed to achieve this.

[0109] By appropriately selecting the averaging window the problem of converging on an erroneous local optima (e.g. what is known as a false-lock in a PLL) can be overcome, with the additional benefit of controlling the noise limiting and tracking behavior of the adaptive estimation and tracking method. In order to estimate the parameters characterizing prominent sinusoidal components comprising the real-valued observed signal xn at a time n, in a preferred embodiment the following parametric windowed (weighted) log-likelihood function is maximized: 6$\begin{array}{cc}L\left(n\right)=\sum _{i=-\infty }^{+\infty }{W}_{n-i}L\left(x_i|{\left\{a_l\left(i\right),{\omega }_l\left(i\right),{\sigma }_N^2\left(i\right)\right\}}_{l=1}^L\right)& \left(3\right)\end{array}$

[0110] In equation (3) Wk is a window function described in detail below, where k is an index for the window function which is set relative to n, i.

[0111] This approach will be referred to herein as the Maximum Weighted Likelihood (MWL) approach. Despite being presented within a system that assumes a stationary environment, it should be noted that this approach is also suitable for non-stationary and adaptive parameter estimation.

[0112] Furthermore, the weighted log-likelihood function can be employed, and adapted for implementation using the apparatus of FIG. 1, by defining the impulse response {hl,n} and {{overscore (h)}l,n} as described below. In another embodiment, the structure of FIG. 1 is provided with no specific constraints on the impulse responses {hl,n} and {{overscore (h)}l,n}.

[0113] Referring to equation (3), the value wn−i is the weight of the information received at time n−i in order to estimate the unknown parameters at time n. In some embodiments, the window function, wk, is selected to satisfy one or more of the following conditions:

[0114] 1. The window is normalized (i.e. Σkwk=1).

[0115] 2. The window is centered at zero (i.e. ΣkWk=0).

[0116] 3. The causality of the adaptive estimation and tracking method relies on ∀k≦D wk=0, where D denotes an acceptable delay for the estimated parameters characterizing the prominent sinusoidal components comprising the real-valued observed signal xn.

[0117] 4. The condition for L(n) to yield a positive integration from the information gathered at different time instances is that: Wk≧0.

[0118] The first condition is not absolutely necessary. It has been included because it simplifies the discussion in regard to the adaptive estimation and tracking method being described. The second condition could also be violated, in such a case, the value of 7$\sum _{k=-\infty }^{+\infty }k{W}_k$

[0119] represents the delay lag of the estimation. Again the second condition has been included because it simplifies the upcoming description of the adaptive estimation and tracking method. More generally, if one or more of the conditions are not satisfied, the performance might degrade depending on the situation.

[0120] Considering more carefully the above conditions, it can be observed that if wk is considered to be the unit impulse response of a linear time-invariant filter, this filter is strictly a lowpass filter having its maximum frequency domain gain for the DC component.

[0121] Additionally, if the window weights wk are all non-negative, L(n) represents a measure of the received information about the parameters in the neighborhood of time n. This window will determine the resolution of the frequency estimation and the extracting filters for the complex amplitude al.

[0122] The maximization of the criterion in equation (3) does not resolve the problem because the number of unknown parameters in equation (3) is more than the number of parameters obtained by observation within the window. So as an approximation, it is assumed that the unknown parameters are approximately constant across the course of the window and this approximation allows the likelihood function to be simplified. Consequently, the following simplified likelihood function is maximized in place of equation (3): 8$\begin{array}{cc}\begin{array}{c}\hat{L}\left(n\right)=\sum _{i=-\infty }^{+\infty }{w}_{n-i}L\left(x_i|{\left\{a_l\left(n\right),{\omega }_l\left(n\right),{\sigma }_N^2\left(n\right)\right\}}_{l=1}^L\right)\approx L\left(n\right)\\ =-\frac{1}{2}\log \left(2\pi {\sigma }_N^2\left(n\right)\right)-\frac{\sum _{i=-\infty }^{+\infty }{\&LeftBracketingBar;x_i-\sum _{l=1}^L\mathrm{Re}\left(a_l{e}^{j{\omega }_{l}i}\right)\&RightBracketingBar;}^2{w}_{n-i}}{{\sigma }_N^2\left(n\right)}\end{array}& \left(4\right)\end{array}$

[0123] Due to the approximation made above an inherent approximation error is introduced resulting from ignoring information about the real-valued observed signal. The approximation error reduces as the variation of the observed parameters reduces as caused by shortening the window length. Yet with a longer window length a longer duration of the observed signal can be used. This has the effect of increasing the estimation accuracy as long as the variance of the observed parameters is sufficiently small. However, it is not always the case that the variances of the parameters characterizing the prominent sinusoidal components of the real-valued observed signal will be sufficiently small and such circumstances can not be guaranteed to occur in practical applications of the invention. This shows a trade-off to be considered when choosing the duration of the time domain window. The length of the time domain window is preferably chosen so that the effects of the varying parameters characterizing the observed (received) signal are balanced against the requirement to observe an adequate portion of the observed signal. Thus, it is clear that a shorter time domain window length is preferable in practical embodiments of the invention. The exact length of the time domain window may be determined for a particular application by empirical methods.

[0124] Maximizing {circumflex over (L)}(n) with respect to ΓN2(n), the result is the following equation that approximates the variance: 9$\begin{array}{cc}{\hat{\sigma }}_N^2\left(n\right)=\sum _{i=-\infty }^{+\infty }{\&LeftBracketingBar;x_i-\sum _{l=1}^L\mathrm{Re}\left(a_l{e}^{{\mathrm{j\omega }}_{l}i}\right)\&RightBracketingBar;}^2{w}_{n-i}& \left(5\right)\end{array}$

[0125] Next, the value of {circumflex over (σ)}N2(n) is substituted from equation (5) into equation (4) and then {circumflex over (L)} is maximized with respect to other parameters. The result of this operation is another estimate of the likelihood function as follows: 10$\begin{array}{cc}\hat{L}\left(n\right)=-\frac{1}{2}\log \left(2e\pi \sum _{i=-\infty }^{+\infty }{\&LeftBracketingBar;x_i-\sum _{l=1}^L\mathrm{Re}\left(a_l{e}^{j{\omega }_{l}i}\right)\&RightBracketingBar;}^2{w}_{n-i}\right)& \left(6\right)\end{array}$

[0126] Maximizing equation (6) yields to minimization of the following Weighted-Least-Square (WLS) criterion for estimation of amplitude and frequency at time n: 11$\begin{array}{cc}J\left(n\right)=\sum _{i=-\infty }^{+\infty }{\&LeftBracketingBar;x_{n-i}-\sum _{l=1}^L\mathrm{Re}\left(a_l{e}^{{\mathrm{j\omega }}_l\left(n-i\right)}\right)\&RightBracketingBar;}^2{w}_i& \left(7\right)\end{array}$

[0127] The window wk has a direct impact on the lock-in range of the adaptive estimation and tracking method, wherein the main-lobe width of the Fourier spectrum of wk defines the lock-in range. In general the lock-in range of a method is the maximum initial frequency deviation that can be tolerated by the method such that the method can acquire and begin to track the input frequency.

[0128] Throughout the rest of this disclosure; the following notations and definitions are adopted in order to provide all of the mathematical details of the present invention:

[0129] the vector of the complex modulated amplitudes: 12$X\equiv {\left[a_l{e}^{{\mathrm{j\omega }}_{l}n},\cdots ,a_L{e}^{{\mathrm{j\omega }}_{L}n}\right]}^T=\left[x_{l,n},\cdots ,x_{L,n}\right]=X_{r,n}+{\mathrm{jX}}_{i,n}$

[0130] the estimation of the clean signal: 13${\hat{x}}_n=\sum _{l=1}^L\mathrm{Re}\left({\hat{x}}_{l,n}\right)$

[0131] the estimation error:

e n ≡x n −{circumflex over (x)} n

[0132] the derivative of a real-valued function f(v) with respect to a complex variable v: 14$\frac{\partial f\left(v\right)}{\partial v}\equiv \frac{1}{2}\frac{\partial f\left(v\right)}{\partial \mathrm{Re}\left(v\right)}-\frac{j}{2}\frac{\partial f\left(v\right)}{\partial \mathrm{lm}\left(v\right)}$

[0133] A window: wi≧0 and W(1)=1: 15$W\left(z\right)=\sum _{i=-\infty }^{\infty }{w}_i{z}^{-i}$$\stackrel{\_}{W}\left(z\right)=\sum _{i=-\infty }^{\infty }{\mathrm{iw}}_i{z}^{-i}=-z\frac{d}{dz}W\left(z\right)\stackrel{\_}{W}\left(1\right)=0$ Ω≡[Ω1, . . . , ΩL]T≡[ejωl, . . . , ejωL]5,

[0134] To begin, a detailed derivation of how the complex envelope estimates {âl,n} are computed by blocks 110, 130, 150 of FIG. 1 will be presented.

[0135] To maximize J(n) set 16$\begin{array}{cc}\frac{\partial J\left(n\right)}{\partial a_l}=0& \end{array}$

[0136] for l=1, . . . , L. Thus, 17$\begin{array}{cc}\begin{array}{c}\frac{\partial J\left(n\right)}{\partial a_l}=\sum _{i=-\infty }^{\infty }\left(\sum _{p=1}^L\mathrm{Re}\left(a_p{e}^{j{\omega }_p\left(n-1\right)}\right)-x_{n-i}\right)e^{j\omega l\left(n-i\right)}{W}_i\\ =\sum _{i=-\infty }^{\infty }\sum _{p=1}^L\mathrm{Re}\left(a_p{e}^{j{\omega }_p\left(n-i\right)}\right)e^{j{\omega }_l\left(n-i\right)}{W}_i-\\ \sum _{i=-\infty }^{\infty }{x}_{n-i}{e}^{j{\omega }_l\left(n-i\right)}{W}_i\\ \sum _{i=-\infty }^{\infty }\sum _{p=1}^L\frac{a_p{e}^{j{\omega }_p\left(n-i\right)}+a_p^{*}{e}^{-j{\omega }_p\left(n-i\right)}}{2}{e}^{j{\omega }_l\left(n-i\right)}{W}_i-\\ e^{j{\omega }_{l}n}\sum _{i=-\infty }^{\infty }{x}_{n-i}{e}^{-j{\omega }_{l}i}{W}_i.\end{array}& \left(8\right)\end{array}$

[0137] Given the values of the frequencies {ωl}l=1L, the above system (equation (8)) is solved to provide estimates of the complex amplitudes {ap}p=1L for each of the prominent sinusoidal components. To simplify the solution of the above system, it is useful to define the following set of linear filters: 18$\begin{array}{cc}{h}_{l,n}\equiv e^{j{\omega }_{l}n}{W}_n{H}_l\left(z\right)\equiv \sum _{i=-\infty }^{\infty }{h}_{l,n}{z}^{-1}=W\left(e^{j{\omega }_l}z\right)=W\left({\Omega }_{l}z\right).& \left(9\right)\end{array}$

[0138] Component estimate yl,n (collectively the vector Yn of FIG. 1) can be computed according to: 19$\begin{array}{cc}{y}_{l,n}\equiv x_n\otimes h_{l,n}=\sum _{i=-\infty }^{\infty }{x}_{n-1}{e}^{-j{\omega }_{l}i}{W}_i=e^{j{\omega }_{l}n}\sum _{i=-\infty }^{\infty }{x}_{n-1}{e}^{j{\omega }_l\left(n-i\right)}{W}_i=e^{-j{\omega }_{l}n}\left[\left(e^{j{\omega }_{l}n}{x}_n\right)\otimes W_n\right]& \left(10\right)\end{array}$

[0139] Referring to FIG. 3A and the filters given by equation (9), it can be seen that yl,n can be computed as the output of a linear bandpass filter 200 that is tuned to the estimated centre frequency {circumflex over (ω)}l, for the lth prominent sinusoidal component of the real-valued observed signal. The CEFs 110 of FIG. 1 include one such band pass filter 200 for each component. FIG. 3B shows an equivalent lowpass filter implementation of the filter shown in FIG. 3A. A multiplier 208 multiplies the input xn by ej{circumflex over (ω)}ln. The result is lowpass filtered with a lowpass filter 210 having an impulse response equal to the weighting function wk. The output is multiplied by e−j{circumflex over (ω)}ln with multiplier 212. Both filters, shown in FIGS. 3A and 3B, attenuate all components except the corresponding component and the corresponding component appears at the output with unity gain. In other words the filter is substantially dimensionless for a particular component (Hl(Ωl)=1) when the estimated frequency is accurate. In this manner, it is easy to understand how the adaptive estimation and tracking method operates to reject additive noise. The filter Hl(z) should be designed to pass through the lth component. The bandwidth of Hl(z) equals the bandwidth of W and should be more than the bandwidth of the amplitudes.

[0140] Each filter is only adjusted by one parameter, namely the previous estimated complex sinusoidal component, and the output of each filter will contain some interference from other components arising from non-zero gain in the stop band of each filter. These cross-interferences can be removed by summarizing (8) in the following equation: 20$\begin{array}{cc}\begin{array}{c}\frac{\partial J\left(n\right)}{\partial a_l}=\sum _{p=1}^L\sum _{i=-\infty }^{\infty }\frac{a_p{e}^{j\left({\omega }_l+{\omega }_p\right)\left(n-i\right)}+a_p^{*}{e}^{j\left({\omega }_l{\omega }_p\right)\left(n-i\right)}}{2}{W}_i-\\ e^{j{\omega }_{l}n}{y}_{l,n}\\ =e^{j{\omega }_{l}n}\left(\sum \frac{a_p{e}^{j{\omega }_{p}n}W\left(e^{j\left({\omega }_l+{\omega }_p\right)}\right)+a_p^{*}{e}^{-j{\omega }_{p}n}W\left(e^{j\left({\omega }_l-{\omega }_p\right)}\right)}{2}-y_{l,n}\right).\end{array}& \text{(11)}\end{array}$

[0141] From 21$\frac{\partial J\left(n\right)}{\partial a_l}=0,$

[0142] the following is derived 22$\begin{array}{cc}2y_{l,n}=\sum _{p=1}^L{a}_p{e}^{j{\omega }_{p}n}W\left(e^{j\left({\omega }_l+{\omega }_p\right)}\right)+a_p^{*}{e}^{-j{\omega }_{p}n}W\left(e^{j\left({\omega }_l-{\omega }_p\right)}\right).& \left(12\right)\end{array}$

[0143] To have a matrix form solution, it is useful to define the following: 23$\begin{array}{cc}{Y}_n\equiv {\left[2y_{l,n},\dots ,2y_{L,n}\right]}^T=2Y_{r,n}+2{{\mathrm{jY}}_{i,n}\left[W\left(\underset{\_}{\Omega }\right)\right]}_{l,p}\equiv {W\left({\Omega }_l{\Omega }_p\right)\left[V\left(\underset{\_}{\Omega }\right)\right]}_{l,p}\equiv W\left({\Omega }_l/{\Omega }_p\right).& \left(13\right)\end{array}$

[0144] Note that there is no need for a computational division operation in calculation of W(Ωl/Ωp) as W(Ωl/Ωp)=W(ΩlΩp*). Thus from equation (12) and the definitions given in (13), the linear system can be solved by solving the following complex set of equations:

Y n =W (Ω)Xn+V(Ω)Xn*  (14)

[0145] where (.)* denotes the conjugate operator. Expanding the L-dimensional complex-valued system of (14) into a 2L-dimensional real-valued system, the following is the result: 24$\begin{array}{cc}\left[\begin{array}{c}{Y}_{r,n}\\ Y_{i,n}\end{array}\right]=\left[\begin{array}{cc}\mathrm{Re}\left\{W\left(\underset{\_}{\Omega }\right)+V\left(\underset{\_}{\Omega }\right)\right\}& -\mathrm{lm}\left\{W\left(\underset{\_}{\Omega }\right)-V\left(\underset{\_}{\Omega }\right)\right\}\\ \mathrm{lm}\left\{W\left(\underset{\_}{\Omega }\right)+V\left(\underset{\_}{\Omega }\right)\right\}& \mathrm{Re}\left\{W\left(\underset{\_}{\Omega }\right)-V\left(\underset{\_}{\Omega }\right)\right\}\end{array}\right]\left[\begin{array}{c}{X}_{r,n}\\ X_{i,n}\end{array}\right]=F\left(\underset{\_}{\Omega }\right)\left[\begin{array}{c}{X}_{r,n}\\ X_{i,n}\end{array}\right].& \left(15\right)\end{array}$

[0146] When the frequencies of the prominent sinusoidal components are far enough apart F(Ω) will not be ill conditioned.

[0147] Thus, by applying the inverse of F(Ω) the following solution determines the cross-interference cancelled component estimates 25$\begin{array}{cc}\left[\begin{array}{c}{\hat{X}}_{r,n}\\ {\hat{X}}_{i,n}\end{array}\right]=F^{-1}\left(\underset{\_}{\hat{\Omega }}\right)\left[\begin{array}{c}{Y}_{r,n}\\ Y_{i,n}\end{array}\right]& \left(16\right)\end{array}$

[0148] F−1({circumflex over (Ω)}) removes the cross-interference of adjacent prominent sinusoidal components and will be referred to as the “cross-interference cancelled matrix”, or CIC matrix. Using (16), the cross-interference cancelled modulated component estimates (i.e. the vector {circumflex over (X)}n={circumflex over (X)}r,n+j{circumflex over (X)}i,n) are obtained. As J(n) is a quadratic function of the amplitudes the above solution is the optimum.

[0149] The complex amplitudes of the prominent sinusoidal components are then estimated by:

 n =[â l,n ,â 2,n , . . . ,â L,n ]T=[e−j{circumflex over (ω)}1n, . . . , e−j{circumflex over (ω)}2n, . . . e−j{circumflex over (ω)}Ln]T*{circumflex over (X)}n  (17)

[0150] where the symbol * denotes element-wise multiplication of the elements of two vectors. Referring now to FIG. 1, the CIC 130 operates to multiply the components output by the CEFs 110 by the CIC matrix to produce the cross-interference cancelled components {circumflex over (X)}n. The CEEs 150 calculate equation (17) and output the complex envelope estimates âl,n. The operation of the CIC 130 and CEEs 150 together is shown in further detail in FIG. 5, which shows the component estimates before CIC Yn={Yl,n} input to a cross-interference cancellation matrix multiplication which multiplies Yn by F−1({circumflex over (Ω)}) The output {circumflex over (X)}n={{circumflex over (x)}l,n} is multiplied by the set of e−j{circumflex over (ω)}ln in the multipliers 152 of the CEEs 150 to produce the complex amplitude estimates.

[0151] The optimization of J(n) to track and to estimate the frequencies of the prominent sinusoidal components involves non-linear time varying filtering. According to a preferred embodiment of the invention, an adaptive method for the estimation and tracking of frequencies and amplitudes is provided.

[0152] It is assumed that an accurate initial value is provided for the frequency of each prominent sinusoidal component. This initial value can be obtained from the previous iteration of the method, by another method or by a simple initialization procedure. Then, once {circumflex over (X)}n is obtained from (16) it is substituted into J(n) as shown in the following: 26$\begin{array}{cc}{e}_{n-i}=x_{n-i}-\sum _{l=1}^L\mathrm{Re}\left({\hat{a}}_l{e}^{j{\hat{\omega }}_l\left(n-i\right)}\right)=x_{n-i}-\sum _{l=1}^L\mathrm{Re}\left({\hat{x}}_{l,n}{e}^{-j{\hat{\omega }}_{l}i}\right)J\left(n\right)=\sum _{i=-\infty }^{\infty }{\&LeftBracketingBar;x_{n-i}-\sum _{l=1}^L\mathrm{Re}\left({\hat{x}}_{l,n}{e}^{-j{\omega }_{l}i}\right)\&RightBracketingBar;}^2{W}_i& \left(18\right)\end{array}$

[0153] Following the above substitution yielding equation (18) J(n) is minimized with respect to the unknown frequencies, each of the unknown frequencies corresponding to each of the prominent sinusoidal components. Once this is completed, an iterative approach may be used to find {circumflex over (X)}n again and so on until the desired level of accuracy is obtained. However, for this example, only a single iteration for each time step n will be demonstrated. The adaptive methods work with only one iteration. Increasing the number of iterations provides some improvement at the expense of increased computational complexity. The frequency estimates obtained by the adaptive method for the estimation and tracking method is then used as the initial value (estimate) for the next time iteration n+1.

[0154] To optimize, one may ignore the dependency between {circumflex over (X)}n and Ω and compute 27$\frac{\partial J\left(n\right)}{\partial \underset{\_}{\Omega }}$

[0155] as follows 28$\begin{array}{cc}\frac{\partial J\left(n\right)}{\partial {\omega }_l}=-\sum _{i=-\infty }^{\infty }\left(x_{n-i}-\sum _{p=1}^L\mathrm{Re}\left({\hat{x}}_{p,n}{e}^{-j{\hat{\omega }}_{p}i}\right)\right)\mathrm{lm}\left({\hat{x}}_{l,n}{e}^{-j{\hat{\omega }}_{l}i}{W}_i\right)& \left(19\right)\end{array}$

[0156] Defining 29$\frac{\partial J\left(n\right)}{\partial {\omega }_l}$

[0157] to be equal to Im(Δl,n), Δl,n equals 30$\begin{array}{cc}{\Delta }_{l,n}=-\sum _{i=-\infty }^{\infty }\left(x_{n-i}-\sum _{p=1}^L\mathrm{Re}\left({\hat{x}}_{l,n}{e}^{-j{\hat{\omega }}_{l}i}\right)\right){\hat{x}}_{l,n}{e}^{-j{\hat{\omega }}_{l}i}{\mathrm{iW}}_i=-{\hat{x}}_{l,n}\sum _{i=-\infty }^{\infty }\left(x_{n-i}{e}^{-j{\hat{\omega }}_{l}i}{\mathrm{iW}}_i-\sum _{p=1}^L\mathrm{Re}\left({\hat{x}}_{p,n}{e}^{-j{\hat{\omega }}_{p}i}\right)e^{-j{\hat{\omega }}_{l}i}{\mathrm{iW}}_i\right)& \left(20\right)\end{array}$

[0158] To simplify the realization of the above system, and to convert it to matrix form, a set of linear filters is defined similar to those that were used to derive estimates for the complex amplitudes of the prominent sinusoidal components. 31$\begin{array}{cc}{\stackrel{\_}{h}}_{l,n}\equiv e^{-j{\omega }_{l}n}{\mathrm{nW}}_n,{\stackrel{\_}{H}}_l\left(z\right)\equiv \sum _{l=-\infty }^{\infty }{\stackrel{\_}{h}}_{l,n}{z}^{-1}=\stackrel{\_}{W}\left(e^{j{\hat{\omega }}_l}z\right)=\stackrel{\_}{W}\left({\Omega }_{l}z\right).& \left(21\right)\end{array}$

[0159] Estimates of the frequency deviations {overscore (Y)}l,n (collectively the vector {overscore (Y)}n of FIG. 1) can be computed according to: 32$\begin{array}{cc}{\stackrel{\_}{y}}_{l,n}\equiv x_n\otimes {\stackrel{\_}{h}}_{l,n}=\sum x_{n-i}{e}^{-j{\hat{\omega }}_{l}i}i{w}_i.& \left(22\right)\end{array}$

[0160] It is now useful to define the following: 33$\begin{array}{cc}\begin{array}{c}{\stackrel{\_}{Y}}_n={\left[{\stackrel{\_}{y}}_{1,n},\dots ,{\stackrel{\_}{y}}_{L,n}\right]}^T,\\ {\left[\stackrel{\_}{W}\left(\widehat{\underset{\_}{\Omega }}\right)\right]}_{l,p}\equiv \stackrel{\_}{W}\left({\hat{\Omega }}_l{\hat{\Omega }}_p\right),\\ {\left[\stackrel{\_}{V}\left(\widehat{\underset{\_}{\Omega }}\right)\right]}_{l,p}\equiv \stackrel{\_}{W}\left({\hat{\Omega }}_l/{\hat{\Omega }}_p\right),\\ {\Delta }_n\equiv {\left[{\Delta }_{1,n},{\Delta }_{2,n},\dots ,{\Delta }_{L,n}\right]}^T.\end{array}& \left(23\right)\end{array}$

[0161] Thus, 34$\begin{array}{cc}\begin{array}{c}{\Delta }_{l,n}={\hat{x}}_{l,n}\left(\sum _{p=1}^L\sum _{i=-\infty }^{\infty }\mathrm{Re}\left({\hat{x}}_{p,n}{e}^{-j{\hat{\omega }}_{p}i}\right)e^{-j{\hat{\omega }}_{l}i}i{w}_i-{\stackrel{\_}{y}}_{l,n}\right)\\ ={\hat{x}}_{l,n}\sum _{p=1}^L\sum _{i=-\infty }^{\infty }\frac{1}{2}\left({\hat{x}}_{p,n}{e}^{-j{\omega }_{p}i}+{\hat{x}}_{p,n}^{*}{e}^{j{\omega }_{p}i}\right)e^{-j{\hat{\omega }}_{l}i}i{w}_i-{\hat{x}}_{l,n}{\stackrel{\_}{y}}_{l,n}\\ ={\hat{x}}_{l,n}\left[\sum \frac{1}{2}\left({\hat{x}}_{p,n}\stackrel{\_}{W}\left(e^{j\left({\hat{\omega }}_l+{\hat{\omega }}_p\right)}\right)+{\hat{x}}_{p,n}^{*}\stackrel{\_}{W}\left(e^{j\left({\hat{\omega }}_l-{\hat{\omega }}_p\right)}\right)\right)-{\stackrel{\_}{y}}_{l,n}\right]\end{array}& \left(24\right)\end{array}$

[0162] From (23) and (24), we find 35$\begin{array}{cc}{\Delta }_n={\hat{X}}_n\otimes \left(\frac{1}{2}\lfloor \stackrel{\_}{W}\left(\widehat{\underset{\_}{\Omega }}\right){\hat{X}}_n+\stackrel{\_}{V}\left(\widehat{\underset{\_}{\Omega }}\right){\hat{X}}_n^{*}\rfloor -{\stackrel{\_}{Y}}_n\right)& \left(25\right)\end{array}$

[0163] The vector Δn can be computed directly using the previous estimated complex sinusoidal components {{circumflex over (Ω)}ln}, the cross-interference cancelled component estimates {circumflex over (X)}n, and the frequency deviation estimates {overscore (Y)}n.

[0164] Different methods of applying the frequency deviation estimates to determine new frequency estimates may be applied. In one embodiment the vector Δn is used to produce new estimates for the frequencies. Referring to FIG. 1, the AFTs 160 compute the Δn from {circumflex over (X)}n, Ŷn and the previous estimated complex sinusoidal components defined by {{circumflex over (Ω)}ln}. The AFTs then apply the values of Δn to compute new frequencies. The following equation defines a simple gradient method that is preferably used in the present embodiment of the invention:

{circumflex over (ω)}l,n+1={circumflex over (ω)}l,n−μlmΔl,n; l=1, . . . , L  (26)

[0165] It is to be understood other methods may be used. In equation (26) μ is a very small-sized positive adaptation step that is usually a constant value, i.e. it does not vary in time. One may apply the result of (26) into (16) and vice versa recursively. However, in the simplest implementation, (16) is evaluated and then (26) is evaluated, each once at each time instant (step) n. The overall method is an adaptive method for the estimation and tracking of small frequency variations and a least square estimation of amplitude of the prominent sinusoidal components.

[0166] Referring to FIG. 4A and the filters defined by equation (23), in some embodiments {overscore (y)}l,n can be implemented as the output of a linear bandpass filter 300 that is tuned to the estimated centre frequency {circumflex over (ω)}l, for the lth prominent sinusoidal component of the real-valued observed signal. FIG. 4B shows the substantially equivalent lowpass filter implementation 310 of the filter shown in FIG. 4A. The structure of the lowpass filter embodiment is similar to that of FIG. 3B, but with a different lowpass filter impulse response. A third alternative to FIGS. 3A, 4A and FIGS. 3B, 4B, is to implement filters in an Intermediate Frequency (IF) (for both digital and the analogue implementation described below). For analogue implementations of the algorithm the IF implementation could be advantageous.

[0167] In order to avoid computational error and to reduce the effect of noise on the generation of the estimation of the prominent sinusoidal components Ωln=ejωln the DCOs 140 are used. The DCOs contain recursive complex sinusoidal signal generators which operate on previous estimated complex sinusoidal components, and the new frequency estimates {circumflex over (ω)}l,n as follows:

{circumflex over (Ω)}ln={circumflex over (Ω)}ln−1ej{circumflex over (ω)}l,n  (27)

[0168] The output of the DCOs 140 is the new set of estimated complex sinusoidal components.

[0169] An example DCO implementation is shown in FIG. 8. The estimated frequency {circumflex over (ω)}l,n is converted to a complex exponential ej{circumflex over (ω)}l,n by a non-linear function 800. The previous output of the DCO is {circumflex over (Ω)}ln−i, and this is subject to one unit delay 802 so as to be made available at the current processing instant. The current complex sinusoidal {circumflex over (Ω)}ln is determined by multiplying ej{circumflex over (ω)}l,n by {circumflex over (Ω)}ln−i with multiplier 804. This structure is repeated for each prominent sinusoidal component.

[0170] It is noted that in a conventional FM radio receiver, the input signal has only one component and the amplitude is assumed to be constant. Consequently, a structure with behavior similar to a single Frequency Deviation Filter (FDF) 120 is used for demodulation of frequency. However, the demodulated signal is proportional to the deviation of its input frequency and is not exactly equal to the actual frequency. In the approach provided by the present invention, this signal is used in a feedback loop as in FIG. 1 and described by equation (22), to adaptively estimate/correct the frequency by minimizing the deviation of the observed signal from its estimate. Because of this, it is expected that the adaptive method for the estimation and tracking will provide improvement even in a FM radio.

[0171] Referring now to FIG. 2, shown is a flow chart that illustrates the adaptive estimation and tracking method implemented by the functional blocks of FIG. 1. Included in the flow chart is a preferred initialization step that may be used to start the process at its very first iteration.

[0172] The adaptive estimation and tracking method has an initialization stage 650 and a steady-state operation stage 660.

[0173] The initialization stage 650 begins at step 2-1 with choosing an incremental step size p to update the instantaneous frequencies of the prominent sinusoidal components. The method continues at step 2-2 with the selection and initialization of a causal window function through which to observe the real-valued signal. The window function is chosen so that it satisfies the aforementioned conditions. An initial estimation of the number of prominent sinusoid components comprising the received signal is made in step 2-3. This is then followed by step 2-4 in which the frequencies are initialized. The first three steps of the initialization stage 2-1, 2-2 and 2-3 may be done in any order; however, step 2-4 can only be done after 2-3 has been completed.

[0174] The first step of the steady state operation 660 of the method is to observe the real-valued signal, as indicated in step 2-5 of FIG. 2. In this case it is assumed that the iterative process begins with the computation of the complex amplitudes. However, it is noted that equivalently, the process could begin with the computation of the frequency estimates. The method continues at step 2-6 with the performance of the component extraction step. At step 2-7 the cross-interference cancellation step is performed. At step 2-8, the complex envelopes are extracted. The output of step 2-8 is a new set of amplitude estimates. Next, at step 2-9 frequency deviation filtering is performed. Step 2-9 occurs on the basis of the same set of observed signals xn as was used in step 2-6. The method continues at step 2-10 with the performance of the adaptive frequency tracking. Finally, the iteration finishes at step 2-11 with the updating of the frequency estimates with the digital control oscillators. Depending upon whether or not the entire process is to be iterated, as determined by step 2-11, all of steps 2-6 through 2-11 can be repeated for the same set of observed signals. In the simplest implementation, only one iteration is performed. The method continues back at step 2-5 with the observation of the next real-valued signal xn. The updating of the complex amplitude estimates is one half-iteration, and the updating of the frequencies is another half-iteration.

[0175] Preferably the invention is further enhanced so as to be able to simultaneously detect the appearance of new prominent sinusoidal components, the presence of already identified and previously tracked prominent sinusoidal components and the disappearance of prominent sinusoidal components comprising the real-valued observed signal.

[0176] Referring to FIG. 6, the method begins at step 6-1 with the assumption that at the previous time instant n−1, L prominent sinusoidal components have been detected and are currently being tracked. The problem is further divided into the following test of hypotheses:

[0177] 1) a previously tracked prominent sinusoidal component is still present or has diminished in energy to the point where it can no longer be considered a prominent component of the real-valued observed signal. In the latter case, the once-prominent sinusoidal component shall be considered to no longer exist within the real-valued observed signal and its corresponding parameters will be ignored (dropped) and no longer updated;

[0178] 2) if a new prominent sinusoidal component has been detected within the real-valued observed signal, its frequency shall be initialized.

[0179] Thus step 6-2 is to compare the estimated energy within bands across the periodgram spectrum with a threshold. To calculate the periodgram, first over a frame of signal samples the Fast Fourier Transform of the input signal is calculated, and then the periodgram equals the squared absolute value of the results. The threshold is preferably proportional with {circumflex over (σ)}N2(n) The estimate of noise energy, {circumflex over (σ)}N2(n), is obtained by applying the signal 36${\&LeftBracketingBar;x_{n-i}-\sum _{l=1}^L\mathrm{Re}\left(a_l{e}^{{\mathrm{j\omega }}_l\left(n-i\right)}\right)\&RightBracketingBar;}^2$

[0180] to the lowpass filter wi, as in (7). This allows the identification of all components in the real-valued signal which have an energy above the threshold.

[0181] In step 6-3 any newly-identified prominent sinusoidal components have their characteristic parameters initialized within the adaptive estimation and tracking method, as previously discussed.

[0182] If the energy of a previously tracked prominent sinusoidal component compared with the noise energy, {circumflex over (σ)}N2(n), is very small, then it means that the corresponding component is no longer a prominent component of the real-valued observed signal and can be ignored by the rest of the method provided by the invention as indicated in step 6-4.

[0183] An additional step 6-5 can be used in the method to improve its performance by estimating the speed of frequency variations. The estimated speeds of the frequencies can be then used to overcome the crossover problem. When the frequencies of two components are not distant enough their separation may not be possible. Consequently, the algorithm may fail to track them properly. For each component, therefore, a test may be made to determine whether it is far enough apart from other frequencies. For each frequency component that is not far enough apart, instead of using the adaptive algorithm, the frequency algorithm can simply assume that the speed of the frequency variation is constant during the crossover.

[0184] The five steps 6-1, 6-2, 6-3, 6-4 and 6-5 of FIG. 6 need to be repeated often enough to keep track of components being added or dropped. The number of signal components might change often or not, and as such the rate of repetition should be selected depending on the application.

[0185] To apply the method to complex observed signals, the cross interference canceller would be a smaller matrix which deals with complex numbers instead, having dimensions half of what are required for real valued observed signals.

[0186] In another embodiment, with slight modifications the above-described methods/apparatuses can be used to estimate the chirp of each sinusoidal. Instead of equation (26) above, the following calculations are performed by the adaptive frequency tracker 160:

{tilde over (ω)}l,n={circumflex over (ω)}l,n−μIm(Δl,n),

{circumflex over (ξ)}l,n=(1−α){circumflex over (ξ)}l,n−1+α({tilde over (ω)}l,n−{tilde over (ω)}l,n−1),

{circumflex over (ω)}l,n+1={tilde over (ω)}l,n+{circumflex over (ξ)}l,n.  (28)

[0187] where ξl,n is the chirp (the speed of variation of the frequency). In the above frequency estimation procedure {circumflex over (ω)}l,n is estimated as in (26). Then the variation of the estimated frequency in two successive time instant, i.e., ({circumflex over (ω)}l,n−{circumflex over (ω)}l,n−1) is applied to a lowpass filter that provides {circumflex over (μ)}l,n as an estimate for the chirp parameter of the corresponding component. Finally, the frequency for next time iteration is predicted by a simple integrator that is the third equation in the above procedure.

[0188] In the case that the frequency varies smoothly with time the above procedure results in considerable improvement in frequency tracking at the expense of negligible computation. Examples of possible applications are wideband FM demodulation, speech frequency tracking, chirp estimation in some tracking radars, and several biomedical applications.

[0189] FIG. 7 is a block diagram of a speech coder provided by an embodiment of the invention. In this example, an input signal is first converted to digital form with analogue-to-digital converter 700 to produce the real-valued observed signal xn. This is processed by an adaptive estimation and tracking function 702 which operates in accordance with one of the previously described embodiments. The output of the adaptive estimation and tracking function 702 is a set of complex amplitudes 704 and a set of estimated frequencies 706. Together, amplitudes 704 and frequencies 706 are fed to a signal coding block 708 which encodes this information either for storage or transmission.

[0190] In the embodiments described, it is assumed that discreet signals are being processed and that the entire implementation is digital, requiring a D/A converter if the original signal is analogue. In another embodiment, the apparatus/method is implemented in the analogue domain directly eliminating any need for A/D conversion. The block diagram of such a system is basically the same as FIG. 1, except that the component extraction filters and frequency deviation filters are continuous time filters having continuous impulse responses. Instead of a discrete time signal function wk, a continuous time function Wt with similar properties should be used. Instead of the digital controlled oscillators 140, quadrature analogue voltage controlled oscillators are employed.

[0191] The adaptive frequency tracker 160 instead of employing a delay unit may use a simple integrator as following instead of (26) 37$\begin{array}{cc}\frac{d}{dt}{\hat{\omega }}_{l,t}=-\mu \mathrm{Im}\left({\Delta }_{l,t}\right).& \left(29\right)\end{array}$

[0192] Other changes are basically to use dual continuous-time elements equivalent to the corresponding digital elements that are described above.

[0193] In the embodiments described, it is assumed that cross-interference cancellation is performed, and this in generally will yield the best performance. In some embodiments of the invention, acceptable results may be realized without including cross-interference cancellation. For example, in processing signals having very little cross-interference, the cross-interference cancellation would provide only a small improvement.

[0194] In the embodiments described, complex envelope extraction is performed at the output to generate new sets of complex amplitude estimates. In some embodiments, complex envelope extraction can be omitted if a complex amplitude output is not required. In these embodiments, all the steps of complex envelope estimation are performed except the last step of extraction, since the output of this step is not fed back into other steps/components of the method/system.

[0195] In the embodiments described, it is assumed the same weighting function is employed for each of frequency tracking and amplitude tracking. More generally, a different weighting function may be employed for each of these purposes in which case a first weighting function is applied for frequency tracking, and a second weighting function is applied for amplitude estimation.

[0196] Simulations and Discussion

[0197] Simulations have been conducted for different signals. Firstly two scenarios are considered in order to study the performance of an example implementation of an embodiment of the invention for a signal with known components. The first scenario deals with the performance in different environments (SNR), for different signals (sinusoids and chirps) and also in the case of a cross-over. The second scenario studies the effect of the initialization and investigates the relationship between LIR and the length of the window. Then, the example implementation is applied to a speech signal as a multi-component stochastic signal, to study the performance for speech signals.

[0198] FIGS. 9A and 9B show the results of the first scenario where a signal with two components is considered. The frequencies of these components are time-dependent with a cross-over. This figure shows the performance of the algorithm in tracking sinusoids and crossing chirp components with constant amplitudes for SNR=0 in FIG. 9A and 20 dB in FIG. 9B. A Hamming window with length 129 was used. Initial points for frequencies are chosen close enough to the true values to ensure initial convergence, and the step-size is μ=10−4. The greater the step-size, the larger the variance of estimated frequencies and the lower the bias of the estimation. FIG. 9 clearly shows that the estimated values converge to the true values. In the chirp section, the algorithm tracks the components with a bias in the estimated frequencies. This bias/lag is a function of the window length; the shorter the window, the smaller the bias. For the chirp signals a shorter window is preferred, while for the sinusoidal signals a longer window works better. Around the cross-over moment, the performance degrades, since the two components are too close to each other so that they pass through the same filter and cannot be separated efficiently. In other words, the CIC is not able to cancel out the interference completely as the matrix F (Ω) becomes ill-conditioned. After the moment of cross-over, when the frequencies are far enough apart, the algorithm recovers the frequencies and tracks them again.

[0199] FIG. 10 supports the same results, where the estimation errors are plotted for SNR=20 dB. The effect of lag in the tracking of the chirp components is reflected in some bias in frequencies and a jump on the level of the MSE in both plots. FIG. 11 shows the magnitude of the estimated amplitudes and the true amplitudes. The estimated complex amplitudes are oscillatory with a very low frequency resulting from the estimation frequency lag, as they are the outputs of demodulators. These oscillations can be effectively smoothed by allying LPFs to the magnitude of the output of the demodulators, assuming a known bandwidth for true amplitudes. FIG. 12 shows the average of the squared frequency estimation error for two different implementations. The first one, the same version used in previous simulations, does the estimation/tracking process once in each time step. The second version employs a second additional iteration for each time step with a smaller step-size of μ2=5×10−5. The algorithm with a second iteration provides a higher accuracy in tracking sinusoidal signals, as the error variance in frequency estimation is as low as 10−7. Due to the tracking lag for the chirp signals, there is a jump in the variance to 10−5. Around the moment of cross-over, the error increases. After the moment of cross-over the error is seen to increase due to the interchange of the frequencies. Using two iterations per time step, one gains shorter convergence period and less lag in chirp tracking at the expense of twice the computation.

[0200] FIGS. 9A, 9B, 10 and 11 show that there is an interchange between the components at the moment of cross-over, which the algorithm does not detect. On the other hand this interchange is reflected in the estimated amplitudes as shown in FIG. 11. Hence, using estimated amplitudes the cross-over moment is detectable and resolvable as long as the amplitudes are different.

[0201] The second scenario is defined in order to study the relationship between the window length and the LIR and was performed in different environments. Extensive simulations show that the LIR does not depend on the SNR (LIR measured with a resolution of 0.0011 fs where fs is the sampling frequency). For a specific window type, a shorter window length causes greater variances in tracking, because fewer signal samples are used in the estimation. At the same time, a shorter window length provides a wider LIR and results in a smaller bias. The tracking bias for the chirp component increases for longer windows, because the assumption of constant frequencies along the support of window becomes invalid.

[0202] For a specific window type, since this window determines the shape of all involved filters, the LIR is inversely proportional to the window length. For instance for a Hamming window the LIR is 0.015 and 0.027, when the window length is 129 and 65, respectively, where the frequencies are normalized with respect to fs. If the initial frequency error is greater than the LIR, then the convergence of the algorithm might take substantial time to fall in the LIR. Once the frequency estimation error is less than the LIR, the algorithm converges with a time constant controlled by the algorithm step-size.

[0203] To study the performance of the algorithm for speech signals, a Hamming window of length 47 samples is used with four components of initial frequencies of 250, 2750, 3500 and 4500 Hz where the sampling frequency was 11025 Hz. The estimation/tracking process is implemented using two iterations per time step and two step sizes μ1,1 and μ2,1=μ1,1/3, respectively. Used step-sizes are μ1,1=0.02 for the first component and 0.04 for all others. FIGS. 13A and 13B show the results obtained for a speech signal. FIG. 13A depicts the tracked frequencies on the spectrum background of the speech signal. Clearly, frequencies are tracking the dominant energy segments of the spectrum. FIG. 13B shows the spectrum of the constructed signal and supports how successful the algorithm is in decomposing stochastic signals.

[0204] What has been described is merely illustrative of the application of the principles of the invention. Other arrangements and methods can be implemented by those skilled in the art without departing from the spirit and scope of the present invention.

Claims

1. A method of tracking amplitude, phase and frequency of a plurality of sinusoidal components in a signal, the method comprising: a) processing the signal to produce a new set of amplitude and phase estimates using a weighted likelihood method; and b) processing the signal to produce a new set of frequency estimates using a weighted likelihood method.

2. A method according to claim 1 further comprising sampling the signal to produce a sequence of real-valued samples, wherein steps a) and b) are performed in the digital domain.

3. A method according to claim 1 further comprising sampling the signal to produce a sequence of complex-valued samples, wherein steps a) and b) are performed in the digital domain.

4. A method according to claim 1 wherein steps a) and b) are performed in the continuous time domain.

5. A method of tracking amplitude, phase and frequency of a plurality of sinusoidal components in a signal, the method comprising: for a current update period: i) processing the signal to produce a new set of complex amplitude estimates by: a) for a first input set of estimated complex sinusoidal components, separating components to produce component estimates; ii) processing the signal to produce a new set of estimated complex sinusoidal components by: b) for each component of a second input set of estimated complex sinusoidal components, estimating a frequency deviation estimate; c) adapting a previous set of frequency estimates taking into account an input set of component estimates and the frequency deviation estimates to produce a new set of frequency estimates; and d) converting the new set of frequency estimates to a new set of estimated complex sinusoidal components.

6. A method according to claim 5 wherein the signal is a sequence of samples and processing is done in the digital domain.

7. A method according to claim 5 wherein the processing is done in the continuous time domain.

8. A method according to claim 5 further comprising: performing cross-interference cancellation on the component estimates to produce a new set of cross-interference cancelled component estimates, and using the new set of cross-interference cancelled estimates as the input set of component estimates in an execution of step c).

9. A method according to claim 5 further comprising: performing complex envelope extraction on the component estimates to produce a new set of complex amplitude estimates.

10. A method according to claim 8 further comprising: performing complex envelope extraction on the cross-interference cancelled component estimates to produce a new set of complex amplitude estimates.

11. A method according to claim 6 wherein: for the first input set of estimated complex sinusoidal components, separating components to produce component estimates is done using a weighted log-likelihood function with a first weighting sequence; for each of the second input set of estimated complex sinusoidal components, estimating the frequency deviation estimate is done using a weighted log-likelihood function with a second weighting sequence.

12. A method according to claim 11 wherein the first and second weighting sequences are the same.

13. A method according to claim 9 wherein step i) comprises a first half-iteration, and step ii) comprises a second half iteration, one first half-iteration and one second half-iteration comprising a complete iteration and wherein for each update period, a plurality of complete iterations are performed to produce the new set of complex amplitude estimates and the new set of estimated complex sinusoidal components.

14. A method according to claim 5 wherein the first input set of estimated complex sinusoidal components and the second set of estimated complex sinusoidal components are initially set to initial values, and thereafter are set to estimated complex sinusoidal components produced by a previous iteration of the method.

15. A method according to claim 10 wherein for each update of the complex amplitude and frequency: the step of processing samples of the sequence of samples to produce a new set of complex amplitude estimates is performed before the step of processing the sequence of samples to produce a new set of estimated complex sinusoidal components; the first input set and the second input set of estimated complex sinusoidal components comprise the new set of estimated complex sinusoidal components determined during a previous update period; wherein the input set of cross-interference cancelled estimates comprises the new set of cross-interference cancelled estimates determined during the current update period.

16. A method according to claim 10 wherein for each update of the amplitude, phase and frequency: the step of processing the signal to produce a new set of estimated complex sinusoidal components is performed before the step of processing the sequence of samples to produce a new set of complex amplitude estimates; the input set of component estimates comprises the set of cross-interference cancelled estimates determined during a previous update period; the first input set of estimated complex sinusoidal components comprises the new set of estimated complex sinusoidal components determined during the current update period and the second input set of estimated complex sinusoidal components comprises the new set of estimated complex sinusoidal components determined during a previous update period.

17. A method according to claim 11 wherein for the first input set of estimated complex sinusoidal components, performing component extraction using a weighted log-likelihood function with the first weighting sequence comprises filtering the samples with a respective component extraction filter tuned to a respective one of the first input set of estimated complex sinusoidal components.

18. A method according to claim 8 wherein performing cross-interference cancellation on the component estimates to produce a new set of cross-interference cancelled component estimates comprises multiplying the component estimates by a cross-interference cancellation matrix.

19. A method according to claim 10 wherein performing complex envelope extraction on the cross-interference cancelled component estimates to produce the new set of complex amplitude estimates comprises multiplying each cross-interference cancelled component estimate by the respective estimated complex sinusoidal component with negative exponent.

20. A method according to claim 11 wherein for each of the second input set of estimated complex sinusoidal components, estimating a frequency deviation estimate using the weighted log-likelihood function with the second weighting sequence comprises filtering the sampled sequence with a respective frequency deviation filter tuned to the estimated complex sinusoidal component.

21. A method according to claim 5 wherein adapting the previous set of frequency estimates taking into account an input set of component estimates and the frequency deviation estimates to produce a new set of frequency estimates comprises applying an adaptation value to each previous frequency estimate, the adaptation value being a function of both the input set of component estimates and the frequency deviation estimates.

22. A method according to claim 21 wherein applying an adaptation value to each previous frequency estimate, the adaptation value being a function of both the input set of component estimates and the frequency deviation estimates comprises: determining a partial derivative with respect to each estimated complex sinusoidal component of a function based on the weighted log-likelihood function; for each frequency estimate, determining the adaptation value from the respective partial derivative.

23. A method according to claim 5 wherein adapting the previous set of frequency estimates taking into account the input set of component estimates and the frequency deviation estimates to produce the new set of frequency estimates comprises: applying an adaptation value to each frequency estimate in the previous set of frequency estimates, the adaptation value being a function of both the component estimates and the frequency deviation estimates to produce an intermediate set of frequency estimates; using the frequency deviation estimates and previous frequency deviation estimates to produce an estimate of chirp for each sinusoidal component; for each sinusoidal component, combining the frequency deviation estimate and the estimate of chirp to produce a new frequency estimate.

24. A method according to claim 5 wherein converting the new set of frequency estimates to new estimated complex sinusoidal components comprises combining previous estimated complex sinusoidal component estimates with the new frequency estimates.

25. A method according to claim 24 wherein combining the previous estimated complex sinusoidal component estimates with the new frequency estimates comprises: multiplying each previous estimated complex sinusoidal component estimate by e{circumflex over ( )}(j×new frequency estimate).

26. One or more ASICs (application specific integrated circuit) adapted to implement a method according claim 1.

27. One or more DSPs (digital signal processors) adapted to implement a method according to claim 1.

28. One or more FPGAs (field programmable gate arrays) adapted to implement a method according to claim 1.

29. One or more general purpose processors adapted to implement a method according to claim 1.

30. A combination of at least two circuits selected from a group consisting of ASIC, FPGA, DSP, and general purpose processor adapted to implement a method according to claim 1.

31. A computer readable medium having executable code embodied therein for causing a processing platform to execute a method according to claim 1.

32. An apparatus for tracking amplitude, phase and frequency of a plurality of sinusoidal components in a signal, the apparatus comprising: a first processing path adapted to process the signal to produce a new set of amplitude and phase estimates using a weighted likelihood method; and a second processing path adapted to process the signal to produce a new set of frequency estimates using a weighted likelihood method.

33. The apparatus according to claim 32 further comprising: a sampler adapted to sample the signal to produce a sequence of real-valued samples, wherein the first and second processing paths perform signal processing in the digital domain.

34. An apparatus according to claim 32 further comprising: a sampler adapted to sample the signal to produce a sequence of complex-valued samples, wherein the first and second processing paths perform signal processing in the digital domain.

35. An apparatus according to claim 32 wherein the first and second processing paths perform signal processing in the continuous time domain.

36. An apparatus for tracking amplitude, phase and frequency of a plurality of sinusoidal components in a signal, the apparatus comprising: at least one component extraction filter adapted process the signal to produce component estimates for each of a first input set of estimated complex sinusoidal components, each component extraction filter being tuned to a respective one of the first input set of estimated complex sinusoidal components; at least one frequency deviation filter adapted to process the signal to produce a frequency deviation estimate for each of a second input set of estimated complex sinusoidal components, each frequency deviation filter being tuned to a respective one of the second input set of estimated complex sinusoidal components; at least one adaptive frequency tracker adapted to produce a new set of frequency estimates by adapting a previous set of frequency estimates taking into account an input set of component estimates and the frequency deviation estimates; and at least one component generator adapted convert the new set of frequency estimates to a new set of estimated complex sinusoidal components.

37. An apparatus according to claim 36 wherein the signal is a sequence of samples and processing is done in the digital domain, and wherein the at least one component generator comprises at least one digital controlled oscillator.

38. An apparatus according to claim 36 further comprising: a cross-interference canceller adapted to perform cross-interference cancellation on the component estimates to produce a new set of cross-interference cancelled component estimates; wherein the new set of cross-interference cancelled estimates are used as the input set of component estimates to the adaptive frequency tracker.

39. An apparatus according to claim 36 further comprising: at least one complex envelope estimator adapted to perform complex envelope extraction on the component estimates to produce a new set of complex amplitude estimates.

40. An apparatus according to claim 38 further comprising: at least one complex envelope estimator adapted to perform complex envelope extraction on the cross-interference cancelled component estimates to produce a new set of complex amplitude estimates.

41. An apparatus according to claim 37 wherein: each component extraction filter implements a weighted log-likelihood function with a first weighting sequence; each frequency deviation filter implements a weighted log-likelihood function with a second weighting sequence.

42. An apparatus according to claim 41 wherein the first and second weighting sequences are the same.

43. An apparatus according to claim 36 wherein the first input set of estimated complex sinusoidal components and the second set of estimated complex sinusoidal components are initially set to initial values, and thereafter are set to previously determined estimated complex sinusoidal components.

44. An apparatus according to claim 38 wherein for each time a new set of complex amplitude estimates is produced by the apparatus: the component extraction filter(s) operate to produce the new set of complex amplitude estimates before the frequency deviation filter(s) operate to produce the new set of estimated complex sinusoidal components; the first input set and the second input set of estimated complex sinusoidal components comprise the new set of estimated complex sinusoidal components determined during a previous update period; wherein the input set of cross-interference cancelled estimates comprises the new set of cross-interference cancelled estimates determined during the current update period.

45. An apparatus according to claim 38 wherein for each time a new set of complex amplitude estimates is produced by the apparatus: the component extraction filter(s) operate to produce the new set of estimated complex sinusoidal components before the frequency deviation filters operate to produce the new set of complex amplitude estimates; the input set of component estimates comprises the set of cross-interference cancelled estimates determined during a previous update period; the first input set of estimated complex sinusoidal components comprises the new set of estimated complex sinusoidal components determined during the current update period and the second input set of estimated complex sinusoidal components comprises the new set of estimated complex sinusoidal components determined during a previous update period.

46. An apparatus according to claim 38 wherein the cross-interference canceller produces the new set of cross-interference cancelled component estimates by multiplying the component estimates by a cross-interference cancellation matrix.

47. An apparatus according to claim 40 wherein the complex envelope estimator(s) produce the new set of complex amplitude estimates by multiplying each cross-interference cancelled component estimate by the respective estimated complex sinusoidal component with negative exponent.

48. An apparatus according to claim 36 wherein the adaptive frequency tracker(s) apply an adaptation value to each previous frequency estimate, the adaptation value being a function of both the component estimates and the frequency deviation estimates.

49. An apparatus according to claim 48 wherein the adaptive frequency tracker(s) determine a partial derivative with respect to each estimated complex sinusoidal component of a function based on a weighted log-likelihood function and for each frequency estimate, determine the adaptation value from the respective partial derivative.

50. An apparatus according to claim 35 wherein the adaptive frequency tracker(s) produce a new set of frequency estimates by applying an adaptation value to each frequency estimate in a previous set of frequency estimates, the adaptation value being a function of both the component estimates and the frequency deviation estimates to produce an intermediate set of frequency estimates, and using the frequency deviation estimates and previous frequency deviation estimates to produce an estimate of chirp for each sinusoidal component, and for each sinusoidal component combine the frequency deviation estimate and the estimate of chirp to produce a new frequency estimate.

51. An apparatus according to claim 36 wherein the component generator(s) convert the new set of frequency estimates to new estimated complex sinusoidal components by combining previous estimated complex sinusoidal component estimates with the new frequency estimates.

52. A computer in combination with a computer readable medium compatible with the computer, cooperatively adapted to implement a method according to claim 1.