Technical Field

This invention relates to frequency estimation for a wide class of signals. In particular it relates to a method for estimating parameters, including the frequency of a single frequency complex exponential signal in additive Gaussian noise or in additive non Gaussian noise and a method for estimating parameters of a signal composed of a plurality of complex exponentials. This invention further relates to a method for estimating the frequency of a bandpass modulated baseband signal, a method for estimating the frequency of a signal which is spectrally even symmetric about its center frequency and a method for estimating the carrier frequency offset for orthogonal frequency division multiplex (OFDM) modulation. The invention also relates to computer hardware programmed to perform the methods.

Background Art Spectral frequency estimation of a complex exponential signal has long been a problem in the area of signal processing.

First consider a set of discrete samples of a single frequency complex exponential signal denoted as r[n], which can be mathematically represented as: . (1) where f _e is the frequency of the complex exponential signal is the phase of the complex exponential signal is the amplitude of the complex exponential signal

T _s is the sampling period and n is the index of the samples and has a range from 0 to N-1. Taking the discrete Fourier transform of r[n] in (1) results in the discrete spectral distribution which can be expressed as: (2) where k is the index of the discrete spectral distribution and has a range from 0 to

N-i .

Taking the magnitude of R \k\ in (2) yields:

where and k ^e = f _e NT _s .

Evaluating |R[k]| in (3) at k = k _e yields:

Evaluating |R[k]| in (4) at k = k _e ±ε where ε is strictly positive and not of an integer value results in the following relationship: |R [k _e + ε]| = |R [k _e - ε]| and | [k _e ± ε]| = 0 if ε is strictly positive integer value.

It can be deduced from the above relationship that the discrete spectral distribution of a complex exponential is even symmetrical about its frequency. In other words, in noiseless conditions, evaluating two discrete spectral components that are - ^ε away from the point k _e in will result in the said evaluated spectral components having equal magnitude as shown in figure 1.

In general, the expression |R[k]| , in (3) is only evaluated at integer values of k's from 0 to N-1. However, there is no guarantee that , which equals to , is of an integer value. In fact, k _e is most likely to lie between two consecutive integer values of k's in most circumstances. In this case, it is not possible to obtain the value of k ^e from searching for the maximum peak of |R[k]| : in (3).

In the case of frequency estimation of a signal composed of multiple complex exponentials, there are two major classes of algorithms that are used. The older of the two is the classical spectral estimation method and the other one is the relatively new parametric frequency estimation method. The classical spectral estimation method has the advantage of being very computationally efficient but the accuracy of the frequency estimates generally suffers from inter-carrier interference (ICI) caused by spectral leakage from the individual complex exponentials components.

Parametric methods, such as the Multiple Signal Classification (MUSIC) and Estimation of Signal Parameter via Rotational Invariance Techniques (ESPRIT), have been the dominant methods for the simultaneous estimation of frequencies of multiple complex exponentials, in recent years. These algorithms are highly accurate and have the ability to mitigate the problems associated with the spectral estimation method. However, the parametric methods suffer from high computational complexity. Hence they are not suitable for real-time applications. In addition they are not very accurate in low signal to noise (SNR) ratio environments.

One class of techniques for estimating the frequency of either unmodulated signals or modulated bandpass signals, rely on calculating the average of the phase difference induced by the time of arrival difference of samples at some carrier frequency.

For complex sinusoidal signals, one technique determines the phase rotation as a function of time, for successive time samples, since the phase rotation is a function of frequency. For instance Fitz [15] uses the phase angle of an un-normalized sampled autocorrelation function of a received signal to obtain a frequency estimate. For OFDM modulation, Schmidl and Cox [16] adopted a technique which involves L complex samples in the first half of a first training symbol. These complex samples are correlated with the complex samples in the second half of the first training symbol by,

where r _k is the received complex sample at time kT _s .

Because the same symbol is transmitted, the frequency error creates a phase difference between the second half of the training symbol and the first half of the training symbol.

Whilst such techniques are adequate, some algorithms are prone to timing errors. Others, whist overcoming the problem of the introduction of timing errors, do so at the expense of additional computational complexity.

Summary of the Invention

In a first aspect the invention is a method for estimating the frequency and phase of a single frequency complex exponential tone in additive Gaussian noise or in additive non Gaussian noise, comprising: performing a fast Fourier transform (FFT) on the tone; estimating the frequency as the frequency corresponding to the largest FFT output coefficient magnitude; computing a discriminant which is proportional to the frequency error in the initial frequency estimate using modified coefficients of the discrete Fourier transform (DFT) with center frequencies plus one half and minus one half of the FFT bin spacing relative to the initial frequency estimate; mapping the value of the discriminant into the estimate of the frequency error in the initial frequency estimate using a mathematically derived function; adding the estimate of the frequency error to the initial frequency estimate to get a first interpolated frequency estimate; computing a further discriminant which is proportional to the frequency error in the first interpolated frequency estimate using modified coefficients of the discrete Fourier transform (DFT) with center frequencies plus one half and minus one half of the FFT bin spacing relative to the first interpolated frequency estimate;

mapping the value of the further discriminant into the estimate of the "frequency error in the first interpolated frequency estimate using the mathematically derived function; adding the estimate of the frequency error in the first interpolated frequency estimate to the first interpolated frequency estimate to get a second interpolated frequency estimate; and computing a phase estimate from the frequency estimate.

The estimate of k _e , denoted as , can be estimated by exploiting the said even symmetry of the discrete spectral distribution about k _e with contraction mapping. The frequency estimation of the co lex exponential can then be estimated with the relationship , where is the estimated frequency of the said complex exponential tive technique of estimating the frequency of a complex exponential that exploits the spectral symmetry i isclosed.

The first interpolated frequency estimate, , may be in a region of relatively low noise induced frequency error and may therefore be quite accurate.

The technique of iteratively deriving an interpolated frequency estimate and then using the frequency discriminant may be continued for a fixed number of iterations or until the value of the discriminant is below a pre-determined threshold value. The technique may be continued until a fixed point occurs, wherein at this fixed point, the discriminant function has zero value.

After recursively obtaining the frequency estimate, , the phase estimate may be obtained from the relation:

The method may further comprise computing an amplitude estimate from the computed phase and frequency estimates. After recursively obtaining the frequency estimate, and obtaining the phase estimate, the amplitude estimate may be obtained from the relation:

The method may be implemented in computer hardware and/or computer software. This includes, but is not limited to, integrated circuits, field programmable gate array devices, microprocessors and digital signal processors. The method is ideally suited to be utilised in a number of communications, signal processing and biomedical applications.

Several functions have been identified to compute the discriminant. In practice, different functions may require a different number of iterations to essentially converge to a fixed-point solution. However, discriminant functions defined by a wide class of functions using two DFT coefficients as the input converge to the same solution and therefore exhibit identical noise performance.

A first example of the discriminant, or distance metric, of frequency estimation error is: (5) where, (6) and, So for the initial frequency estimate using the FFT, and

Other examples of the discriminant having the properties required for the algorithm include:

where Re[.] is the real part and * denotes the complex conjugate.

In these equations, β and α are the modified DFT coefficients defined by,

It is also possible to define discriminant functions which use more than two DFT coefficients to obtain further improvements in frequency estimation performance in additive Gaussian noise relative to discriminants that use only two DFT coefficients. An example where 2M+2 coefficients are used, where and the FFT coefficients are used in the discriminant with optimal weighting coefficients obtained by using the concept of matched filtering is,

where, , mod N indicates modulo N, and * denotes complex conjugate.

_{J .} and, are the modified DFT coefficients given by,

The discriminant using more than two DFT coefficients may be used in the last iteration to obtain additional frequency accuracy. In a similar manner, discriminant functions may be formulated which use more than two DFT coefficients and less or equal to all N FFT coefficients. Additional frequency accuracy may be obtained by computing the frequency discriminant recursively until convergence for the frequency estimate is reached. Convergence for the frequency estimate may be reached after zero to three iterations, depending upon the specific discriminant used and the signal to noise ratio.

In any iteration, the frequency discriminant may be computed using any one of the functional forms:

where Re[.] denotes the real part and denotes the complex conjugate, γ may vary on each iteration. In general, the frequency incremental shift Δf _m (r) is related to the previously defined frequency discriminant, D, by,

The frequency discriminant may be driven to zero input and output values by either modifying the frequency of the DFT coefficients or frequency translating the signal. Signal frequency translation may be achieved by multiplication of the signal by a locally generated complex exponential signal. The advantage of frequency multiplication of the signal is that the algorithm may be implemented with a standard hardware, software, or combination hardware/software FFT. This FFT may be highly optimized for one or a multiplicity of processors operating as a system.

The process for obtaining additional frequency accuracy may be scaled to save multiplies by scaling the frequency estimate during recursion. The process may involve a final step of multiplying the scaled frequency estimate f _m+1 τ _s with the sampling frequency f _s to remove the scaling from the frequency estimate.

In a second aspect, the invention is a method for estimating the frequency of a bandpass modulated baseband signal comprising: multiplying the complex bandpass frequency modulated signal with its complex conjugate baseband self to produce a signal z[n] that is spectrally symmetric about its center frequency; and estimating the center frequency of z[n] using a frequency estimation algorithm which relies on spectral symmetry.

In one embodiment the frequency estimation algorithm comprises: performing a fast Fourier transform (FFT) on the symmetric signal; estimating the frequency as the frequency corresponding to the largest FFT output coefficient magnitude; computing a discriminant which is proportional to the frequency error in the initial frequency estimate using modified coefficients of the discrete Fourier transform (DFT) with center frequencies plus epsilon and minus epsilon of the FFT bin spacing relative to the initial frequency estimate; mapping the value of the discriminant into the estimate of the frequency error in the initial frequency estimate using a mathematically derived function; adding the estimate of the frequency error to the initial frequency estimate to get a first interpolated frequency estimate; computing a further discriminant which is proportional to the frequency error in the first interpolated frequency estimate using modified coefficients of the discrete

Fourier transform (DFT) with center frequencies plus epsilon and minus epsilon of the FFT bin spacing relative to the first interpolated frequency estimate; mapping the value of the further discriminant into the estimate of the frequency error in the first interpolated frequency estimate using the mathematically derived function; and adding the estimate of the frequency error in the first interpolated frequency estimate to the first interpolated frequency estimate to get a second interpolated frequency estimate.

At each step of the iteration, the modified DFT coefficients may be calculated at

either frequencies, where, is the ith recursive frequency estimate, and, some situations, it may be

advantageous to selected a value of such that, .

The value of epsilon depends on the type of bandpass signal the method is applied to.

The spectral observation of the spectrally symmetric signal class, such as signals which are bandpass in nature, will yield a zero discriminant value under a noiseless condition when the two spectral components are either, plus and minus half a DFT bin for the case of complex exponential, or plus and minus "epsilon" DFT bin/s (where epsilon need not be an integer value) for the general case. In general when the signal has an even symmetric power spectral density (PSD) about its center frequency and the discrete spectrum is observed at "epsilon" bin/s above and below the center frequency the discriminant described will have zero value in the noiseless case. In the presence of AWGN, the discriminant value will be non-zeros, in general. Using recursion in an iterative algorithm, the estimate of the center frequency is successively adjusted in each iteration, which will drive the discriminant to a zero value. It is shown that the steady state frequency estimate is close to the center frequency of the signal. Due to the contraction mapping theorem, the iterative algorithm will always converge to a final value of the frequency estimate. Embodiments of the method in accordance with the second aspect may incorporate a discriminant, or distance metric, as described in embodiments relating to the first aspect of the invention.

In other embodiments, traditional frequency estimation algorithms which rely on spectral symmetry are used. Examples of frequency estimation algorithms which are incorporated herein by reference include [1], [13] and [14].

The method in accordance with the second aspect has the advantage of yielding results close to the Cramer Rao lower bound on frequency estimation error in the additive white Gaussian noise (AWGN) and may be implemented with very lower computational complexity. The performance and computational complexity enable real time digital signal processing implementation of a large number of applications requiring center frequency determination.

The method in accordance with the second aspect may further include computing a phase estimate from the frequency estimate and may further include computing an amplitude estimate from the computed phase and frequency estimates. In a third aspect the invention is a method for estimating the frequency of a bandpass signal having a baseband component, the method comprising: performing an autocorrelation on the bandpass signal; performing an autocorrelation on the baseband component of the bandpass signal; multiplying the auto correlated bandpass signal with the complex conjugate of the auto correlated baseband component of the signal to produce a signal z[n] that is spectrally symmetric about its center frequency; and estimating the center frequency of z[n] using a frequency estimation algorithm which relies on spectral symmetry. The frequency estimation algorithm may be any one of those described in relation to embodiments of the second aspect of the invention.

