Bandlimited Interpolation
Robert Gilchrist Huenemann, M.S.E.E.
January 6, 2006
Introduction
Have you seen ads for compact disc players? They are full of references to oversampling. What does this mean? Do they somehow give you something for nothing?
Of course they do not. Modern CD players all use band-limited interpolators. Band limited interpolation can be performed in either the frequency or the time domain, or as part of the process of conversion from one domain to another. I hope that this article will help you to understand these useful digital filters.
This article is not intended to be an exhaustive review of digital filtering, or of Fourier analysis. If you have some background in these areas, I hope that you will be able to follow this discussion more easily. If you exposure to these subjects is limited, I hope that I can stimulate your curiosity.
Digital filters are applied in sampled data systems. In these systems, an electrical signal (which may represent a physical phenomenon such as a sound wave) has been measured precisely at a series of equally spaced times. The digital filter operates on a series of measured values, or numbers. Let’s assume that the measurement itself is perfect – there is no jitter and no quantization error.
We also assume that the signal is band-limited. This means that there are no frequencies present which are greater than half the sampling rate. This of course is the well-known Nyquist criterion, as it applies to baseband signals. The more general case applies to passband signals, in which case the number of complex samples equals the sample rate. These cases can be approximated by using suitable analog filters ahead of the measurement process.
If the measurement is quantized, or if the Nyquist criterion is not strictly met, the signal will contain errors. These errors have been studied extensively, but we will ignore them for purposes of this discussion except to note that they can be made arbitrarily small, for a price.
We have a series of samples which perfectly represent a waveform at equally spaced times, and the waveform was band-limited to less than half the sample rate. This means that we have all the information needed to compute the spectrum of the waveform. But if we know the spectrum, we should be able to reconstruct the signal for every point in time – not just at the original sample points.
It turns out that we can do exactly that, and rather easily. In fact, we can do it in a couple of ways. Let’s consider each.
Frequency Domain Interpolation
Time Domain Frequency Domain
We can perform interpolation in the frequency domain. This approach is described in Figure 1. We compute the spectrum with a discrete Fourier transform, presumably by using one of the fast algorithms. If we have N samples, we then have N pieces of spectral information.
This relationship bears further comment. There are always N independent pieces of spectral information that can be computed from N time domain samples. The independent spectral lines fall at integer submultiples of the sampling rate. The so called fast Fourier transforms are carefully designed to give only those N pieces of information. Any additional computation would be redundant. But when we interpolate, we are deliberately computing redundant information. Another way to view this is to recognize that the fast algorithms have been deliberately designed to do no interpolation.
In any case, once we have computed the real and imaginary components of the spectrum, we can do something called “zero-extending”. We take the original arrays, containing the reals and imaginaries, and place them at the ends of arrays which are M times as long as the originals. The middles of these arrays are filled with zero values. Figure 1 illustrates this process.
We now perform an inverse discrete Fourier transform (IDFT) on the new spectrum, which is M times as long as the original. This gives us M * N samples in the time domain. The samples all lie on the original time span, at equally spaced intervals. Thus we now have M time samples for each original sample.
These samples must lie exactly on the original unsampled waveform. If they did not, the interpolated waveform would not have the same spectrum as the original. Either the spectrum below the Nyquist frequency would be different, or there would be spectral components above the Nyquist frequency. But we constructed the new samples from a spectrum which was identical to the original below the Nyquist frequency, and zero above. Thus, the new samples fall exactly on the original unsampled waveform.
Band-limited interpolation is a name which describes this process nicely. Note that any other interpolation scheme, no matter how sophisticated, will generate points which lie slightly off of the original unsampled waveform, which means the spectrum of the interpolated waveform will be slightly different. If you find yourself considering parabolic curve fitting, spline functions, or the like, you probably will find that band-limited interpolation will do a better job with less computation.
Because we zero-extended the spectrum, there are many steps in the inverse Fourier transform which consist of multiplication by zero. The computation can be shortened by eliminating these steps. This process has become known as “pruning”.
Time Domain Interpolation
Your CD player isn’t doing Fourier transforms. Band-limited interpolation can be implemented in the time domain by using a class of FIR (finite impulse response) filters known as polyphase filters. Typical FIR filters have symmetrical coefficients. Polyphase filters are filter banks. Each individual filter in the bank has asymmetrical coefficients.
Please refer to the five parts of Figure 2 for this discussion. The process of designing a polyphase filter can be visualized by imagining that a series of zero values are inserted between each measured input. This combined series of values is fed into an FIR filter with the desired frequency response. This is known as zero-packing in the time domain, and it is the equivalent of zero-extending the spectrum.
The resultant filter computation contains many multiplications by zero. If these steps are eliminated or pruned, the filter becomes a set of shorter filters with asymmetrical coefficients which are a subset of the original symmetrical coefficients. Each filter computes one output value for every input value. When the outputs are arranged in proper order they correspond exactly to the outputs from a frequency domain interpolation.
After band-limited interpolation, the results are fed to an analog-to-digital converter and an analog low pass filter. The signal coming out of the converter looks like a set of stairs. If the steps are smaller and closer together, the analog filter will do a better job of converting them to one continuous waveform. This is why band-limited interpolators are useful in CD players. They make it possible to use a simple analog reconstruction filter with less phase shift to accomplish a given level of harmonic distortion in the output signal.
If course, all these harmonics lie above your range of hearing so don’t be too impressed by claims of eight times oversampling. On the other hand, certain measurement instruments use band-limited interpolators to achieve a required level of accuracy.
Spectral Analysis
I want to briefly mention one method which is not an interpolation technique. One of the challenges of spectral analysis is to achieve adequate resolution of each spectral component. The Fourier transform gives good resolution if we have many samples in the time domain. Normally this requires a long Fourier transform algorithm. One way to overcome this limitation is to transform a small number of time domain samples, save the results, and repeat the operation on the next group of samples.
The output of the Fourier transform can be viewed as a bank of filter outputs, and the data sequence out of each filter represents the time history of a certain range of frequencies. Therefore, the sequence of outputs from any given filter can again be transformed to the frequency domain. Arbitrarily large sets of samples can be transformed in this way, and the resolution increases with the total size of the sample set. This technique has been called the ZOOM transform.
Now back to interpolation. What can we do if we want to increase the apparent spectral resolution from a small number of time domain samples? We zero-extend the data in the time domain, by arbitrarily adding zeroes after the measured data. After Fourier transformation, the zero-extended data yields a spectrum with interpolated points between those for the uninterpolated spectrum. If no window is used, each spectral line will have interpolated sidelobes with a sin(x)/x amplitude. The sidelobes can be controlled by windowing the time domain data. More information about windows can be found in the first two references.
The Chirp z-Transform
A more general way to get arbitrary apparent spectral resolution is to use the input data to compute the discrete Fourier Transform at a single basis frequency. This involves multiplying the data by sine and cosine waveforms with the desired basis frequency, and forming a running sum of the products, to get the amplitude and phase of the spectrum at the basis frequency.
Now imagine that we vary the basis frequency in arbitrarily small steps. We can now compute the spectrum with any desired resolution. Again, we are not creating any new information. We are just providing a smoothly interpolated line between the N independent data points.
The method just described can be used to compute an arbitrarily large number of points of the spectrum, over an arbitrary frequency range, from a small number of time domain samples. However, the process is slow. Lawrence Rabiner, of Bell Labs, recognized that the fast Fourier transform algorithms could be used to speed up this process.
To understand this last statement, remember that the fast Fourier transform algorithms have been carefully designed to return information about the independent spectral lines which fall at submultiples of the sampling rate. They start with N points in the time domain, and compute exactly N complex spectral products.
Rabiner started with a deceptively simple relationship first noted by Bluestein, and from it developed the chirp z-transform. This transform has been published as part of the IEEE digital signal processing library.
The chirp z-transform takes its name from the classical z-transform, and from the fact that one of the functions involved has a linearly changing frequency. Radar sets which use linear frequency sweeps are known as chirp radars, by analogy with a bird call.
The chirp z-transform starts with a setup process which takes into consideration the number of samples, input time range and output frequency range. Based on these parameters, a complex array is computed which contains the linear frequency sweep, or chirp. The setup computation includes a Fourier transformation.
The measured time domain data is next zero-extended, based on the number of spectral points desired. The array is Fourier transformed, multiplied by the chirp array, and Fourier transformed a second time. The resulting array contains the desired portion of the spectrum.
The Inverse Chirp z-Transform
If we have a series of measurements in the frequency domain (such as those from a network analyzer) and we want to convert them to an arbitrary range in the time domain, the inverse chirp z-transform can be used. This is a relatively simple modification to the algorithm, and it is used in time domain options to network analyzers. An inverse chirp z-transform can be realized by passing a negative value for one of the setup parameter values to the standard chirp z-transform algorithm.
Wavelets
It is my impression that wavelets may be related to band limited interpolators.
Fast Transform Algorithms
The most familiar fast transform algorithm is the power of two, or Cooley-Tukey, transform.
The prime factor, or Winograd, algorithm is approximately twice as fast as the Cooley-Tukey algorithm.
If the number of points to be transformed is small (approximately 32 or less), the discrete Fourier transform is faster than any of the ‘fast’ algorithms. And if the application requires the computation of only a small number of spectral lines as opposed to the full spectrum, the discrete Fourier transform may also be faster than a ‘fast’ algorithm.
Testing of DSP algorithms
There is a simple, obvious tool that I have found to be of enormous value when developing and implementing various kinds of DSP algorithms. Simply use a very small number of data points for the test case. It is common to test algorithms with special data sets such as impulse functions or sinusoids. Whatever the test data may be, problems that would be invisible with thousands of data points become painfully obvious with a small number such as eight or 16. The importance of splitting end point data, and the question of whether to use n versus n-1 in an algorithm, can easily be understood if the number of data points is deliberately made quite small.
Conclusion
The band-limited interpolators are elegant algorithms. Using
them has been a challenge, and a source of satisfaction. I am especially
grateful to fred harris of
References
Elliott, HANDBOOK OF DIGITAL SIGNAL PROCESSING ENGINEERING APPLICATIONS, Academic Press, 1987. Note in particular chapter 3 (Multirate FIR filters for Interpolating and Desampling) and chapter 8 (Time Domain Signal Processing with the DFT), by frederic j. harris.
f.j. harris, “On the use of windows for harmonic analysis with the discrete Fourier transform,” Proc. IEEE, vol. 66, pp. 51-83, Jan. 1978.
Rabiner and Gold, THEORY AND APPLICATION OF DIGITAL SIGNAL PROCESSING, Prentice-Hall, 1975. Note the discussion of Bluestein’s algorithm and the chirp z-transform on pages 392-399.
Rabiner, Schafer and Rader, “The Chirp z-Transform Algorithm and its Application,” Bell Sys. Tech. J., vol. 48, no. 5, pp. 1249-1292, May-June 1969.
IEEE Digital Signal Processing Library. This library contains a FORTRAN implementation of the chirp z-transform. It is available in hard copy and machine readable form.