The general topic for this article is sampling, which most readers likely feel they're familiar with, but I want to cover it from a different perspective. In particular, the interpolation formula provides all the information of the Nyquist sampling theorem, but its form allows some new interpretations and insights.

Something for nothing?

First, though, let's straighten out some terminology. Although sampling is the basic tenant for DSP, it seems strange that the "D" stands for digital as opposed to discrete. In the DSP world we quantize a continuous signal using an A/D converter, which also samples the signal's index (normally time). Thus, most of us actually work on DVDTSP (digital-value discrete-time signal processing), but that acronym just doesn't sound as cool as DSP. A key distinction between the two is that while you can't reverse the digitization process, under certain circumstances you can reverse sampling with the aid of the interpolation formula.

Let's also review some fundamentals, starting with the Nyquist sampling theorem (NST), which tells how often you should sample a signal. As most everyone knows, it states that if you sample a continuous signal at a frequency at least twice as great as the highest frequency it contains (assuming a dc component exists), the digital representation contains exactly the same information as the continuous signal.

Amazingly, that theorem says you can represent an infinite number of points with only a discrete number of points. The skeptic in you should immediately question such a claim -- after all, you can't get something for nothing. The kicker is that the signal must be bandlimited, and so in a formal sense you can't apply it directly to real-world signals because nature always generates noise that's not bandlimited.

The good news is that the practical side of us overtakes the theoretician because we realize that most signals are approximately bandlimited. Thus you either filter a signal to force bandlimitedness, or you intentionally forget that a few high frequencies exist. Either way, some error creeps into the sampling process, but this error is practically nonexistent. In practice we usually choose a digitizing rate sufficiently high to more than account for the bandwidth of interest and thus allow for good representation. For instance, compact audio disks sample at 44.1 kHz, more than twice the 20-kHz frequencies for which our ears can detect sound (in fact, most people can only hear to about 16 kHz).

To get a mathematical feel for the NST's validity, consider the signal x(t) = Asin(2πft) where f is unknown. You know the signal is bandlimited to f, but without the NST, how would you determine how often to sample it yet retain all information? Sampling x(t) twice during one period at two arbitrary points in time generates the following equations:

x(t₁) = Asin(2πft₁)

x(t₂) = Asin(2πft₂).

You can solve them for the two unknowns, A and f, thus showing it takes only two samples per period to resolve both amplitude and frequency -- feedback that proves the NST actually makes sense.

Fig 1 -- When you take the FFT of a squarewave such as that in (a), the magnitude plot shows frequencies that go to infinity. Also note in (a) that you can completely define the waveform with just two samples.

An interesting aside is that the NST doesn't always produce a minimal set of samples. Consider a continuous squarewave (Fig 1a) and its Fourier transform (Fig 1b), which shows that the signal contains components for all frequencies. According to the NST, to retain all signal information the sampling rate must be infinite. However, by sampling the square-wave at the two points in Fig 1a you can capture it with a digitization rate of two points per period. This methodology, though, requires you to make an assumption about the signal's shape, as opposed to the NST's ability to handle arbitrary signals.

Is f_s = 2f_max sufficient?

A related question concerns how to take discrete sampled signals and reconstruct a continuous signal. Again look to the NST. Nyquist showed that

which is typically called the reconstruction formula; and because it interpolates values between samples, it's also known as the interpolation formula.

It's worthwhile studying this formula for a moment. First, you can rewrite it as x(t) = x[n] ⊗ h(n) where h(n) = sinc{π(t - nT) / T}. Thus, the continuous values are just the convolution of the discrete values with a sinc function. A better appreciation of the interpolation formula in turn forces you to understand the sinc function, whose generic form is sinc(πx / T). Fig 2 shows a sinc function for n = 0, and its shape gives an idea of how the reconstruction formula functions.

Fig 2 -- The shape of the sinc function has an important role in the interpolation function.