The method in accordance with the third aspect has the advantage that the averaging performed by the autocorrelation function minimises the effect of that the fading channel has on the performance of the frequency estimator. In addition it reduces the adverse effect that sample and frame timing offset has on the accuracy of the frequency estimate.

In a fourth aspect, the invention is a method for estimating the frequency of a signal which is spectrally even symmetric about its center frequency comprising: performing a fast Fourier transform (FFT) on the signal with a number of sample points N =N ₀ , where N ₀ is a reduced number of sample points; estimating the frequency as the frequency corresponding to the largest FFT output coefficient magnitude; setting Ni = N ₀ ; computing a discriminant which is proportional to the frequency error in the initial frequency estimate using modified coefficients of the discrete Fourier transform (DFT) with center frequencies plus epsilon and minus epsilon of the FFT bin spacing relative to the initial frequency estimate;

mapping the value of the discriminant into the estimate of the frequency error in the initial frequency estimate using a mathematically derived function; adding the estimate of the frequency error to the initial frequency estimate to get a first interpolated frequency estimate; computing a further discriminant which is proportional to the frequency error in the first interpolated frequency estimate using modified coefficients of the discrete Fourier transform (DFT) with center frequencies plus epsilon and minus epsilon of the FFT bin spacing relative to the first interpolated frequency estimate; mapping the value of the further discriminant into the estimate of the frequency error in the first interpolated frequency estimate using the mathematically derived function; and adding the estimate of the frequency error in the first interpolated frequency estimate to the first interpolated frequency estimate to get a second interpolated frequency estimate. setting N,=N ₀₊ i and repeating the steps of iteratively deriving a first and second interpolated frequency estimate.

No degradation to the error variance performance will result when operating at SNR higher or equal to η if the following condition is met:

where η is lowest expected signal to noise ratio of the operating environment.

If the above condition is not met but it is necessary to keep N ₀ as low as possible and N relative large for the require accuracy, intermediate set of iterations can be inserted to bridge between the first and last set of iterations. Denoting N ₁ as the number of samples used in the intermediate set, the frequency estimate can be found in three sets of iterations by the said frequency estimator, with the first set of iterations operated on N ₀ samples points and the second set of iterations operated on N ₁ sample points and the last set of iterations operated on N sample points, without any degradation to the error variance performance when operating at SNR higher or equal to η if the following condition is met:

Notice the number of intermediate sets of iteration can be increased to further reduce the number of samples, No, required for computing the FFT for the initial

frequency search. In general, at any particular SNR, is it possible to determine a sufficiently large N ₀ , N ₁ , N ₂ such that the described procedure may be implemented with negligible degradation in error variance performance compare to the algorithm implemented with N samples. The steps of iteratively deriving an interpolated frequency estimate and using the frequency discriminant may be repeated for each N, until either the number of iterations has reached a pre-defined number or the output of the discriminant is below a predefined threshold, or until a fixed point solution occurs, where at this fixed point the discriminant function has zero value. The signal may be a bandpass modulated signal and the method may further include the initial step of multiplying the bandpass modulated signal with the complex conjugate of the baseband signal to produce a signal that is symmetric about its center frequency.

The value of epsilon depends on the type signal to which the method is applied. The signal may be a complex exponential signal.

The method in accordance with the fourth aspect may further include computing a phase estimate from the frequency estimate and may further include computing an amplitude estimate from the computed phase and frequency estimates.

The method of the second, third or fourth aspect, or any of its embodiments may further include filtering, linear filtering, median filtering, or any combination of time filtering and order statistical filtering of a time sequence of estimates to obtained improved estimation accuracy.

The method of the second, third or fourth aspect may be used in a feedback control loop to vary the frequency of a local oscillator (either voltage controlled oscillator, numerically controlled oscillator, or digital output numerically controlled oscillator), in such a manner as to align the FFT output frequencies to the frequencies of the received carrier signals, such as in the case of OFDM (orthogonal frequency division multiplexed) demodulation.

The method of the second, third or fourth aspect may be used in a feedback control loop to vary the frequency of a local oscillator (either voltage controlled oscillator, numerically controlled oscillator, or digital output numerically controlled oscillator), in such a manner as to align the frequency down converted signal center frequency to the center frequency of a demodulator.

