Shehrzad is a member of the Technical Staff for Picoliter Inc. He can be contacted at [email protected].
Time-delay estimation of waveforms is of fundamental importance in a variety of signal-processing applications. The problem boils down to computing the location of a digitized "signature" waveform residing within a larger time-slice. In Figure 1, for instance, I'm searching for the location of the triangular pulse waveform within the larger signal. In a typical signal-processing application such as radar or sonar, an "echo pulse" is emitted from a transmitter, and a receiver device records the (digitized) response. The transmitted ping bounces off objects, resulting in something like Figure 1(b). After the time-delay estimation (TDE) is calculated, the distance of the object can be derived by a simple formula, assuming the velocity of the wave field is known a priori.
In the absence of noise, or in a relatively high signal-to-noise ratio (SNR) environment, this problem can easily be solved. Consider two simplistic approaches, one where the estimator searches for the maximum signal value in the digitized response, and another where the estimator thresholds the response, and uses the first location where the response rises above the threshold as the time-delay. While easy to implement, both methods suffer from significant amounts of "jitter" in the presence of noise. In this article, I describe a commonly employed signal-processing tool—the matched filter—which solves the TDE problem, even in the presence of significant amounts of line noise. The matched filter is implemented using a mathematical formulation known as "cross-correlation."
As it turns out, Intel's Integrated Performance Primitives (IPP) includes just such a cross-correlation routine. IPP (http://www.developer.intel.com/software/ products/perflib/spl/index.htm) is a software API that provides a multitude of functions for multimedia, audio and video codecs, image processing, signal processing, speech compression, plus math support for lower-level routines like cross-correlation. In addition, IPP exposes a set of data types, data structures, and associated naming conventions so that these routines can be tailored for different types of input data (bytes, single-precision floating-point numbers, double-precision floating-point numbers, and so on). IPP runs on Intel processors such as the Pentium 4, Itanium architecture, Xeon, SA-1110, and PCA, and supports both Linux (32- and 64-bit) and Windows (32- and 64-bit). In this article, I use IPP and The Mathworks' Matlab programming environment (http://www.mathworks.com/) to illustrate how a matched filter, implemented using cross-correlation, can extract a signature waveform from a given signal—even when the signature waveform is buried within background noise. The line noise described in this simulation may not always be of the ambient variety; you can imagine a military setting where an enemy combatant is intentionally jamming a radar or sonar system, and the matched filter's job is to search the signal space, throw away (or ignore) the jamming noise, and extract whatever the system is interested in.
In addition, if you effectively model the noise, then the matched filter can be augmented with some form of adaptive filtering, further enhancing its effectiveness. In certain situations, matched filters can also be used in image processing (for example, to perform image registration) or to implement certain machine-vision primitives.
Noise and Estimation Jitter
The 64-sample waveform in Figure 1(b) is an example echo response. Figure 2 is a plot of three noisy echo responses, where I have added gaussian white noise (GWN) to the original, pristine echo response. The magnitude of the additive noise can be quantified by the standard deviation of the normal distribution added to the original signal; in this case I used standard deviations of 0.1, 0.2, and 0.3, respectively. I then computed time-delay estimations for 10 noisy input signals, for each example noise scenario, using both the max and threshold methods. Table 1 lists the simulation results.
As the level of noise increases, the uncertainty in the time-delay estimate increases. In particular, the thresholder performs terribly in the noisiest scenario. However, even the "max method" jitters around quite a bit in this scenario. Fortunately, there is another tool in our arsenal—the matched filter.
Cross-Correlation and Matched Filters
In its simplest incarnation, matched filters are built directly on top of the cross-correlation. The cross-correlation of two signals is a measure of similarity between the two sequences. The cross-correlation function (ccf) is maximized when the two signals have similar frequency content and are in phase with one another. In the continuous domain, the cross-correlation of two input signals x(t) and y(t) is defined as in Example 1(a).
Computers operate in the discrete domain, with finite sequences; hence, the cross-correlation of two digitized input sequences x[n] and y[n] is defined as in Example 1(b). These equations may look ominous, but are actually simple to visualize. The ccf is essentially the same operation as convolution, sans the time reversal of y(t) (or y[n]). Figure 3 is a graphical interpretation of cross-correlation. Performing the cross-correlation involves computing a series of shifted vector dot products, where the parameter l in Example 1(b) is typically referred to as the "lag" parameter.
If the input sequences y[n] and x[n] correspond to a cross-correlation "kernel" such as that given in Figure 1(a), and an input sequence (any of the three in Figure 2, for example), then a matched filter is nothing more than maximizing the ccf rx,y(l). That is, compute the ccf as given in Example 1(b), find the maximum correlation value, and then the output of the matched filter is the lag that this maximum corresponds to. This lag is the time shift where the two signals are the most similar, thus providing a robust method for solving the TDE problem.
The increase in performance when using a matched filter to solve the TDE problem is profound. To illustrate its robustness in the presence of significant amounts of noise, I reran the TDE experiment, this time adding GWN with standard deviations of 0.1, 0.2, 0.3, and 0.4. Table 2 shows the results. The column marked "xcorr" is the TDE as computed using the matched filter; note the stability of its results. The constant delta between the max/threshold estimates (24 versus approximately 40) and the matched filter estimate is to be expected, and has to do with the fact that the matched filter approach utilizes the entire "signature" waveform for its calculations, rather than a single point or threshold. This delta is always present and is dependent on the support (length) of the cross-correlation kernel. The salient point is that the matched filter is able to determine the location of the signature waveform every single time, even in relatively noisy cases. Figure 4 is a sample noisy waveform with GWN of standard deviation 0.4.
Implementation Issues
Not surprisingly, there are highly efficient algorithms for computing cross-correlations, in particular for input sequences of large length. Recall that cross-correlation is similar to convolution. Example 2 shows the relationship between cross-correlation and convolution, where the "*" symbol is the standard symbol denoting convolution.
Implementing convolution in software by extending the smaller input to be the same length as the larger input (by zero-padding), shifting y(l) for each lag, and then performing the vector dot product for each lag l is a brute-force and computationally prohibitive algorithm. Rather, DSP libraries and DSP hardware take advantage of one of the fundamental tenets of signal-processing theory—the time-convolution theorem. This theorem states that convolution in the time domain is equivalent to multiplication in the frequency domain. If the inputs to the cross-correlation are lengths that are powers of two (for example, the length of the sequence x[n] is 2k1 samples and the length of the sequence y[n] is 2k2 samples), then the fast fourier transform (FFT) algorithm can be used to transform the sequences from the time domain to the frequency domain. After applying the FFT to x[n] and y[n], thereby yielding X[f] and Y[f], you can multiply X[f] and Y[f] together, perform an inverse FFT, and the result is the cross-correlation. Example 3 is pseudocode for the fast cross-correlation algorithm.
The IPP Library includes cross-correlation routines for both signal-processing and image-processing. Listing One is a function that utilizes the signal-processing API function to compute the cross-correlation of two sequences. Listing Two is an example matched filter implementation, using the function in Listing One and a couple of other functions contained in the signal-processing API. Listing Three is a Matlab function that performs a cross-correlation of two sequences using Matlab's built-in convolution functionality. All of the Matlab source used to implement the simulation described here is available electronically; see "Resource Center," page 5. The Matlab Signal-Processing Toolkit (http://www.mathworks.com/) includes the function xcorr(), which can be used to perform cross-correlations.
DDJ
Listing One
#include "ipp.h"
bool xcorr(Ipp32f *x, int nx, // x[n]
Ipp32f *y, int ny, // y[n]
Ipp32f **r, int *nr, // rxy
int *startLag) // 1st lag we compute
{
// this function assumes that x & y are of power-of-2 lengths!
// the lags we care about
int lagLo = -(nx-1),
lagHi = ny-1;
// # of elements in the output array
*nr = lagHi-lagLo+1;
// allocate output array (use Intel SP function for mem allocation)
*r = ippsMalloc_32f(*nr);
IppStatus st = ippsCrossCorr_32f(x, nx, y, ny, *r, *nr, lagLo);
*startLag = lagLo;
return (st==ippStsNoErr);
}
Listing Two
int matched_filter(Ipp32f *x, int nx, // x[n]
Ipp32f *y, int ny) // y[n]
{
// matched filter, return value is the lag where
// the maximum cross-correlation by sequences x and y occurs.
// both of the following to be initialized by the function xcorr()
Ipp32f *crossCorrelation = 0;
int nxcorr = 0,
startLag = 0;
// perform the cross-correlation of x[n] & y[n]
xcorr(x, nx, y, ny, &crossCorrelation, &nxcorr, &startLag);
// | crossCorrelation |
ippsAbs_32f_I(crossCorrelation, nxcorr);
// max value of | crossCorrelation |
int maxLag = 0;
Ipp32f maxCrossCorr = 0.0;
ippsMaxIndx_32f(crossCorrelation, nxcorr, &maxCrossCorr, &maxLag);
return (maxLag+startLag);
}
Listing Three
function [l, rxy] = cross_corr(x, y) % % CROSS_CORR compute the cross-correlation of two sequences % [lags, rxy] = CROSS_CORR(x,y) returns the cross-correlation % of x[n] and y[n] in rxy. The output l contains the % lag indices. % % time-reverse x[n] x = fliplr(x); % construct lag indices l = 1-length(x):length(y)-1; % use built-in Matlab function conv to perform cross-correlation rxy = conv(y, x);