To see how the sinc function comes into play, I generally examine the interpolation formula and determine what each x[n] is multiplied by before summing, giving an indication of how much each discrete value contributes to the overall x(t). To do so, realize that the generic sinc function equals 1 when its argument equals 0, and it equals 0 for integer multiples of the sampling rate (note that the sinc function's value lies between 0 and 1 for all other arguments). As an example, suppose you want to know x(t) for t = 0. The interpolation formula indicates that x(t = 0) = x[0] due to the sinc function's values. For other values of t you can see that each x[n] contributes to x(t) because the sinc function is nonzero at those points.

The interpolation formula also helps answer a perennial DSP question concerning the sampling rate (f_s): must f_s > 2f_max, or is f_s = 2f_max sufficient? You'll find one of these two forms in every DSP book, but which is correct? To dig into the answer, suppose you sample a sinewave twice per period (thus adhering to NST, where f_s = 2f _max). Also suppose the plausible case where you sample the sinewave exactly where it crosses the time axis. Now, if you try and find an intermediate value, say x(1.5), the interpolation formula gives an answer of zero because all x[n] equal zero. This result is clearly incorrect, so we've shown that f_s > 2f _max to guarantee exact reconstruction.

As you can see, the interpolation formula is a convenient way to understand the NST. Most people stress the theoretical aspects, as I've done here. However, a practical side does exist. Given a sequence of discretely sampled values, how do you determine the exact value at some arbitrary time t? The interpolation formula is the place to get that value, as long as you sampled according to Nyquist.

The exact reconstruction works only if lots of data exist (practically an infinite amount so that the formula's sinc function dies out to near zero). But what happens for a finite-length sequence, which is what you mostly find in practice? Assume you sample x(t) according to NST and store x[0],...,x[M]; assuming that T = 1 sec, that fact means you sampled it from t = 0 to M sec. What does x(1.5) equal? As just established, you can find the answer with the interpolation formula, but you require x[n] for practically all n.

The next question then becomes what values to assign to x[n] for n < 0 and n > M. Because you know nothing about these values and don't care what they are, a first-pass guess might be zero, which means you're probably assuming that x(t) = 0 for t < 0 and t > MT. The problem is that you've just violated the NST. Why? The answer comes from the Fourier transform property that finite-length continuous signals have infinite-length spectra (and the opposite is also true). Thus, in this case, x(t) has a highest frequency of infinity, which means you can't sample x(t) fast enough to adhere to the NST.

How can you get around this knotty problem? Well, basically you can't because you just don't know what lies in the future or what has happened in the past. However, for interpolation problems I've found that assuming x(t) = 0 outside the region of interest is fine for most applications and doesn't lead to many problems.

Filling between samples

While thinking about how the interpolation relates to the real world, consider another example that shows the practical side of this topic. Suppose you encounter a digitized signal x[n]; how can you double the number of samples between any two time values? Your first answer should be to use the interpolation formula, and while that response is theoretically correct, it's impractical because the computational load is large.

A quick alternative to the interpolation formula that produces identical results uses the FFT. Suppose you have x[n] for n = 0,...,15 (Fig 3a) and want to double the number of samples to 32.

Fig 3 -- Given one cycle of a sinewave that you sample 16 times (a), padding it with an extra 16 zeros and taking the inverse FFT of the 32-sample signal produces the signal represented by the Xs in (b). Note how the ones at the extreme edges lack the accuracy of those in the center.

First find the FFT, X(k), and pad it with 16 zeros, yielding X_e [k]. Taking the inverse FFT of X _e [k] produces the desired 32 values of x[n] as in Fig 3b, whose results are interesting. Notice that this procedure correctly interpolates points that already existed. As for the new points, a closer examination shows that those close to the center of the waveform (n = 7.5 sec) are the most accurate, while those at the outer edges lose some accuracy, an effect especially pronounced near n = 16 sec. The reason for this discrepancy from the true values arises because the sinc function that's part of the interpolation formula has the most effect near the center of the waveform, which becomes apparent if you imagine the plot of Fig 2 overlaid on Fig 3b.

This example also employs a little-known DSP trick: frequency domain zero-padding, where you append zeros on the end -- but be careful what you mean by "end." To appreciate the complexities of zero padding in the frequency domain, it's easiest to recognize what each frequency value represents. The first value is the dc term, X(0). Next come the positive frequencies X(1), X(2),..., X(?). The last value, X(?), is unclear because the sequence must also include negative frequencies, so also add X(-1), X(-2),..., X(??). Putting all this information together in the correct order produces X(k) = [ X(0), X(1) ... X(?) X(-??)...X(-2), X(-1) ].

The problem then becomes determining the values for the two indices "?" and "??". The entire sequence is N samples long, but X(0) isn't replicated, so for an even value of N (true for many practical situations) the number of positive frequency components doesn't equal the number of negative frequency components, so both X(?) and X(??) can't be in the sequence. For now assume this single value is denoted X(???). Then the Fourier transform has the form of X(k) = [ X(0) X(1)...X((N/2)-1) X(???) X(-((N/2)-1))...X(-2) X(-1)].

As for the middle point X(???), should you use X(N/2) or X(-N/2)? FFT theory says X(N/2), and notice that the negative component of this frequency value doesn't occur in the sequence. However, had you chosen an odd-length sequence, then X(-N/ 2) would occur. Based on this analysis, you now know that padding in the frequency domain means inserting the zeros between X(N/2) and X(-((N/2)-1)).

Naturally discrete events

Before concluding this discussion of sampling, let's review one final aspect: what happens if you encounter samples but have no information about the original digitization rate or the continuous signal's bandwidth? Even worse, what if the data are naturally discrete? These scenarios actually arise often in many fields. For example, financial analysts deal with discrete data such as the Dow Jones Industrial Average, currency values and economic indicators, and these data sequences aren't based on continuous signals. However, I've seen much DSP-like analysis on financial data, and problems sometimes crop up due to analysts forgetting about the NST, which means that high-frequency components might be highly aliased, which leads to misinterpreting data. For those who enjoy playing with DSP algorithms, analysis of such discrete time series offers unique challenges.

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, Sept 1994, pgs 63-67. Reprinted with permission of PEC Inc; all rights reserved.

		Click here to get your listing up.