The method of the second, third or fourth aspect may be used to track the carrier frequency of an ultra wideband signal in a feedback control loop to vary the frequency of a local oscillator (either voltage controlled oscillator, numerically controlled oscillator, or digital output numerically controlled oscillator, in such a manner as to

align the signal at the output of the complex mixer to the demodulator center frequency (whether the center frequency is zero or non-zero).

The method of the second, third or fourth aspect may be incorporated into a digital demodulator which involves functions of carrier frequency tracking, channel estimation, and/or digital demodulation. The method of the second or third aspect may be incorporated into a pre-processor or a co-processor for pattern recognition.

Broader applications of the method of the second third or fourth aspect, or any of its embodiments, include the estimation of one or more of frequency amplitude and phase, of radar signals, sonar signals, ultra sound signals, tracking systems and tracking loops for demodulation for time division multiple access communication systems and burst mode digital communication systems, audio signals with a low pass reference signal, biomedical signals including, electroencephalogram (EEG) ₅ electrocardiogram (EKG), ultrasound imaging, ultrasound Doppler measurements, and other biomedical signal types, and bandpass signals including the Doppler shift of a signal. An advantage of embodiments of the second, third or fourth aspects of the invention is that the method is able to be applied to each of multiple antenna outputs to estimate the center frequency, amplitude, and phase of the signal. The estimates can then be used to estimate the angle of arrival of the signal to the array of antennas.

An advantage of embodiments of the second, third or fourth aspects of the invention is that the method is able to be applied to each of multiple antenna outputs to estimate the center frequency, amplitude, and phase of the signal. The angle of arrival of the signal may be obtained from these estimates and the antenna gain pattern can be electronically directed to the direction of the angle of arrival in either a closed loop or an open loop control configuration. An advantage of embodiments of the second, third or fourth aspects of the invention is that the method is able to be applied to each of multiple antenna outputs to estimate the center frequency, amplitude, and phase of the multiple carrier signals transmitted over a multiple input multiple output (MIMO channel).

In a fifth aspect, the invention is a frequency and phase estimation program for estimating the frequency and phase of a single frequency complex exponential tone in additive Gaussian noise or additive non-Gaussian noise, wherein the estimation program has functionality to perform the method in accordance with the first aspect, or an embodiment of the first aspect.

In a sixth aspect, the invention is computer hardware programmed to perform the method in accordance with the first, second, third or fourth aspects, or an embodiment of the any of these aspects. The hardware may comprise a DSP processor chip, or any other programmed hardware.

In a seventh aspect the invention is a method for estimating the carrier frequency offset for orthogonal frequency division multiplex (OFDM) modulation comprising: performing coarse frequency estimation on a received OFDM symbol which is subject to a fading channel and a carrier frequency offset; compensating the received OFDM symbol with the coarse carrier frequency offset estimate; estimating a channel characteristic from the compensated symbol; equalizing the compensated symbol with the channel characteristic; detecting the modulation of data carriers from the equalized symbol; forming a function from detected data modulation and the channel characteristic; multiplying the compensated symbol with the complex conjugate of the function to form a signal z[n] that is spectrally symmetric about its center frequency; and estimating the carrier frequency offset estimate of z[n].

The step of estimating the carrier frequency offset estimate of the signal z[n] may utilise a frequency estimation algorithm as described in relation to any one of the embodiments of the second aspect of the invention.

The received OFDM symbol which is subject to a fading channel and a carrier frequency offset may be represented mathematically as:

where d _k is the modulation of the k ^th data carrier of the transmitted OFDM signal, h _k is the k ^x DFT coefficient of the channel frequency response (CFR)

Δ/ is the frequency of the CFO

Aθ is the phase of the CFO f _k is the k ^th frequency corresponding to the k ^th DFT bin; and w _h [n] is the sequence of noise at the output of the channel.

The step of performing coarse frequency estimation on a received OFDM symbol may be performed with any integer CFO estimation algorithm such as that disclosed in references [11], [12], which are hereby incorporated by reference. The method may comprise estimating the channel characteristic in the form of one of a channel frequency response (CFR) and a channel impulse response (CIR) using a channel estimation technique.

The method may further comprise forming the data symbols by detecting the modulation of data carriers from the equalized compensated symbol. In one embodiment the function r _est [n] is formed from one of the relations, and

where represents the modulation of the data carriers, is the known or estimated coefficients of the CFR at the k ^th DFT bin, and h[n] is the set of known or estimated coefficients of the CIR.

The method may comprise iteratively deriving the channel characteristic, detected data modulation and carrier frequency offset estimate and updating the channel characteristic, detected data modulation and carrier frequency offset with each iteration for a predetermined number of times. The total number of iterations can be adjusted according to bit error rate requirement at the demodulator. Lower bit error rate would imply less error on the frequency estimate hence more iterations and visa versa.

The method may comprise iteratively deriving the channel characteristic, detected data modulation and carrier frequency offset estimate and updating the channel characteristic, detected data modulation and carrier frequency offset with each iteration until the output of the discriminant of the frequency estimation is below a predetermined threshold value.

In an eighth aspect the invention is a method for estimating the carrier frequency offset for orthogonal frequency division multiplex (OFDM) modulation when both the channel characteristic and transmitted data symbols are known, comprising: forming a function from the data symbols and the channel characteristic; multiplying a received OFDM symbol which is subject to a fading channel and carrier frequency offset with the complex conjugate of the computed function to form a signal z[n] that is symmetric about its center frequency and; estimating the carrier frequency offset estimate of the signal z[n]. The step of estimating the carrier frequency offset estimate of the signal z[n] may utilise a frequency estimation algorithm as described in relation to any one of the embodiments of the second aspect of the invention.

The computed function r _est [n] may be formed from one of the relations, and

where h _k is the known DFT coefficients of the CFR at the k ^th DFT bin, and h[n] is the set of known coefficients of the CIR.

The function formed from the data symbols and the channel characteristic may be represented mathematically as,

In a ninth aspect the invention is a method for estimating the carrier frequency offset for orthogonal frequency division multiplex (OFDM) modulation in additive white Gaussian noise, comprising the steps of: performing coarse frequency estimation on a received OFDM symbol which is subject to a fading channel and a carrier frequency offset; compensating the received OFDM symbol with the coarse carrier frequency offset estimate; detecting the modulation of data carriers from the compensated symbol; forming a function from detected data modulation; multiplying the compensated symbol with the complex conjugate of the formed function to form z[n] and estimating the carrier frequency offset estimate of z[n].

The step of estimating the carrier frequency offset estimate of the signal z[n] may utilise a frequency estimation algorithm as described in relation to any one of the embodiments of the second aspect of the invention.

The received OFDM symbol which is subject to a carrier frequency offset may be represented mathematically as,

The method may comprise iteratively deriving the detected data modulation and carrier frequency offset estimate and updating the detected data modulation and carrier frequency offset with each iteration for a predetermined number of times. Optionally, the method may comprise iteratively deriving the detected data modulation and carrier frequency offset estimate and updating the detected data modulation and carrier frequency offset until the output of a discriminant of the frequency estimation is below a pre-determined threshold value.

The step of performing coarse frequency estimation on a received OFDM symbol may be performed with any integer CFO estimation algorithm such as that disclosed in references [11], [12], which are hereby incorporated by reference.

In a tenth aspect the invention is a method for estimating the carrier frequency offset for orthogonal frequency division multiplex (OFDM) modulation in additive white Gaussian noise and the transmitted data symbols are known, comprising: forming a function r _est [n] from known data, where

multiplying the received symbol with the complex conjugate of r _est [n] and estimating the carrier frequency offset estimate.

The step of estimating the carrier frequency offset estimate of z[n] may comprise applying a method in accordance with any one of the previous aspects, or with any one of their embodiments. In another embodiment traditional frequency estimation techniques may be used. In an eleventh aspect the invention is a method for estimating the carrier frequency offset for an ultrawide band signal assuming the transmitted signal is known, comprising: multiplying the signal with the complex conjugate of the known transmitted signal and estimating the carrier frequency offset estimate using any of the above mentioned frequency estimation techniques or an existing frequency estimation technique.

In a twelfth aspect, the invention is computer hardware programmed to perform the method in accordance with the sixth, seventh, eight, ninth, or tenth aspects, or an example of the any of these aspects. The hardware may comprise a DSP processor chip, or any other programmed hardware.

In a thirteenth aspect the invention is a method for estimating the frequency, amplitude and phase of a signal composed of a plurality of complex exponentials, the method comprising: performing a fast Fourier Transform on the signal; estimating the frequency as the frequency corresponding to the largest FFT output coefficient magnitude; refining the initial frequency estimate by: computing a discriminant which is proportional to the frequency error in the initial frequency estimate using modified coefficients of the discrete Fourier transform (DFT) with center frequencies plus epsilon and minus epsilon of the FFT bin spacing relative to the initial frequency estimate; mapping the value of the discriminant into the estimate of the frequency error in the initial frequency estimate using a mathematically derived function; adding the estimate of the frequency error to the initial frequency estimate to get a first interpolated frequency estimate; computing a further discriminant which is proportional to the frequency error in the first interpolated frequency estimate using modified coefficients of the discrete Fourier transform (DFT) with center frequencies plus epsilon and minus epsilon of the FFT bin spacing relative to the first interpolated frequency estimate; mapping the value of the further discriminant into the estimate of the frequency error in the first interpolated frequency estimate using the mathematically derived function; and

adding the estimate of the frequency error in the first interpolated frequency estimate to the first interpolated frequency estimate to get a second interpolated frequency estimate; using the refined frequency estimate to estimate the phase and amplitude; simplifying the signal by subtracting from the signal the complex exponential defined by the refined estimate of frequency and estimates of the amplitude and phase; and repeating the steps by substituting the simplified signal as the signal until a refined estimate of frequency and estimates of the amplitude and phase are obtained for each complex exponential.

In accordance with the twelfth aspect, the step of simplifying the signal by subtracting from the signal the complex exponential defined by the refined estimate of frequency, amplitude and phase may be performed either in the time domain or the frequency domain. The method in accordance with the thirteenth aspect may include further refining of the frequency, amplitude and phase estimates for each of the complex exponentials. This may be obtained via an iterative process which continues until the estimates for each parameter for each of the respective complex exponentials substantially converge. Further refining of the frequency, amplitude and phase estimates for each of the complex exponentials may be referred to as a method of fine estimation of each of the parameters. In a preferred embodiment the method of fine estimation may include, for each complex exponential, subtracting from the original signal all estimated complex exponentials formed by their current frequency, phase and amplitude estimates, except the complex exponential of interest. The frequency of the complex exponential of interest may then be refined by: computing a discriminant which is proportional to the frequency error in the frequency estimate using modified coefficients of the discrete Fourier transform (DFT) with center frequencies plus epsilon and minus epsilon of the FFT bin spacing relative to the frequency estimate; mapping the value of the discriminant into the estimate of the frequency error in the frequency estimate using a mathematically derived function; adding the estimate of the frequency error to the frequency estimate to get a first interpolated frequency estimate; computing a further discriminant which is proportional to the frequency error in the first interpolated frequency estimate using modified coefficients of the discrete Fourier transform (DFT) with center frequencies plus epsilon and minus epsilon of the FFT bin spacing relative to the first interpolated frequency estimate;

mapping the value of the further discriminant into the estimate of the frequency error in the first interpolated frequency estimate using the mathematically derived function; and adding the estimate of the frequency error in the first interpolated frequency estimate to the first interpolated frequency estimate to get a second interpolated frequency estimate.

The fine estimates of frequency may then be used to obtain fine estimates of the phase and amplitude of the complex exponential of interest.

The frequency, phase and amplitude of the other complex exponentials can be found by repeating these steps.

An advantage of at least an embodiment of the twelfth aspect of the invention is that the problem of estimating the frequency, phase and amplitude of multiple complex exponentials is reduced to estimating the said parameters of a single complex exponential. The method in accordance with the thirteenth aspect or an embodiment of the twelfth aspect may be implemented in a sequential manner where the parameters of M individual complex exponentials are estimated one at a time. Alternatively the method may be implemented in a parallel manner. For instance parameters of the M individual complex exponentials can be produced as a separate process, either each parameter implemented on a separate processor or multi-threaded on several processors if M processors are not available. Optionally, or in addition, the further refinement of parameters for each of the complex exponentials may occur on a separate processor or multi-threaded on several processors. Parallelisation has the advantage of reducing the computational time required to obtain the estimates. The method in accordance with the thirteenth aspect or an embodiment of the twelfth aspect may further comprise a preliminary step of determining the number of complex exponentials present in the signal. Such a step may comprise determining whether local maxima in the signal spectrum are frequency components or whether the local maxima in the spectrum are fluctuations due to spectrum noise. In an embodiment such a step may comprise determining the magnitude of each complex exponentials FFT spectral coefficients and comparing the respective magnitudes to an estimated value of the magnitude of the FFT coefficient for frequency bins that do not contain signal. In such an example if any particular FFT coefficient magnitude is sufficiently larger than an estimated noise only coefficient magnitude, then it may be determined that the frequency bin associated with the particular coefficient magnitude contains a complex exponential. However if any particular FFT coefficient magnitude is not sufficiently larger than the estimated noise only coefficient magnitude, then it may be determined

that the frequency bin associated with the particular FFT coefficient magnitude does not contain a complex exponential.

In a fourteenth aspect the invention is a parameter estimation program for estimating the frequency, amplitude and phase of a signal composed of a plurality of complex exponentials, wherein the program has functionality to perform the method in accordance with the twelfth aspect.

In a fifteenth aspect the invention is computer hardware programmed to perform the method in accordance with the twelfth aspect. The hardware may comprise a digital signal processor, a field programmable gate array device, an integrated circuit, a microprocessor, a computer, an embedded processor, or other signal processing devices using electronic hardware, software, or a combination of hardware and software.

Embodiments of the thirteenth, fourteenth and fifteenth aspects of the invention in which the frequency of each of the complex exponentials are reasonably far apart from each other, demonstrate superior performance in terms of complexity and accuracy than the previous known state of the art parameter estimator such as MUSIC and

ESPRIT.

A further advantage of at least an embodiment of the thirteenth, fourteenth and fifteenth aspects of the invention is that the likelihood of spectral leakage or spectral smearing affecting the accuracy of the estimation of one or more parameters is reduced. An advantage of an implementation of at least one embodiment of the thirteenth, fourteenth and fifteenth aspects of the invention is that the implementation has low computational complexity and results in extremely accurate estimates of the multiple frequencies, phases, and amplitudes.

Brief Description of the Drawings

Examples of the invention will now be described with reference to the accompanying drawings, in which:

Fig. 1 is a graph illustrating the FFT Coefficients, where the discrete spectral distribution is even symmetrical about its center frequency; Fig. 2 is a graph illustrating the FFT Coefficients, where the signal frequency is closer to the lower FFT frequency than the higher FFT frequency;

Fig. 3 is a graph illustrating the FFT Coefficients, there are two equal peak coefficients and the signal frequency is halfway between;

Fig. 4 is a graph illustrating the FFT Coefficients, where the signal frequency is closer to the upper FFT frequency than the lower FFT frequency; Fig. 5 is a flow diagram for an embodiment of the invention;

Fig. 6 is a graph that illustrates the ratio of the variance of the normalized frequency error, , to Cramer-Rao Bound variance in dB as a function of the FFT length, N; and

Fig. 7 is a graph that illustrates the variance of the normalised estimator frequency error estimate against the frequency error for the first interpolation.

Simulations of the invention show the rms frequency error performance of the algorithm vs SNR in dB, for different values of N. Figures 8-13 include curves for a single interpolation, two interpolations, and the Cramer-Rao Bound, where:

Fig. 8 is a graph showing RMS normalised frequency error vs SNR for N=2; Fig. 9 is a graph showing RMS normalised frequency error vs SNR for N=4;

Fig. 10 is a graph showing RMS normalised frequency error vs SNR for N=I 6;

Fig. 11 is a graph showing RMS normalised frequency error vs SNR for N=64;

Fig. 12 is a graph showing RMS normalised frequency error vs SNR for N=256;

Fig. 13 is a graph showing RMS normalised frequency error vs SNR for N=I 024;

Fig. 14 is a flow diagram for an embodiment of the invention using a fixed number of iterations stopping rule;

Fig. 15 is a flow diagram for an embodiment of the invention using a magnitude of the frequency error discriminant stopping rule; Fig. 16 is a graph that illustrates the comparison of the frequency estimation error variance versus SNR in dB between the modified SSI frequency estimation algorithm and the frequency estimation algorithm in accordance with the first aspect of the invention;

Fig. 17 is a graph illustrating a comparison of the phase estimate error variance verus SNR in dB between the described phase estimator using the frequency estimate from the modified SSI frequency estimation algorithm and the phase estimator using the frequency estimate from the first aspect of the invention;

Fig. 18 is a graph illustrating a comparison of the amplitude estimate error variance verus SNR in dB between the described amplitude estimator using the frequency estimate from the modified SSI frequency estimation algorithm and the said amplitude estimator using the frequency estimate from the frequency estimation algorithm the first aspect of the invention;

Fig. 19 is a graph illustrating a comparison of the CFO estimate error variance verus SNR in dB between the described "OFDM algorithm IV" employing the modified SSI frequency estimation algorithm and the "OFDM algorithm IV" employing the SSI frequency estimation algorithm;

Fig. 20 is a graph illustrating the frequency estimate error variance verus SNR in dB of the described "OFDM algorithm 1" employing the SSI frequency estimation algorithm;

Fig. 21 to Fig. 26 illustrate an example of the effect of frequency and phase offset has on UWB signals and the effectiveness of the SSI frequency estimation algorithm and the phase estimator of estimating the said frequency and phase offset:

Fig. 21 illustrates the real component of the UWB pulse in the noiseless case;

Fig. 22 illustrates the imaginary component of the UWB pulse in the noiseless case; Fig. 23 illustrates the real component of the same UWB pulse subjected to a frequency and phase shift

Fig. 24 illustrated the imaginary component of the same UWB pulse subjected to a frequency and phase shift;

Fig. 25 illustrates the real component of the recovered UWB pulse by compensating the frequency and phase shifted version of the pulse with the frequency and phase estimate from the said SSI frequency estimation algorithm and phase estimator;

Fig. 26 illustrates the imaginary component of the recovered UWB pulse by compensating the frequency and phase shifted version of the pulse with the frequency and phase estimate from the said SSI frequency estimation algorithm and phase estimator;

Fig. 27 to Fig 31 illustrate an example of estimating and compensating for the frequency and phase offset of a noise corrupted UWB pulse;

Fig. 27 illustrates the real component of the UWB pulse at SNR equal to 5 dB; Fig. 28 illustrates the imaginary component of the UWB pulse at SNR equal to 5 dB;

Fig. 29 illustrates the real component of the noise corrupted UWB pulse subjected to the same frequency and phase shift as in the example above;

Fig. 30 illustrates the imaginary component of the noise corrupted UWB pulse subjected to the same frequency and phase shift as in the example above;

Fig. 31 illustrates the real component of the recovered noise corrupted UWB pulse by compensating the frequency and phase shifted version of the pulse with the frequency and phase estimate from the said SSI frequency estimation algorithm and phase estimator; Fig. 32 illustrates the imaginary component of the recovered noise corrupted

UWB pulse by compensating the frequency and phase shifted version of the pulse with the frequency and phase estimate from the said SSI frequency estimation algorithm and phase estimator;

Fig. 33 to Fig 38 illustrate another example of estimating and compensating for the frequency and phase offset of a noise corrupted UWB pulse as described above except the SNR is not equal to 0 dB;

Fig 39 shows the frequency estimate error variance verus SNR in dB of the SSI frequency estimation algorithm; and

Fig 40 shows the phase estimate error variance verus SNR in dB of the phase estimator employing the frequency estimate from SSI frequency estimation algorithm.

Each of figures 41 to 52 illustrate a line that represents the Cramer-Rao lower bound (CRLB) for frequency estimation of a signal composed of a plurality complex exponentials and a simulated result of the error variance performance;

Figure 41 is a graph that illustrates the frequency error variance performance against SNR in dB for a signal having a frequency spacing of 3/(2NT _s );

Figure 42 is a graph that illustrates the phase error variance performance against SNR in dB for a signal having a frequency spacing of 3/(2NT _s ); Figure 43 is a graph that illustrates the amplitude error variance performance against SNR in dB for a signal having a frequency spacing of 3/(2NT _s );

Figure 44 is a graph that illustrates the frequency error variance performance against SNR in dB for a signal having a frequency spacing of 2/(NT _s );

Figure 45 is a graph that illustrates the phase error variance performance against SNR in dB for a signal having a frequency spacing of 2/(NT _s );

Figure 46 is a graph that illustrates the amplitude error variance performance against SNR in dB for a signal having a frequency spacing of 2/(NT _s );

Figure 47 is a graph that illustrates the frequency error variance performance against SNR in dB for a signal having a frequency spacing of 3/(NT _s ); Figure 48 is a graph that illustrates the phase error variance performance against

SNR in dB for a signal having a frequency spacing of 3/(NT _s );

Figure 49 is a graph that illustrates the amplitude error variance performance against SNR in dB for a signal having a frequency spacing of 3/(NT _s );

Figure 50 is a graph that illustrates the frequency error variance performance against SNR in dB for a signal having a frequency spacing of 4/(NT _s );

Figure 51 is a graph that illustrates the phase error variance performance against SNR in dB for a signal having a frequency spacing of 4/(NT _s );

Figure 52 is a graph that illustrates the amplitude error variance performance against SNR in dB for a signal having a frequency spacing of 4/(NT _s ); Figure 53 is a graph that illustrates a comparison performance of the frequency error variance performance between the multiple complex exponential parameter estimation algorithm, the root MUSIC algorithm and the ESPRIT-TLS algorithm; and

Figure 54 is a schematic flow diagram of an alternative example of the invention which incorporates a detector to detect the number of complex exponentials present in a signal.

Best Modes of the Invention

Referring first to Figs. 2 to 5, the received signal r[n] is given by: (7) where:

is a set of independent, complex, zero mean, Gaussian random variables η _I [n] = Imag{η[n]},

f is the frequency of the tone, T _s is the sampling period,

A is the signal amplitude and the phase, θ _e = 0.

A fast Fourier transform is performed and the sampling frequency, f _s , is given by, (8)

Then, an initial frequency estimate is taken as the frequency corresponding to the largest FFT output coefficient magnitude. A discriminant which is proportional to the frequency error in the initial frequency estimate is computed using modified coefficients α ₀ , β ₀ of the discrete Fourier transfor FT) with center frequencies plus one half and minus one half of the FFT bin spacing relative to the initial frequency estimate . The value of the discriminant is then mapped into the estimate of the frequency error in the initial frequency estimate using a mathematically derived function. The estimate of the frequency error i ded to the initial frequency estimate to get a next interpolated frequency estimate The process is then repeated, using the next interpolated frequency estimate and computing a new frequency discriminant to produce a next, more precise, frequency estimate f ₂ .

