A previous article (see archives) examined two ways to calculate the discrete autocorrelation function (ACF). This article builds on that concept by illustrating the ACF's usefulness. Indeed, an ACF has many applications in DSP but doesn't receive the same attention as the Fourier transform.

As a starting point, this discussion considers fractals. A fractal is an object or image that repeats itself at different scales; when you examine it in greater detail, you see a pattern similar to that visible at a coarser level. Nature uses fractals frequently for structures such as bronchial tubes or mountain textures. A common computer-related fractal that many readers are likely familiar with is the Mandelbrot set often used as a demo or benchmark for graphics cards.

In terms of signal processing, a fractal's attributes become interesting because a small piece of a signal exhibits similar characteristics to the entire fractal signal. Thus a small fractal can describe unbelievably complex signals. Examples of fractals that occur in DSP include communication-channel pattern errors and earthquake distributions.

To study fractals in more detail, let's quickly generate one based on the fact that colored noise generates fractals. But first a discussion of various types of noise should help some readers. Noise is either white or colored, and "coloredness" describes a noise signal's correlation. More specifically, white noise has no correlation from one sample to the next, is completely random and thus has a flat spectrum. In contrast, colored noise exhibits some correlation from one sample to the next and has a spectrum you can describe with an equation based on frequency. For instance, scientists often run into 1/f noise, whose spectrum falls off monotonically to zero as frequency increases. As a further note, don't confuse a noise signal's probability distribution with coloredness; those two aspects are independent. For example, either white or colored noise can have a Gaussian, uniform or other distribution.

Two popular methods

Given these definitions of colored noise, realize that many methods exist for generating it. This article concentrates on the two most popular: The first bases itself on Fourier representation, while the second uses the ACF. In the Fourier method, you start by considering the general shape of the spectrum's magnitude, for example 1/f. Given the magnitude in the frequency domain, you can take an inverse FFT to create a time-domain signal (colored noise), but before performing an inverse FFT, you also must specify the phase. You can set the phase to an arbitrary shape (the only limit is that its values lie between ±π/2, and by definition dc phase equals zero), perform the inverse transform and get colored noise -- each different phase waveform generates a different 1/f noise waveform in the time domain.

To see what 1/f noise looks like, compare the white-noise sequence in Fig 1a with the 1/f version based on the same phase values in Fig 1b. You might also be interested in examining other colored noises, specifically those that follow the formulas (1/f)² , (1/f)¹⁰ and (1/f)¹⁰⁰ as shown in Figs 1c-e, respectively. Comparing these plots, the 1/f version has the most self-similarity and thus is the most fractal-like.

Fig 1 -- Taking the phase characteristics of the white noise in (a), you can create colored noise using various formulas based on frequency (b-e). Notice how 1/f noise exhibits the most self-similarity and thus is the most fractal-like.

One problem with generating colored noise using the FFT method concerns signal length -- the more values you want in the time-domain waveform, the longer the FFT must run. In fact, sequences of moderate length going beyond several thousand values prohibit using this technique. Thus engineers sometimes turn to another popular method that applies random-signal analysis to the problem. Based on ACF, this alternative scheme produces coefficients for a special type of FIR filter such that when you feed that filter a series of random numbers (typically Gaussian with mean = 0 and arbitrary standard deviation), it outputs a sequence whose correlation characteristics can make it colored noise. This technique is possible because the ACF is directly related to the Fourier transform, which allows generating signals with a desired frequency spectrum such as the 1/f shape discussed earlier.

The first task in obtaining filter coefficients is to find the correlation matrix R that describes desired correlation characteristics. A trivial case exists with white noise where R(0) = 1 and every other correlation value equals 0. To get Rs for 1/f noise, you take magnitude values from the frequency domain and square them (in this case, resulting in (1/f)²). Taking the inverse FFT of these values generates ACF values that fill each column of the correlation matrix R according to the following scheme:

This correlation matrix R has special significance. In DSP, it forms the basis of many algorithms such as lattice and least-squares filters; mathematically, its inverse always exists if ACF is positive definite. Assuming the inverse exists, you can always decompose the matrix. Thus, taking advantage of some DSP theory and using a Cholesky decomposition (one efficient method of decomposing matrices), it's possible to rewrite the autocorrelation matrix into an upper and lower matrix as in:

The important fact is that a_k,l = a_0,k and h(n) = a_n,N. When you run a sequence of white-noise values through the FIR filter defined by h(n), the result is colored noise with the desired correlation and spectrum. Note that this method generates as many values as desired.

A major problem

One problem with this method is finding an ACF and the corresponding ACF matrix that's positive definite and that produces accurate results. Another way to state this difficulty is trying to find a time-domain sequence whose Fourier transform contains only real and non-negative values.