THE FREQUENCY INTERPOLATION DISCRIMINANT It is assumed that the signal to noise ratio (SNR) is sufficiently high such that the largest magnitude FFT coefficient corresponds to a frequency closest to the signal frequency. This assumes that the signal to noise ratio is sufficiently high that the

probability of the statistical outlier event of a noise only FFT bin magnitude being larger than a FFT bin containing both signal and noise is negligible. Define,

Then the discriminant, or distance metric, of frequency estimation error is defined as,

and,

For the initial frequency estimate using the FFT,

In the noiseless case,

O(ε, ε) is a monotonically increasing function of ε - ε . Therefore, each there is a unique inverse mapping to Clearly, may be used as a discriminant for fine frequency interpolation between FFT bin center frequencies. There exists some functional relationship such that,

where, ψ( .) is a monotone increasing function. ψ(.) is called the frequency interpolation function and f _x is the first interpolated frequency estimate.

The requirement that f _j has zero error in the noiseless case is,

THE FREQUENCY INTERPOLATION FUNCTION

Without loss of generality, assume that έ = 0. Also assume the noise free case.

The FFT output coefficients are given by,

The discriminant can be expressed as ,

After some trigonometric simplification,

This inverse mapping from D(ε, ε) to ε — ε can be obtained as,

Then the first interpolated frequency estimate, f _x ,may be obtained, where,

The implication is that a two point (N=2) FFT is sufficient to obtain zero error frequency determination in the noiseless case. However, the Cramer-Rao lower bound is relatively large for N=2 and the SNR threshold is relatively large. There is a motivation to use larger N to reduce the rms frequency estimation error. However, the implication of larger N is increased computational complexity and longer delay time to

obtain the transform results. It is desirable to obtain low frequency estimation error with the smallest possible N.

For large N or for any N and small |D|, the function ψ ₁ (D) closely approximates Ψ(D) .

FREQUENCY ESTIMATION BY ITERATION

For the case of r[n] consisting of signal plus noise, the noise will cause a perturbation of D(^, ε), and some error in the frequency estimate will result. Although an algorithm for the exact frequency determination in the noiseless case has been presented, it will be shown that the noise performance of this algorithm improves substantially when | ε — ε | is close to zero. Since the discriminant O(ε,έ) can be used to get an interpolated frequency estimate with less than one half FFT bin size error, it then follows that the algorithm can be iterated to move \ ε -έ\ towards zero and I D(S-, ε) I towards zero. In this way the variance of the frequency estimator output can be reduced.

The iterative algorithm is defined below.

Define,

Define a monotone increasing function of

The iterative algorithm is defined by,

which implies tha . Then,

The steady state frequency estimate at the end of the iteration is . Define Then, and, the normalized frequency error is,

The iteration may be viewed as a convergence to a fixed point of the equation,

[Theorem 1 below is referenced 8];

Let g(x) be a continuous function on and

Furthermore, assume that there is a constant 0 < λ < 1, with, , Then

x = g(χ) has a unique solution x _∞ in [a,b], also, the iteration will converge to x _∞ for any choice of X ₀ ∈ [a,b], and,

In the situation under analysis, for fixed and ε, ^?(.)is a function of

Using Theorem 1, it follows that the iteration will always converge to a fixed point under the appropriate conditions. Also, Φ(D _∞ ) = 0,and from the properties of Φ(.), D _∞ = 0.

The fixed point solution, ^satlsfies >

The two previously defined functions Ψ(D) and ψ ₁ (D) fulfil the requirements of Φ(D) and may be used in the iteration to obtain . While Ψ(D) iteration will tend to converge more rapidly than ψ _x (D) iteration, both will yield identical values of . However, evaluation of ψ _x (D) has lower computational complexity than evaluation of ψ(D) . There is performance advantage in using ^"(D) when the computation is limited to a few iterations.

ADDITIONAL FREQUENCY ACCURACY Assuming that the SNR is sufficiently high, it is highly probable that A fine interpolation is obtained to improve the frequency accuracy.

The normalized frequency estimate is computed recursively in order to save computational complexity. After computation o is obtained by

First Algorithm

Referring to Figures 14 and 15, a first algorithm is provided to improve the accuracy of the frequency estimation. At step 1, the N point complex FFT is computed. The FFT output coeffficents are Y(k), 0≤k<N-1.

At step 2, the peak search to find k _max is: k _max =max ^-1 [|Y[k] |: 0 < k ≤ N -1] . At step 3, the initial frequency estimate is computed by:

At step 4, recursion is started at m=0.

At step 5, the DFT coefficients for the m:th frequency estimate are computed:

The frequency discriminant is then computed:

The m+1 :th frequency estimate is computed:

At step 6, convergence for the frequency estimate is reached if is sufficiently small. If there is convergence, then the frequency estimate is If convergence for the frequency estimate has not been reached, then m is incremented by 1 and step 5 is repeated. In practice, the algorithm will converge for m=2. Therefore, onl for m=0, 1, and 2 need to be computed (two iterations) which means that the frequency estimate is

Second Algorithm

To handle large N, a modification to the first algorithm is made. The modification is made to the step of computing the frequency discriminant:

Third Algorithm

A third algorithm is provided to improve the frequency accuracy. At step 1, the the N point complex FFT is computed. The FFT output coefficients are Y(k), O≤k≤N-1. At step 2, the peak search to find k _max is: k _max =max ^-1 [| Y[k] |: 0 < k ≤ N - 1]

At step 3, the initial frequency estimate is computed by:

At step 4, the DFT coefficients for the initial frequency estimate are computed:

The frequency discriminant is then computed:

The first interpolated frequency estimate is computed:

At step 5, recursion is started at m=l At step 6, the DFT coefficients for the m:th frequency estimate are computed:

The frequency discriminant is then computed:

The m+1 :th frequency estimate, , is computed:

At step 7, convergence for the frequency estimate is reached i is sufficiently small. If convergence has been reached, then the frequency estimate is . If convergence for the frequency estimate has not been reached, then m is incremented by 1 and step 6 is repeated. In practice, the algorithm will converge after f ₂ is evaluated. Therefore, only for m=0, 1, and 2 need to be computed. This algorithm is less computationally complex than the first algorithm and has essentially the same convergence properties in the recursion.

Fourth Algorithm

A fourth algorithm is provided to improve the frequency accuracy. At step 1, the N point complex FFT is computed. The FFT output coefficients are Y(k), 0<k<N-1. At step 2, the peak search to find k _max is: k _max =max ^-1 [| Y[k] |: 0 ≤ k < N -1]

At step 3, the initial frequency estimate is computed by:

At step 4, the DFT coefficients for the initial frequency estimate are computed:

The frequency discriminant is then computed:

The first interpolated frequency estimate is computed:

At step 5, recursion is started at m=l

At step 6, the DFT coefficients for the m:th frequency estimate are computed:

The frequency discriminant is then computed:

The m+l:th frequency estimate, , is computed: , convergence for the frequency estimate is reached if is sufficiently small. If convergence has been reached, then the frequency estimate i

If convergence for the frequency estimate has not been reached, then m is incremented by 1 and step 6 is repeated, hi practice, the algorithm will converge after f ₂ is evaluated.

Therefore, only for m=0, 1, and 2 need to be computed. This algorithm is less computationally complex than the first, second or third algorithms and has essentially the same convergence properties in the recursion. There is reduced computational complexity in the computation of Δf ₁ (r) because of the elimination of the need to compute square roots in the evaluation of the absolute value of complex variables.

Fifth Algorithm

A fifth algorithm is provided to improve the frequency accuracy. At step 1 , the N point complex FFT is computed. The FFT output coefficients are Y(k), 0≤k≤N-1. At step 2, the peak search to find k _max is: k _max =max ^-1 [| Y[k] |: 0 < k < N - 1]

At step 3, the initial frequency estimate is computed by:

At step 4, recursion is started at m=0

At step 5, the DFT coefficients for the m:th frequency estimate are computed:

The frequency discriminant is then computed: , where . γ _m is a set of non-negative constants. γ _m is a cons tan t, γ _m > 0

The m+l:th frequency estimate, is computed:

At step 6, convergence for the frequency estimate is reached if is sufficiently small. If convergence has been reached, then the frequency estimate is If convergence for the frequency estimate has not been reached, then m is incremented by 1 and step 5 is repeated.

Scaled Computational Version of an algorithm

The frequency scaled frequency estimate can be computed and then multiplied by f _s . This process is described using the example of the first algorithm but can similarly be done for all the algorithms. The scaled computational version of the first algorithm is more computationally efficient as it saves multiplies.

At step 1, the N point complex FFT is computed. The FFT output coefficients are Y(k), 0≤k≤N-l. At step 2, the peak search to find k _max is: k _max =max ^-1 [| Y[k] |: 0 ≤ k ≤ N -1]

At step 3, the initial frequency estimate is computed by:

At step 4, recursion is started at m=0

At step 5, the DFT coefficients for the m:th frequency estimate are computed:

The frequency discriminant is then computed:

The m+1 :th frequency estimate, , is computed: At step 6, convergence for the frequency estimate is reached if is sufficiently small. If convergence has been reached then the scaled frequency estimate is f _m+1 T _s . If convergence for the frequency estimate has not been reached, then m is incremented by 1 and step 5 is repeated.

As a final operation, one additional multiply is require to remove the scaling from the frequency estimate.

^smce T _s f _s =1, where is the m+l:th frequency estimate.

Sixth Algorithm

A sixth algorithm is provided to improve the frequency accuracy. The sixth algorithm uses any of the previously defined functional forms for Af _1n (r) for any step.

The difference between the sixth algorithm and the other algorithms types is that the frequencies of the two modified DFT coefficients are not changed. Instead, the center frequency of the signal is modified by multiplying the previously defined signal by

_{tQ obtain} . The effect of this multiplication is to frequency translate the signal by - Δf _m Hertz. The frequency discriminant is driven to zero input and output values by either modifying the frequency of the DFT coefficients or frequency translating the signal. The principle of driving the discriminant to zero by recursion is the same. The frequency error performance of the algorithm as a function of signal to noise ratio is the same whether the signal is frequency translated or whether the DFT coefficients are frequency shifted.

There are situations where standard FFT or DFT functions are available in hardware, software, or combined hardware and software configurations. These FFT and DFT functions are highly optimised for their respective signal processors and are run at very high computational efficiency. Often, parallel processing for multiple processors is utilised extremely effectively. In these cases, the technique of frequency translation of the signal is of considerable implementation benefit. Very efficient computation is achievable. Frequently, an optimised large N point FFT runs faster on a parallel processor than the computation of two DFT coefficient.

For the sixth algorithm, Z ₀ (n)= r(n), n = 0,1,2,3,... is initialised.

At step 1, the N point complex FFT is computed. The FFT output coefficients are Y(k),

O≤k≤N-1.

At step 2, the peak search to find k _max is: k _max =max ^-1 [|Y[k] |: 0 < k < N -1] At step 3, the initial frequency estimate is computed by:

At step 4, recursion is started at m=0

At step 5, the DFT coefficients for the m:th frequency estimate are computed:

The frequency discriminant, Af _m (r), is then computed for any of the functional forms as a function of α _m and β _m .

The m+1 :th frequency estimate, , is computed:

At step 6, convergence for the frequency estimate is reached if is sufficiently small. If convergence has been reached, then the frequency estimate is .

If convergence for the frequency estimate has not been reached, then the signal is frequency translated (by complex time domain multiplication): , for n = 0, 1, 2, ..., N-1. Then m is incremented by 1 and step 5 is repeated.

For convergence, all algorithms have the same steady state performance.

In the noiseless case, the first algorithm gives the exact frequency for one iteration, which is of great benefit in some applications. The exact frequency is obtained for any N, including N=2. The reason is that an exact functional mapping from the magnitudes of the two DFT coefficients to the frequency was analytically derived and used in the algorithm. This is a new analytical result and forms the basis of the algorithm.

PERFORMANCE ANALYSIS For the case of signal in additive noise, D _∞ is a random variable.

In general, |α| and |β| are both Rician distributed random variables. However, under high signal to noise ratio, both |α| and |β| are essentially Gaussian distributed random variables . ε _∞ , which is a random variable, is found by the constraint,

ε _∞ will be perturbed by the noise component in D _∞ . Even though D _∞ is constrained to be zero, the constraint and noise induce randomness in ε _∞ . The noise perturbation in D _∞ induces the perturbation in ε _∞ .

The approach taken is the computation of the variance of D from the point of view of the creation of D from noisy observations and then to find the corresponding perturbation o

For high signal to noise ratios, may be represented by the first three terms of the Taylor Series expansion about, μ _a and μ _β , which are the means of |α| and

|β|, respectively.

Assuming then and for high signal to noise ratio, Then, E[D]=0, it then follows that,

Then, for high signal to noise ratio, the normalized frequency error may be computed. The largest part of the probability density function of D is in the region of where the atan(x) ≈ x . Therefore,

PERFORMANCE COMPARISON TO THE CRAMER-RAO LOWER BOUND

The Cramer-Rao Lower Bound on the variance of the frequency error of any unbiased frequency estimator is given by,

The performance of the DFT based estimator may be compared to the Cramer- Rao Lower Bound.

For high SNR and large N, the performance frequency estimation variance is . above the Cramer-Rao Lower Bound.

Figure 6 shows in dB verses N, where N is the length of the FFT. JUSTIFICATION FOR ITERATION TOWARD A ZERO VALUE OF THE

DISCRIMINANT

The reason for the performance improvement of the proposed class of algorithms relative to prior algorithms is the first frequency interpolation allows the computation of two DFT coefficients, which are ½ DFT bin spacing above the first interpolated frequency and ½ DFT bin space below the first interpolated frequency. While the first

interpolation may still have significant error, which is dependent on the relationship of the true frequency relative to the FFT coefficient frequencies, the error discriminant evaluated for the first interpolated frequency will have a value close to zero. The variance of the frequency error is relatively low in the region of small values of the frequency discriminant. Therefore, the second interpolated frequency will have small error variance. There is significant noise performance advantage in using the first interpolation to allow a low error variance second interpolation. The interpolation may be iterated, with diminishing improvements of estimation accuracy, until convergence to a fixed point solution is obtained. Figure 6 shows the variance of the normalised estimator frequency error estimate vs the frequency error for the first interpolation. For the figure, N=64 and the signal to noise ratio is 6 dB. There is a very sharp reduction the rms error of the frequency estimator in the region of the frequency being close to the center frequency of the frequency discriminator. This indicates that tremendous improvement in performance obtained by iteration.

ALGORITHM SIMPLIFICIATONS RESULTING IN LARGE REDUCTIONS IN

COMPUTATIONAL COMPLEXITY

Simulation results have verified a single FFT and two iterations involving the computation of the discriminant function and are sufficient to obtain RMS frequency error performance close to the computation of . The first iteration moves the discriminant towards zero and decreases of . The estimate resulting from the second iteration therefore results in small error variance of .

The algorithm has complexity of O(N log ₂ (N) +0(8N)=O[N log ₂ (N)] The algorithm has the same order of complexity as the original FFT for performance, which is very close to the Cramer-Rao Lower Bound.

The algorithm has the same order of complexity as the original FFT for performance, which is very close to the Cramer-Rao Lower Bound.

Iteration with ψ ₁ (D) = D yield results very close to ψ(D) and saves considerable computational complexity. The fixed point solution for the two cases is identical. Since the algorithm comes very close to convergence with two iterations, two iterations are sufficient for most applications. The performance improvement with additional iterations is small.

OTHER FREQUENCY DISCRIMINANTS WITH THE SAME NOISE PERFORMANCE

There are a number of discriminants, which have the same performance, when used iteratively to obtain the fixed point solution, as the previously introduced

discriminants. The noise performance is identical, for iteration, because the fixed point solution is identical.

This class of discriminants includes functional forms,

and in particular, and, where * denotes complex conjugate, and where Re[.] is the real part.

SIMULATION

Figures 8 to 13 show the rms frequency error performance of the algorithm vs SNR in dB, for N=2,4, 16,64,246, and 1024, respectively. Both the cases of one interpolation and two iterative interpolations are shown. The two interpolation case is essentially achieves the performance of the infinite interpolation case.

SPECTRAL SYMMETRY ITERATIVE (SSI) FREQUENCY ESTIMATION

ALGORITHM

This technique of exploiting the said symmetry discrete spectral distribution to estimate the frequency of a complex exponential is further applicable to a wide class of signals that are bandpass in nature. Consider a general bandpass signal s[n] which is represented by:

where d _k is the corresponding k ^th discrete Fourier coefficients of the signal s[n] which in general can be expressed as d _k = a _k exp(jθ _k ) where a _k and θ _k corresponds to the magnitude and phase of the k ^th discrete Fourier coefficients respectively. k is the index of the discrete spectral distribution and has a range from 0 to

N-1 k _e is the value of k that corresponds to center frequency of the signal s[n] . / _e is the center frequency of the signal s[n]

f _k is the frequency corresponding to the discrete Fourier coefficients of the signal s [»] θ _e is the phase of the signal s[n] .

Now the baseband representation of the signal s[n], denoted as x[n], can be expressed without the knowledge of the phase of the carrier as:

A new signal, denoted as z[n] , is formed by multiplying s[n] in (24) and the complex conjugate of x[n] in (25) which can be mathematically expressed as:

where ψ can be any real value constant.

Denote the difference between any two θ _k 's as θ _α,β = θ _α -θ _β = -(θ _β -θ _α )= -θ _β,α , then multiplying any two of the said d _k 's yields:

First, consider the case where f _e is a integer multiple of 1/(NT _s ) , then k _e , which is the value of k that corresponds to frequency f _e , will be of an integer value. Taking the DFT

of z[n] in (26) and substituting the definition in (27) of multiplying any two of the said d _k 's yields the discrete spectral distribution Z[k] , which can be expressed as:

Substituting the definition of k _α,β and θ _α,β into Z[k] in (28) yields:

Taking the magnitude of Z[k] in (29) yields:

From the above equation, it can be seen the magnitude of the said discrete spectral distribution Z [k] in (29) is even symmetric about the frequency correspond to DFT bin k _e . In this example, k _e was chosen to be an integer value for illustrative purpose only. In general, k _e can be of an integer or non-integer value and the described even symmetric property of the magnitude of the discrete spectral distribution Z [k:]| in (30) the will hold whether k _e is of an integer or a non-integer value.

Given a signal that has a center frequency of f _e and a copy of the baseband version of the said signal is known, z[ri\ can be formed and the frequency estimation algorithm disclosed can be applied on z[n] in (26) to estimate the frequency f _e in the signal. Whilst the estimation algorithm disclosed above has ε , the distance at which to evaluate the modified DFT coefficient from the initial estimate of k _e , equal to plus and minus FFT bin. This value of ε was proven to be optimal in both accuracy and convergence rate for the case of estimating the frequency of a complex exponential.

However this value of ε may not apply for other passband signal types and maybe adjusted to other values depending on the type of passband signal the algorithm is applied to in order to obtain the frequency estimate.

The SSI frequency estimation algorithm can be seen as a generalisation of the algorithm in accordance with the first aspect if one considers s[n] in (24) as r[n] in

(1), x[n] in (25) to have the value of 1 for all n 's and ψ in (26) equals to 1 when constructing z[n] in (26). This generalisation allows frequency estimation of any general passband signal about a center frequency.

Another way of obtaining a spectral symmetric (both in phase and magnitude) spectrum from a bandpass signal, y[n] having a baseband component is as follows: z[n]= R _y [n] x R _x * [n] where,

R _y [n] is the autocorrelation of the received signal y[n] , R _x ^* [n] is the autocorrelation of the complex conjugate of the baseband signal x [n]

and y[n] and x[n] has the following relationship:

where Δf is the relative carrier frequency offset between the transmitter and receiver, Δθ is the relative carrier phase offset between the transmitter and receiver. The autocorrelation R _y [n] and R _x [n] is mathematically defined as:

Computing the said autocorrelation of y[n] and x[n] can be computationally intensive. An example of a fast method of computing the autocorrelation of y[n] and x[n] is as follows:

Assuming the signal set y[n] and x[n] has a length of N.

Step 1) Zero pad the respective signal set to length 2N by appending a vector of zeros of length N to the end of the said signal set.

Step 2) Perform a length 2N DFT of the zero padded signal set. Denote the output of the DFT as S _2p [k] .

Step 3) Calculat

Step 4) Perform a length 2N inverse DFT (IDFT) on M [k] . Denote the output as m[n].

Step 5) The autocorrelation of the respective signal set is then the first N samples of m[n] . That is

Using z[n] formed by the autocorrelation R _y [n] and R _x [n] for carrier frequency estimation has the advantage that the averaging performed by the autocorrelation of y[n] and x[n] minimizes the effect of fading channel has on the performance of the SSI frequency estimator. It also reduces the adverse effect that sample and frame timing offset has on the accuracy of the frequency estimate.

EXAMPLES OF APPLICATIONS OF THE SSI FREQUENCY ESTIMATION ALGORITHM

Phase Estimator

Once the frequency estimate has been established by the SSI frequency estimation algorithm the phase of the signal can be estimated by the relationship:

where θ - denotes the estimated phase of the said complex exponential given the said frequency estimate and Re[.] denotes the imaginary and real components of respectively. For the case of the signal being a complex exponential, replace z[n] with r[n] as described.

MODIFIED SPECTRAL SYMMETRY ITERATIVE (MSSI) FREQUENCY

ESTIMATION ALGORITHM