Several methods exist to get around this problem. One is to guess ACF values until you find a positive definite sequence. You can automate this procedure by using an optimization algorithm such as the simplex method. A second method assumes a sample noise signal is available, and you compute its ACF directly. Note that computed ACFs are automatically positive definite but might be quite long, which makes Cholesky decomposition more prone to error.

As an example, consider a typical correlation sequence I've encountered in practice: R(0) = 1, R(1) = -0.5, R(2) = 0.5, R(3) = -0.25 and R(4) = 0.125. (Note that these ACF values don't describe 1/f noise and illustrate that other colored noises are popular.) As described earlier, these values form a correlation matrix. Now the question is, how do you perform Cholesky decomposition? Rather than crank it out by hand, turn to a matrix-processing program such as Matlab from The Mathworks, which contains a built-in command for that operation. Entering

>> c = chol(R)

generates the following matrix:

From it you can read the following values, which represent the FIR filter coefficients: h(0) = 0.125, h(1) = -0.216, h(2) = 0.459, h(3) = -0.289 and h(4) = 0.801.

To check how good these filter coefficients are at generating desired ACF values, Matlab offers some effective tools. First, to check if the ACF values are positive definite, the easiest way is to compute R's eigenvalues. The command

>> e = eigen(r)

generates eigenvalues of 0.3750, 0.3521, 0.5971, 1.000 and 2.6758, which are all positive and thus make the ACF positive definite.

However, even though the positive-definite condition is true, Cholesky decomposition might contain slight errors due to numerical instabilities. To investigate this possibility, check matrix R's condition number. Recall from numerical analysis that a matrix's condition number indicates the amount of error possible when inverting it. Because Cholesky decomposition requires similar numerical procedures as found in matrix inversion, the condition number also indicates a possible error induced by this technique. Matlab offers several condition number routines, my favorite being

>> rcond(R)

which produces a value between 0 and 1 depending on whether the matrix is singular (0) or fully invertible (1). For this example, the condition number equals 0.123, which indicates that some error exists in the Cholesky decomposition (because 0.123 is closer to 0 than 1) but probably not enough to warrant finding new ACF values in my experiments.

To see how well the matrix under consideration works in generating colored noise, start with white noise (Fig 2a) and pass it through the FIR filter to arrive at the plot of colored noise in Fig 2b. It might be difficult to see any correlation in this signal, so to verify that the process has actually generated colored noise, take the FFT of the two noise signals with the resulting spectrum plots in Figs 2c and 2d.

Fig 2 -- To test this column's example FIR filter for how well it generates colored noise, consider the white noise in (a), which has the spectrum in (c) and ACF in (e). Compare it to the colored noise the filter generates (a), which has a spectrum that varies with frequency (d) and shows correlation between samples as shown in its ACF plot (f).

Note that the white noise has a relatively flat spectrum, while the colored noise shows a non-flat spectrum with more high-frequency than low-frequency components. In terms of autocorrelation, Fig 2e shows the white noise ACF, and as expected all values are near zero except at Lag 0. Finally, Fig 2f shows the ACF of the colored noise, which generates values similar to those that result if you take the ACF of the FIR coefficients and also shows some correlation among the first few lags. It thus appears that this example functions correctly. Also note that a relatively low condition number does induce some error in higher lag values, but it's for the most part negligible.

A promised solution

This article shows one useful application for the ACF. Along these lines I still owe you a solution for a problem I posed in a previous article, specifically, how to efficiently compute the ACF using Fourier transforms.

Suppose you want the ACF for a signal x(n) of length N. The process follows these steps:

1. Zero pad x(n), creating an extended sequence x_e (n) so that the new length and N_e = 2^K for some integer K. This step prepares the signal to handle the fact that the ACF is 2N - 1 long and attempts to minimize computational time by making sure that the procedure can use a radix-2 FFT.

2. X_e(k) = FFT(x_e(n)) where X_e(k) denotes the FFT of the extended sequence.

3. A(k) = X_e(k)X^*e(k) where the asterisk denotes a complex conjugate and A(k) contains power spectral density.

4. a(l) = IFFT(A(k)) where a(l) represents the ACF values.

5. These previous four steps equate time convolution to multiplication in frequency. The next step unravels the output into a more usable form:

6. R_UNBIASED = [ 1/ (N_e - abs(l)]R(l).

This procedure can also produce cross correlation by replacing Step 3 with A(k) = Xe(k)Y^*e(k).

I use this procedure often to compute an ACF, but a warning is in order. I've seen many FFT-based ACFs that ignore zero padding, so make sure in all convolution-based FFT procedures to include zero padding. Otherwise resulting convolutions are erroneous.

Matthew Perry is VP and general manager of Cirrus Logic's Embedded Processors Div (www.cirrus.com); Matt was with Motorola Inc when he wrote this article.

This article originally appeared in Personal Engineering & Instrumentation News, March 1994, pgs 65-71. Reprinted with permission of PEC Inc; all rights reserved.

		Click here to get your listing up.