In the frequency estimation algorithm as disclosed in the first aspect, the number of samples points, N, used to obtained the frequency estimate was kept constant for all iterations of the said frequency estimation algorithm. Assuming that the signal to noise (SNR) ratio is sufficiently high, the number of samples used to obtain the initial discrete spectral distribute for the peak search and the first iteration of the said frequency estimation algorithm can be reduced. However, the Cramer Rao Lower Bound (CRLB) which is the lowest possible bound for estimation error variance, monotonically decreases with the increasing number of samples used for the estimation at a fixed SNR. Hence, reducing the number of samples used for the frequency estimate may not produce an estimate that will meet the required accuracy.

Since frequency estimation algorithm as disclosed in the first aspect requires performing a FFT operation on the signal sample for the initial frequency estimate, it is advantages to keep the size of the FFT as small as possible. Denoting N ₀ as the reduced number of sample points used to compute the FFT for the initial frequency estimate, and N as the number of samples required producing an estimate that meets the required accuracy, the frequency estimate can be found in two set of iterations of said frequency estimator with each set of iterations consist of at least two iterations. The first set of iterations will operate on N ₀ samples points and the last set of iterations will operate on N sample points. No degradation to the error variance performance will result when operating at SNR higher or equal to η if the following condition is met:

where η is lowest expected signal to noise ratio of the operating environment.

If the above condition is not met but it is necessary to keep N ₀ as low as possible and N relative large for the require accuracy, intermediate set of iterations can be inserted to bridge between the first and last set of iterations. Denoting N ₁ as the number of samples used in the intermediate set, the frequency estimate can be found in three sets of iterations by the said frequency estimator, with the first set of iterations operated on

N ₀ samples points and the second set of iterations operated on N ₁ sample points and the last set of iterations operated on N sample points, without any degradation to the error variance performance when operating at SNR higher or equal to η if the following condition is met:

Notice the number of intermediate sets of iteration can be increased to further reduce the number of samples, No, required for computing the FFT for the initial frequency search. In general, at any particular SNR, is it possible to determine a sufficiently large N ₀ , N ₁ ,

N ₂ such that the described procedure may be implemented with negligible degradation in error variance performance compare to the algorithm implemented with N samples.

To summarise the steps of the modified SSI frequency estimation algorithm. Denote z[n] as set the samples from which to estimate the frequency from.

Ste N ₀ point complex FFT is computed. The FFT output coefficients are d ,enot .ed , as w ,here

Step 2) Perform a peak search o o find which defined as:

Step 3) Set the initial frequency estimate as .

Step 4) Recursion starts at i = 0

Step 5) Perform the first set of iterations which consist of the following steps:

(ii) Recursion is started at m = 0 (iii) Compute the modified DFT coefficients for the m:th frequency estimate:

(iv) The frequency discriminant is then computed for any of the functional form as a function of and as presented in reference [I]. (v) The m+1 :th frequency estimate, f _m , is computed: f _m+1 = f _m + Δf _m (r ) .

Iterate steps (iii) to (v) until either the number of iterations has reached a predefined number or the output of the discriminant is below a predefined threshold.

Step 5) Set the value of equal to the last value of f _m in step 4)(v). Repeat step 4) for . The number of times step 4) is repeated depends upon the chosen value of ° and N _as described in the above description.

This modified version of the SSI frequency estimation algorithm can be applied to the frequency estimation of a complex exponential if one replace with in

(1). Figure 16 shows a comparison of the frequency estimation error variance verus

SNR in dB between the modified SSI frequency estimation algorithm against the frequency estimation algorithm as described in the first aspect of the invention. Both of the said frequency estimation algorithms were used to estimate the frequency of a complex exponential. N ₀ = 128 for the modified frequency estimator and N = 512 for both frequency estimators. The frequency of the complex exponential was randomly generated in each trial from a uniform distribution that spans from -1/(2T _s ) to 1/(2T _s ) with its phase randomly generated from a uniform distribution that spans from 0 to 2π and it has an amplitude of one. The results presented were averaged across 5000 trials. The solid line represents the CRLB for frequency estimation of a complex exponential, the line with the '*' is the simulated frequency error variance performance of the modified SSI frequency estimator, and the line with the is the simulated frequency error variance performance of the frequency estimator as described in the first aspect of the invention.

Figure 17 illustrates a comparison of the phase estimate error variance verus SNR in dB between the described phase estimator using the frequency estimate from the modified SSI frequency estimation algorithm and the said phase estimator using the frequency estimate from the frequency estimation algorithm as described in the first aspect of the invention. Both of the said phase estimation algorithms were used to estimate the phase of a complex exponential. N ₀ = 128 for the modified frequency estimator and TV" = 512 for both frequency estimators. The frequency of the complex exponential was randomly generated in each trial from a uniform distribution that spans

from -1/(2T _s ) to l/{2T _s ) with its phase randomly generated from a uniform distribution that spans from 0 to 2π and it has amplitude of one. The results presented were averaged across 5000 trials. The solid line represents the CRLB for phase estimation of a complex exponential, the line with the '*' is the simulated phase error variance performance using the frequency estimate from the modified SSI frequency estimator, and the line with the ' D ' is the simulated phase error variance performance using the frequency estimate from the frequency estimator as described in the first aspect of the invention.

Figure 18 illustrates a comparison of the amplitude estimate error variance verus SNR in dB between the described amplitude estimator using the frequency estimate from the modified SSI frequency estimation algorithm and the said amplitude estimator using the frequency estimate from the frequency estimation algorithm as described in the first aspect of the invention. Both of the said amplitude estimation algorithms were used to estimate the amplitude of a complex exponential. N ₀ =128 for the modified frequency estimator and N = 512 for both frequency estimators. The frequency of the complex exponential was randomly generated in each trial from a uniform distribution that spans from -1/(2T _s ) to 1/(2T _s ) with its phase randomly generated from a uniform distribution that spans from 0 to 2π and it has an amplitude of one. The results presented were averaged across 5000 trials. The solid line represents the CRLB for amplitude estimation of a complex exponential, the line with the '*' is the simulated amplitude error variance performance using the frequency estimate from the modified SSI frequency estimator, and the line with the 'G' is the simulated amplitude error variance performance using the frequency estimate from the frequency estimator as described in the first aspect of the invention.

CHANNEL AND FREQUENCY ESTIMATION FOR OFDM SYSTEMS One of the major impairment wireless communication system such as orthogonal frequency division multiplex (OFDM) encounters whilst operating in urban or indoor environment is multipath fading effect. Multipath fading causes distortion in the spectrum of the OFDM signal and causes significant degradation to the bit error rate (BER) performance.

In some OFDM standards such as the ETSI standard for digital video broadcasting (DVB-T) [see reference 2], pilot symbols are inserted amongst data symbols for channel estimation and carrier frequency offset (CFO) acquisition purposes. There are number of techniques in the literature that can acquire the integer component of the CFO (i.e. CFO that are integer multiples of 1/(NT _s ) ) using the said pilot symbols

[see references 3 and 4].

In other OFDM standards, such as the IEEE 802.11 series, training symbols are transmitted for framing synchronisation purposes and there exists a number of methods that can acquire the integer CFO from these training symbols [see references 5 and 6].

In either cases, once the integer component of the CFO has been estimated accurately, the channel characteristic and the fraction CFO (i.e. the remaining CFO that is bounded between -1/(2NT _s ) and 1/(2NT _s )) can be found using the following iterative algorithm.

ORTHOGONAL FREQUENCY DIVISION MULTIPLEXING MODULATION (OFDM) ALGORITHM 1

Denoting r _CF0 [n] as the received OFDM symbol subjected to a fading channel and CFO. r _cpo [n] can be mathematically expressed as:

where d _k is the modulation of the k ^th data carrier of the transmitted OFDM signal h _k is the k ^th DFT coefficient of the channel frequency response (CFR)

Af is the frequency of the CFO Aθ is the phase of the CFO f _k is the k ^th frequency corresponding to the k ^th DFT bin w _h [n] is the sequence of noise at the output of the channel.

Step 1) Perform a coarse frequency estimate on r _CF0 [n] using any integer CFO estimator.

Step 2) Compensate r _CF0 [n] with the available CFO estimate. Denote r _comp [n] as the output of the compensator.

Step 3) Estimate the channel characteristic from r _comp \n\ either in the form of the channel frequency response (CFR) or the channel impulse response (CIR) using existing channel estimation techniques. If there is prior knowledge on the channel characteristics, this step can be skipped.

Step 4) Equalization is then performed on r _comp [n] with the known or estimated channel characteristic. Assigning r _eq [n] as the output of the equaliser.

Step 5) If there is no prior knowledge on the data modulation, d _k , of the transmitted OFDM signal, detect the modulation of the data carriers, denoted as d _k , from r _eq [n] . Form r _est [n] with d _k and the known or estimated channel characteristics either by:

or

where is the known or estimated coefficients of the CFR at the k ^th DFT bin is the set of known or estimated coefficients of the CIR

* denotes the discrete convolution operator.

The symbols may be detected by making threshold decisions on the r _est [n] sequence. The sequence may then go into an error detection decoder, which will detect a valid code word. The decoded data may then be re-encoded to form a higher reliability sequence of r _es t[n] . This enhances the reliability of r _est [n] that is used in step 6.

If there is prior knowledge on the data modulation, d _k , of the transmitted signal

(like the case of transmitting a training symbol), Form r _est [n] with the known data symbols, d _k , and the known or estimated channel characteristics either by:

^or

Step 6) Form z[n] = r _comp [n]r _est * [n] where r _comp [n] is the same as the output of step 2).

Step 7) Apply the SSI frequency estimation algorithm or its variant on z[n] to find the CFO estimate

Step 8) Iterate steps 2 - 7 until either the number of iterations has reached a predefined number or the output of the discriminant of the said frequency estimation algorithm is below a predefined threshold. Update the estimated CFO in step 2) with the CFO estimate in step 7) in each iteration. Figure 20 shows the frequency estimate error variance verus SNR in dB of the described "OFDM algorithm 1" employing the SSI frequency estimation algorithm. The said CFO estimation algorithm was used to estimate the fractional CFO of the received OFDM symbol. The structure of, the OFDM symbol follows the DVB-T standard as described in [2]. The data symbols are randomly modulated using QPSK modulation scheme. The channel model follows the fixed Rayleigh channel model as stated in [2]. No prior knowledge of the channel characteristics or the modulation of the data carriers used in the estimation of the CFO in this simulation. The CFO was randomly generated in each trial from a uniform distribution that spans from -1/{2NT _s )

to l/(2NT _s ) with its phase randomly generated from a uniform distribution that spans from 0 to 2π . The results presented were averaged across 5000 trials. The solid line represents the CRLB for frequency estimation of a complex exponential, the line with the ,'*' is the simulated result of the CFO error variance performance using the SSI frequency estimation algorithm.

If there is prior knowledge on the channel characteristic and transmitted data symbol, the following algorithm can be used to estimate the integer and fractional CFO.

OFDM ALGORITHM 2 Step 1) Form r _est [n] with the known data modulation, d _k , of the k ^th data carrier and the known channel characteristics either by:

where h _k is the known coefficients of the CFR at the k ^th DFT bin h[n] is the set of known or estimated coefficients of the CIR

* denotes the convolution operator.

Step 2) Form z[n] = r _CFO [n]r _est * [n] where r _CFO [n] as defined in OFDM Algorithm I.

Step 3) Denoting Z[k] as the FFT z[n] . Perform a peak search on |Z[k]|to find k _max which defined as: k _max = max[|Z [k]| : 0 < k ≤ N ₀ - 1] . Step 4) Set the initial CFO estimate as where N is the number of samples in z[n] .

Step 5) Using as the initial CFO estimate, apply the SSI frequency estimation algorithm or its variant on z [n] to find the CFO estimat

If channel is of additive white Gaussian noise (AWGN) nature, and the transmitted data symbols are unknown, the following algorithm can be used to estimate the fractional CFO.

OFDM ALGORITHM 3 Since the channel characteristics of AWGN channel is flat, r _CFO [n] can be mathematically expressed as:

where d _k , Af , Aθ and f _k are as defined in OFDM Algorithm I and w[n] is the sequence of the AWGN.

Step 1) Perform a coarse frequency estimate on r _CF0 [n] using any integer CFO estimator.

Step 2) Compensate the r _CF0 [n] with the available CFO estimate. Denote ^r _comp [ ⁿ ] ^{as me} output of the compensator. Step 3) Detect the modulation of the data carriers, denoted as , from r _comp [n] .

Form r _est [n] with by:

The symbols may be detected by making threshold decisions on the r _est [n] sequence. The sequence may then go into an error detection decoder, which will detect a valid code word. The decoded data may then be re-encoded to form a higher reliability sequence of r _est [n]. This enhances the reliability of r _est [n] that is used in step 4.

Step 4) Form z [n] = r _comp [n] rj [n] .

Step 5) Apply the SSI frequency estimation algorithm or its variant on z[n] to find the CFO estimate

Step 6) Iterate steps 2 - 5 until either the number of iterations has reached a predefined number or the output of the discriminant of the said frequency estimation algorithm is below a predefined threshold. Update the estimated CFO in step 2) with the CFO estimate in step 5) in each iteration. If channel is of AWGN nature, and the transmitted data symbols are known, the following algorithm can be used to estimate the integer and fractional CFO.

OFDM ALGORITHM 4 Step 1) Form r _est [n] with the known data modulation, d _k , of the k^ data carrier by:

Step 2) Form z[n] = r _CFO [n]r _est * [n] where r _CF0 [n] is as defined in OFDM Algorithm 3.

Step 3) Denoting Z[k] as the FFT z[n] . Perform a peak search on |Z[k]| to find k _max which defined as: k _max = max [|z [k]| : 0 ≤ k ≤ N ₀ - 1] .

Step 4) Set the initial CFO estimate as , where N is the number of samples in z[n] . Step 5) Using as the initial CFO estimate, apply the SSI frequency estimation algorithm or its variant on z [n] to find the CFO estimate

Figure 19 illustrates a comparison of the CFO estimate error variance verus SNR in dB between the described "OFDM algorithm IV" employing the modified SSI frequency estimation algorithm and the "OFDM algorithm IV" employing the SSI frequency estimation algorithm. Both of the said CFO estimation algorithms were used to estimate the integer and fractional CFO of the received OFDM symbol. N ₀ = 32 for the modified frequency estimator and N = 64 for both frequency estimators. The

OFDM symbol used is the long training symbol as described in the IEEE 802.11a standard [7]. The CFO was randomly generated in each trial from a uniform distribution that spans from -1/(2T _s ) to 1/(2T _s ) with its phase randomly generated from a uniform distribution that spans from 0 to 2π . The results presented were averaged across 5000 trials. The solid line represents the CRLB for frequency estimation of a complex exponential, the line with the '*' is the simulated result of the CFO error variance performance using the modified SSI frequency estimator, and the line with the ' D' is the simulated result of the CFO error variance performance using the SSI frequency estimation algorithm.

FREQUENCY AND PHASE ESTIMATION OF ULTRAWIDE BAND (UWB)

SIGNAL The SSI frequency estimation algorithm can be applied to estimating the frequency offset of the received UWB signal given a known training signal was transmitted. The said phase estimator can also be applied to estimate the phase offset of the received UWB signal.

Figure 21 to figure 26 illustrates an example of the effect of frequency and phase offset has on UWB signals and the effectiveness of the SSI frequency estimation algorithm and the phase estimator of estimating the said frequency and phase offset.

Figure. 21 and figure 22 show the real and imaginary component of the UWB pulse in the noiseless case. Figure 23 and figure 24 show the real and imaginary component of the same UWB pulse subjected to a frequency and phase shift. In this illustration, the number of samples of the UWB pulse N = 24, the frequency shift amounts to 3/{2NT _s ) Hz and the phase shift is equal to π/2.5 radians. Figure 25 and

figure 26 show the real and imaginary component of the recovered UWB pulse by compensating the frequency and phase shifted version of the pulse with the frequency and phase estimate from the said SSI frequency estimation algorithm and phase estimator. Figure 27 to figure 32 illustrate example of estimating and compensating for the frequency and phase offset of a noise corrupted UWB pulse. Figure 27 and figure 28 show the real and imaginary component of the UWB pulse at SNR equal to 5 dB. Figure 29 and figure 30 show the real and imaginary component of the noise corrupted UWB pulse subjected to the same frequency and phase shift as in the example above. Figure 31 and figure 32 shows the real and imaginary component of the recovered noise corrupted UWB pulse by compensating the frequency and phase shifted version of the pulse with the frequency and phase estimate from the said SSI frequency estimation algorithm and phase estimator.

Figure 33 to figure 38 illustrate another example of estimating and compensating for the frequency and phase offset of a noise corrupted UWB pulse as described above except the SNR is not equal to 0 dB.

Figure 39 illustrates the frequency estimate error variance verus SNR in dB of the SSI frequency estimation algorithm. The said frequency estimation algorithm was used for estimating the frequency offset of a frequency and phase shifted UWB pulse. Prior knowledge on the transmitted pulse is assumed in this simulation. The frequency offset was randomly generated in each trial from a uniform distribution from -1/(2 NTj to l/(2NT _s ) with the phase offset randomly generated from a uniform distribution from 0 to 2π . The results presented were averaged across 5000 trials. The solid line represents the CRLB for frequency estimation of a complex exponential, the line with the '*' is the simulated result of the frequency offset error variance performance using the SSI frequency estimation algorithm.

Figure 40 illustrates the phase estimate error variance verus SNR in dB of the phase estimator employing the frequency estimate from SSI frequency estimation algorithm. The said phase estimator was used to estimate the phase offset of a frequency and phase shifted UWB pulse. Prior knowledge on the transmitted pulse is assumed in this simulation. The frequency offset was randomly generated in each trial from a uniform distribution that spans from to l/(2NT _s ) with the phase offset randomly generated from a uniform distribution that spans from 0 to 2π . The results presented were averaged across 5000 trials. The solid line represents the CRLB for phase estimation of a complex exponential, the line with the '*' is the simulated result of the phase offset error variance performance using said phase estimator.

FREQUENCY ESTIMATION OF A SIGNAL COMPOSED OF A PLURALITY OF

COMPLEX EXPONENTIALS

In the following embodiment it is assumed that the number of complex exponentials in the signal are known. Denote x[n] as the set of samples from which to perform estimation of each of the parameters, frequency, phase and amplitude. Mathematically x[n] can be expressed as:

where m is the index of the m ^th complex exponential;

M is the number of complex exponentials; A ₁₁₁ is the amplitude of the m ^th complex exponential; f _m is the frequency of the m ^th complex exponential (in Hertz); θ _m is the phase of the m ^th complex exponential (in Radians); w[n] is the sequence of a complex additive white Gaussian noise (AWGN);

T _s is the sampling period (in seconds); and n is the index of the time samples and has a range from 0 to N-I where N is the number of samples.

In this embodiment a first estimation referred to as a coarse estimation is performed. This is followed by a fine estimation to refine the parameters of each complex exponential.

Coarse estimation:

Step 1) Recursion is started at m = 0.

Step 2) Set x _m [n] = x [n] where m = 0.

Step 3) An N point complex FFT is computed on x _m [n] . The FFT output coefficients are denoted as X _m [k] where 0 ≤ k ≤ N-l .

Step 4) A peak search is performed on X _m [k]| to find k _m which defined as: k _m = max[|X _M [k]|: 0 ≤ k ≤ N-1] . Step 5) The initial frequency estimate is set a

Step 6) A frequency estimation algorithm (as later described) is applied on x _m [n] to refine on the initial frequency estimat . This refined frequency estimate is used to update .

Step 7) A phase estimation algorithm (as later described) is applied on x _m [n] using the refined frequency estimate from step 6). The estimated phase is denoted as .

Step 8) An amplitude estimation algorithm (as later described) is applied on x _m [n] using the refined frequency estimate from step 6). The estimated amplitude is denoted a Step 9) Compute are the outputs from step 4), step 5) and step 6) respectively. Step 10) Iterate step 3) to step 9) for m = 1 to M -1.

Fine estimation:

Step 1) Recursion is started at q = 0 . Step 2) Set and for m = 0 to M -1 where and are the frequency, phase and amplitude coarse estimates of the M complex exponentials.

Step 3) Recursion is started at m = 0.

Step 4) Compute . Step 5) Using f _{q m} as the initial frequency estimate, the frequency estimation algorithm (as later described) is applied on x _m [n] to refine on the initial frequency estimate The output of the said frequency estimation algorithm is denoted as Step 6) The phase estimation algorithm (as later described) is applied on x _m [n] using the frequency estimate from step 5). The output of the said phase estimation algorithm is denoted as .

Step 7) The amplitude estimation algorithm (as later described) is applied on x _m [n] using the frequency estimate from step 5). The output of the said amplitude estimation algorithm is denoted as .

Step 8) Iterate step 4) to step 7) for m = 1 to M -1. Step 9) Iterate step 3) to step 8) for q = 1 to Q where Q is the number of iterations required for all the parameter estimations to reach convergence. Convergence is defined as the value of the residual error, is below of predefined threshold.

In a further embodiment the above method of estimating the frequency of a signal composed of a number of complex exponentials is parallelised so that the fine estimation of parameters of the M complex exponential are spawn as separate processes.

The advantage here is that computational time required to obtain the estimates is reduced. Parallelisation involves the following steps:

Step 1) Wait for all the processes to be free.

Step 2) Fetch the incoming x[n] from the main queue and store it into a shared memory bank of which all processes have access.

Step 3) Perform the coarse estimation as described above to determine the estimates for the parameters .

Step 4) Store the initial coarse parameter estimates into a shared memory bank of which all processes have access.

Step 5) Indicate to the individual processes that the initial coarse parameter estimates and the new x[n] is ready. Step 6) Mark all processes to be busy.

For each of the M individual process do the following:

Step 1) Wait for the ready flag from the coarse estimation process. When ready, fetch a copy of x[n] from the shared memory bank and store it into own memory bank.

Step 2) Recursion is started at q = 0.

Step 3) Wait for all the parameter estimates to be updated. the all become available, fetch a copy of all the updated parameter estimates and for m = 0 to M -1 from the share memory bank and store then into own memory ban Step 4) Set and where m denotes the index of the m complex exponential this process is estimating the parameters for. are the frequency, phase and amplitude estimates obtained in step 3).

Step 5) Compute .

Step 6) Using as the initial frequency estimate, apply the frequency estimation algorithm (as later described) on x _m [n] to refine on the initial frequency estimate f _{q m} . Denote the output of the said frequency estimation algorithm as

Step 7) Apply the phase estimation algorithm (as later described) on x _m [n] using the frequency estimate from step 6). Denote the output of the said phase estimation algorithm as . Step 8) Apply the amplitude estimation algorithm (as later described) on x _m [n] using the frequency estimate from step 6). Denote the output of the said amplitude estimation algorithm as .

Step 9) Update and , in the shared memory bank with the values ^{and obtained in ste} P ⁶ )> ^ste P ⁷ ) ^{and ste} P ⁸ ) respectively.

Step 10) Indicate to other process that the parameter estimates corresponding to the m ^th complex exponential has been updated.

Step 11) Iterate step 3) to step 10) for q = 1 to Q where Q is the number of iterations required for all the parameter estimations to reach convergence. Convergence is defined as the value of is below of predefined threshold.

Step 12) Indicate to the coarse estimation process the fine estimation process for the current m ^th complex exponential is complete and ready for the next set of samples.

Frequency estimation algorithm

Denot as the current frequency estimate to be refined and z[n] as the set of signal samples from which to estimate the frequency from. Step 1) recursion is started at r = 0. Step 2) the DFT coefficients for the r ^th frequency estimate are computed:

Step 3) The frequency discriminant is then computed:

Step 4) The r ^{th +1} frequency estimate, f _m+1 , is computed:

Convergence for the frequency estimate is reached if is sufficiently small. If there is convergence, then the frequency estimate is . If convergence for the frequency estimate has not been reached, then r is incremented by 1 and ste to 4 are repeated. In practice, the algorithm will converge for r = 2. Therefore, only for r = 0, 1, a need to be computed (i.e. two iterations) which means that the frequency estimate is .

Phase estimation

Denote as the current frequency estimate and z[n] as the set of signal samples from which to estimate the frequency from.

The phase of the signal can be estimated by the relationship:

where denotes the estimated phase of the complex exponential given the frequency estimate Im[.] and Re[.] denotes the imaginary and real components of respectively.

Amplitude estimation

The amplitude of the complex exponential can be estimated using the following relationship:

Denote and as the current frequency estimate and phase estimate respectively. z[n] as the set of signal samples from which to estimate the frequency from.

where denotes the estimated amplitude of the said complex exponential given the said frequency estimate Im[.] and Re[.] denotes the imaginary and real components of respectively.

PERFORMANCE ANALYSIS BY SIMULATION

The signal to noise ratio (SNR) for all the results illustrated in figures 41 to 53 is defined as the ratio of the total power of all complex exponentials present in the signal against the total noise in the signal.

The disclosed method relies on subtracting estimated exponential components from the original signal to mitigate the ICI caused spectral leakage components of the individual complex exponentials. Assuming the amplitudes of the individual complex exponentials are approximately equal, this implies the resulting signal vector used for

the fine estimation process will have a SNR degradation that is approximately proportional to the number of complex exponentials present in the signal.

Denoting SNR _S!M as the SNR used to generate the AWGN which can be mathematically defined as:

Then the above statement implies:

where SNR _EFF is the effective SNR the fine estimation process sees in the signal and M is the number of complex exponentials in the original signal.

The error variance performance of the method for signals having different frequency spacing are illustrated in Figures 41 to 53. The SNR used for calculating the

Cramer-Rao lower bound (CRLB) for frequency, phase and amplitude estimation is defined by SNR _EFF . The channel model is assumed to be of flat non-fading in nature.

With reference to figure 41, the signal contains seven complex exponentials and their respective frequencies are equally spaced at 3/{2NT _s ) Hz. The number of samples./V = 512. A random frequency was generated in each trial from a uniform distribution that spans from ~1/(2NT _s ) to 1/{2NT _s ) and all the complex exponentials are frequency shifted by this random frequency. A random phase was generated in each trial from a uniform distribution that spans from 0 to 2π for each of the individual complex exponentials. The amplitudes of the complex exponentials are equal to 1. The results presented were the averaged frequency error variance of the individual complex exponential. The results were averaged over 1000 trials. The solid line represents the

Cramer-Rao lower bound (CRLB) for frequency estimation of a complex exponential, the line with the '*' is the simulated result of the frequency error variance performance.

With reference to figure 42, the simulation parameters are identical to those described for figure 2. The solid line represents the CRLB for phase estimation of a complex exponential, the line with the '*' is the simulated result of the phase error variance performance.

Referring now to figure 43, the simulation parameters are identical to those described for figure 41. The average power of the complex exponentials were taken for calculating the CRLB for amplitude estimation. The solid line represents the CRLB for amplitude estimation of a complex exponential, the line with the '*' is the simulated result of the amplitude error variance performance.

Referring now to figure 44, the simulation parameters are identical to those described for figure 41 except the frequencies of the complex exponentials are equally spaced at 2/(NT _s ) Hz. The solid line represents the CRLB for frequency estimation of a complex exponential, the line with the '*' is the simulated result of the frequency error variance performance.

With reference to figure 45, the simulation parameters are identical to those described for figure 41. The solid line represents the CRLB for phase estimation of a complex exponential, the line with the '*' is the simulated result of the phase error variance performance. Referring now to figure 46, the simulation parameters are identical to those described for figure 44. The average power of the exponentials were taken for calculating the CRLB for amplitude estimation. The solid line represents the CRLB for phase estimation of a complex exponential, the line with the '*' is the simulated result of the amplitude error variance performance. Referring now to figure 47, the simulation parameters are identical to those described for figure 41 except the frequencies of the complex exponentials are equally spaced at 3/(NT _s ) Hz. The solid line represents the CRLB for frequency estimation of a complex exponential, the line with the '*' is the simulated result of the frequency error variance performance. With reference to figure 48, the simulation parameters are identical to those described for figure 47. The solid line represents the CRLB for phase estimation of a complex exponential, the line with the '*' is the simulated result of the phase error variance performance.

With reference to figure 49, the simulation parameters are identical to those described for figure 47. The average power of the exponentials was taken for calculating the CRLB for amplitude estimation. The solid line represents the CRLB for phase estimation of a complex exponential, the line with the '*' is the simulated result of the amplitude error variance performance.

Referring now to figure 50, the simulation parameters are identical to those described for figure 41 except the frequencies of the complex exponentials are equally spaced at 4/(NT ₁ ) Hz. The solid line represents the CRLB for frequency estimation of a complex exponential, the line with the '*' is the simulated result of the frequency error variance performance.

Referring to figure 51, the simulation parameters are identical to those described for figure 50. The solid line represents the CRLB for phase estimation of a complex exponential, the line with the '*' is the simulated result of the phase error variance performance.

Referring now to figure 52, the simulation parameters are identical to those described for figure 50. The average power of the complex exponential was taken for calculating the CRLB for amplitude estimation. The solid line represents the CRLB for phase estimation of a complex exponential, the line with the '*' is the simulated result of the amplitude error variance performance.

Figure 53 illustrates a comparison performance of the frequency error variance performance between the multiple complex exponential parameter estimation algorithm, the root MUSIC algorithm and the ESPRIT-TLS algorithm. The signal contains seven complex exponentials with frequencies given by the following vector [-123.4353 -24.4353 -20.343 28.4354 35.4343 128.32 134.4343] in Hz. The number of samples N = 512. A random phase was generated in each trial from a uniform distribution that spans from 0 to 2π for each of the individual complex exponentials. The amplitudes of the complex exponentials are given by the following vector [1.34243 1.4343 1.54354 1.6766 1.4343 1.2545 1.4656]. The results presented were the averaged frequency error variance of the individual complex exponential. The results were averaged over 1000 trials. The solid line represents the CRLB for frequency estimation of a complex exponential, the line with the '*' is the simulated result of the frequency error variance performance, the line with the ' D ' is the simulated result of the root MUSIC algorithm and the line with the 'Δ' is the simulated result of the ESPRIT-TLS algorithm.

Spectral Peak Detector

In the above embodiments it was assumed that the number of complex exponentials in the signal were known. In this embodiment a spectral peak detector is utilised to statistically determine whether the local maxima in the spectrum are frequency components or if they are simply small fluctuations due to noise in the spectrum.

There exist a number of algorithms in the literature for spectral peak determination such as spectral null hypothesis [9] and spectral deconvolution [10]. The purpose of these algorithms is to statistically determine whether the local maxima in the spectrum are frequency components or if they are simply small fluctuations due to noise in the spectrum.

We have incorporated these spectral peak determination techniques with the fine estimator in an iterative way to detect all the complex exponentials and estimate their respective frequencies, amplitude and phase as shown below:

Let M _άet be the number of detected complex exponentials.

Referring to figure 54, recursion begins with the peak search process 100, where using existing spectral peak determination techniques, the location of the spectral peaks

corresponding to the frequencies of the complex exponential and the number of complex exponentials are statistically determined.

If the peak search algorithm returns M _det > 0 , a coarse parameter estimation 200, can be performed using the frequencies estimates obtained by the peak search process with the following steps:

Step 1) Recursion is started at m = 0 , where m denotes the index of the m ^th detected complex exponential.

Step 2) Set x _m [n] = x[n] where x[n] is the current set of samples of the received signal. Step 3) A frequency estimation algorithm (as previously described) is applied on x _m [n] to refine on the initial frequency estimate of the m ^th complex exponential obtained from peak detection process. Update with the said refined frequency estimate.

Step 4) A phase estimation algorithm (as previously described) is applied on r I x _m \n\ using the refined frequency estimate from step 3). Denote the estimated phase as .

Step 5) An amplitude estimation algorithm (as previously described) is applied r i ^Λ on x _m \ n\ using the refined frequency estimate f _m from step 3). Denote the estimated amplitude as A _n . Step β) Compute and

A ₁₁₁ are the outputs from step 3), step 4) and step 4) respectively. Step 7) Iterate step 2) to step 6) for m = 1 to M _det -1.

These coarse parameter estimates are then sorted 300, in descending order according to the coarse amplitude estimates. The estimated parameters of complex exponentials with the largest amplitude will be listed first, than the second largest and so on. Since complex exponential with higher amplitude is statistically more likely to exist rather than caused by fluctuation in the noise, it is logical to have the parameter estimates of those components refined by the fine estimation process 400, first. These sorted estimated are then fed to the fine parameter estimation process 400, as described in the earlier embodiment. It is also possible to incorporate another spectral peak statistical test into the fine parameter estimation process to eliminate any false complex exponentials that was detected by the peak search process.

This process of peak search, coarse and fine parameter estimation can be repeated for L times. At the end of each iteration, the samples of the received signal, x[n], is updated

(500) according to the following rule:

where are the refined parameter estimates of the M _det complex exponentials. This updated version of x[n] is then fed back to the peak search so that complex exponentials with lower amplitudes that might have been missed by the peak search process in the first iteration can be redetected without interference from complex exponentials with higher amplitudes.

Note the stated parallel implementation can also be applied to this case with additional control logic between peak search 100, coarse parameter estimation 200, sorting the parameter estimates 300 and fine parameter estimation, to control the access of shared memory and the process flow.

It should be appreciated that embodiments of the invention can be applied in a variety of areas, including, but not limited to the fields of telecommunications, radar, sonar, electronic measurements and instrumentation, biomedical measurements, speech processing, image processing. It will be appreciated by persons skilled in the art that numerous variations and/or modifications may be made to the invention as shown in the specific embodiments without departing from the spirit or scope of the invention as broadly described. The present embodiments are, therefore, to be considered in all respects as illustrative and not restrictive.

References:

[1] Reisenfeld S., and Aboutanios E., "Frequency Estimation", Patent Application No. WO 2004/005945 Al, July 4 ^th , 2003.

[2] ETSI EN 300 744 Vl.1.2 (1997-08): "Digital broadcasting systems for television, sound and data services; framing structure, channel coding and modulation for digital terrestrial television," ETSI, 1997.

[3] Li M., and Zhang W., "A novel method of carrier frequency offset estimation for OFDM systems," IEEE Trans, on Consumer Electronics, Vol. 49, No. 4, pp 965- 972, December 2003.

[4] Speth M., Fechtel S., Fock G., and Meyr H., "Optimum Receiver Design for OFDM-Based Broadband Transmission— Part II: A Case Study," IEEE Trans, on Communications, VOL. 49, NO. 4, pp 571-578, April 2001.

[5] Dlugaszewski Z., and Wesolowski K., "Simple coarse frequency offset estimation schemes for OFDM burst transmission," IEEE International Symposium on Personal, Indoor and Mobile Radio Communications, Vol.2, pp 567 - 571, Sept. 2002.

[6] Chen C, and Li J., "Maximum likelihood method for integer frequency offset estimation of OFDM systems," Electronics Letters , Vol. 40, Issue. 13, pp 813 - 814, June 2004.

[7] IEEE 802.1 la-1999, Part 11: "Wireless LAN medium access control (MAC) and physical layer (PHY) specifications: High-speed physical layer in the 5 GHz Band," September 1999.

[8]. Huang D., "Approximate maximum likelihood algorithm for frequency estimation," Statitϊca Sinica, VoI 10, pp 157-171, 2000.

[9] Pervical D. B., and Walden A. T., "Spectral Analysis For Physical Applications: Multitaper and Conventional Univariate Techniques, " Cambridge University Press, Cambridge, 1993.

[10] Miekina A., Morawski R. Z., and Barwicz A., "The Use of Deconvolution and Iterative Optimization for Spectrogram Interpretation", IEEE Trans, on Instrumentation and Measurement, Vol. 46, No. 4, pp 1049-1053, August 1997.

[11] Li M., and Zhang W., "A novel method of carrier frequency offset estimation for OFDM systems," IEEE Trans, on Consumer Electronics, Vol. 49, No. 4, pp 965- 972, December 2003.

[12] Du Z., and Zhu J., "Improved coarse frequency synchronization algorithm with extended differential detection," WCNC 2003, Vol. 1 , pp 470-477, March 2003

[13] Macleod, M. D., "Fast Nearly ML Estimation of the Parameters of Real or Complex Single Tones or Resolved Multiple Tones," IEEE Transactions on Signal Processing, Vol. 49, No. 1, pp 141-148, January 1998.

[14] Quinn, B. G., "Estimating frequency by interpolation using Fourier coefficients," IEEE Transactions on Signal Processing, Vol. 42, No. 5, pp 2264-2268, May 1994.

[15] Fitz M.P., "Planar Filtered Techniques for Burst Mode Carrier Synchronization," IEEE Trans, on Instrumentation and Measurement, pp 365-368, August 1991.

[16] Schmidl T.M., and Cox D.C., "Robust frequency and timing synchronization for OFDM," IEEE Transactions on Communications, VoI 45, No. 12, December, 1997.