Low leakage technique for determining power spectra of non-coherently sampled data

[0002] This invention relates generally to equipment and methods for testing electronic devices, and more particularly to a technique for extracting from an electronic test signal the amplitudes of the test signal's frequency components.

[0003] Test programs for automatic test systems commonly require a tester to measure the power spectrum of a signal sampled from a device under test (DUT). In a conventional testing scenario, an automatic test system generates a stimulus for exercising an input of a DUT and samples the output of the DUT as the DUT responds to the stimulus. Tester software computes the power spectrum of the sampled output signal by performing a Discrete Fourier Transform (DFT) on the samples acquired.

[0004] As is known, an error called “leakage” manifests itself in power spectra produced by a DFT whenever the sample clock is not “coherent” with the sampled signal. The sample clock is “coherent” if its frequency is a precise integer multiple of each frequencies present in the signal being sampled. Leakage is the mathematical consequence of performing a DFT on truncated frequencies—i.e., frequencies that do not complete a full cycle within the sample window. Leakage can be observed as an erroneous broadening of spectral lines, a creation of false peaks and troughs (lobes), and a general elevation of the power spectrum's noise floor.

[0005] Several methods have been devised to reduce leakage. One method is to increase sampling rate. In general, the higher the sampling rate, the smaller the amount of truncation in the frequency range of interest and the smaller the leakage error. Although effective, increasing sampling rate reduces leakage only in proportion to the magnitude of the increase. It also tends to greatly increase the cost of the sampling equipment used.

[0006] Another common technique for reducing leakage is to multiply the sampled data sequence by a windowing function. The windowing function has the effect of tapering the sampled data sequence around its endpoints, eliminating discontinuities that can give rise to leakage errors. Different windowing functions can be used, such as Blackman, Hanning, or Hamming windowing functions, each with its own particular characteristics. Windowing functions tend to diminish leakage errors distant from the peaks in a power spectrum, but also tend to create wider peaks. Thus, they have the effect of redistributing rather than completely eliminating leakage. Also, because windowing functions actually change the data on which a DFT is performed, they tend to slightly distort frequency spectra.

[0007] Yet another technique is to “resample” the waveform data at a sampling rate that is coherent with the frequencies of the signal being sampled. Resampling works by interpolating between actual points sampled at one rate to mathematically construct a series of points that appear to have been sampled at a different rate. Although resampling can be extremely effective for reducing leakage, it is computationally intensive and its accuracy can suffer from interpolation errors.

[0008] Still another technique for reducing leakage is to vary the rate of the sampling clock so that it precisely equals an integer multiple of every frequency found in the sampled signal. This technique is extremely effective, but requires expensive hardware. This approach is particularly expensive when a tester includes a large number of sample clocks, as is often the case.

[0009] Manufacturers of automatic test equipment (ATE or “testers”) commonly seek to improve their products by providing less costly solutions to conventional testing problems. Great benefits can be gleaned by increasing tester performance while decreasing tester cost. To this end, there is a strong need to provide an inexpensive technique for reducing leakage in the spectra of signals sampled by automatic test systems.

[0010] With the foregoing background in mind, it is an object of the invention to reduce leakage in sampled signals without requiring a significant increase in tester cost.

[0011] To achieve the foregoing object, as well as other objectives and advantages, a technique for analyzing the frequency content of a sampled waveform includes assembling a list of N frequencies expected to be found in the sampled waveform. The sampled waveform is assumed to conform to a waveform model that mathematically corresponds to a sum of N sinusoids. Each of the N sinusoids has unknown amplitude and phase, and a frequency that equals a different one of the N frequencies in the list of frequencies. The technique solves for the unknown amplitudes and/or phases that best fit the model to the sampled data.

[0012] According to one illustrative embodiment, the technique described above is also used when the frequencies of a sampled waveform are not known in advance. According to this variation, a Fourier Transform is computed on a sampled waveform to produce a rough power spectrum. Peaks in the rough power spectrum are identified, and their frequencies are compiled to form the list of N frequencies. Other factors may be considered in compiling the list of N frequencies, for example, known stimuli to a device from which the sampled data is acquired, and other attendant circumstances. The technique described above is then conducted on the resulting list of N frequencies, to determine the precise amplitudes and/or phases of each of the N sinusoids.

[0013] Additional objects, advantages, and novel features of the invention will become apparent from a consideration of the ensuing description and drawings, in which

[0014] FIG. 1 is a highly simplified schematic of an automatic test system stimulating the input of a device under test and sampling a signal from the output of the device under test;

[0015] FIG. 2 is a high level flowchart depicting a process according to the invention for determining the amplitudes of frequency components of a waveform sampled from the automatic test system of FIG. 1, when the frequencies of the frequency components are known in advance;

[0016] FIG. 3 is a flowchart depicting a process according to the invention for computing a power spectrum of a waveform sampled by the automatic test system of FIG. 1 when the frequency components are not known in advance; and

[0017] FIGS. 4A-4C are simulated power spectra comparing the performance of different techniques for reducing leakage in response to different simulated errors.

[0018] Technique

[0019] FIG. 1 is a highly simplified block diagram that shows a conventional configuration for testing a device under test (DUT) 120 with an automatic test system 110. The automatic test system includes a host computer 112. The host computer 112 includes tester software for running test programs. A test program controls the tester's hardware resources for testing a DUT. For example, a test program can control a frequency synthesizer 114 to apply a stimulus to the input of the DUT, and can control a digitizer 116 to sample a response produced at an output of the DUT.

[0020] The host computer generally stores sampled data acquired by the digitizer 116 in memory for analysis. The test program, or software routines available to the test program, operates on the sampled data stored to analyze its contents. Conventionally, the test program would direct the tester software to perform a Discrete Fourier Transform (DFT) on the sampled data. The test program would then test the results of the DFT.

[0021] FIG. 2 shows the general process according to the invention for sampling and analyzing test signals in an ATE environment. At step 210, the automatic test system 110 applies a stimulus to an input of the DUT 120.

[0022] At step 212, a list of N frequencies is accessed. The N frequencies represent frequencies of the sampled waveform about which amplitude and/or phase information is desired to be known. The number N can be any positive integer. Not all of the N frequencies need actually be present in the sampled waveform. In fact, this technique can be used to test for the presence or absence of any particular frequency component. Preferably, the list of N frequencies is known in advance and stored in the test program.

[0023] At step 214, the sampled waveform is computer modeled. The computer model representing the sampled waveform consists of N sinusoids summed together to approximate the actual sampled waveform. Each of the N sinusoids takes the form

[0024] A_ksin(ω_ki)+B_kcos(ω_ki),   (EQ 1)

[0025] where

[0026] “A_k” and “B_k” are unknown coefficients,

[0027] “k” is an index that ranges from 1 to N and represents one of the N frequency components,

[0028] ω_kcorresponds to the k^thfrequency component (specifically, ω_k=2πF_k, where F_kis the k^thfrequency), and

[0029] “i” is an index that identifies a particular sample and corresponds to time.

[0030] Although EQ1 appears to be the sum of two sinusoids, it mathematically equates to a single sinusoid having a frequency equal to ω_k/2π, an amplitude that equals {square root over (A²k+B²k)}, and a phase that equals the 2-argument arctangent of B_kand A_k.

[0031] Given that EQ1 represents each of the N sinusoids in the sampled waveform, the overall sampled waveform can be modeled by the expression $\begin{array}{cc}\sum _{k=1}^N\left(A_k\cos {\omega }_{k}i+B_k\sin {\omega }_{k}i\right).& \left(\mathrm{EQ}2\right)\end{array}$

[0032] At step 216 the waveform model of EQ2 is computer-processed to achieve a best fit between the model and the actual sampled waveform. The preferred embodiment employs a linear least-squares technique to fit the model to the data. In particular, step 216 attempts to minimize the following least-squares estimator: $\begin{array}{cc}{\sum _{i=0}^M\left[y_i-\sum _{k=1}^N\left(A_k\cos {\omega }_{k}i+B_k\sin {\omega }_{k}i\right)\right]}^2,& \left(\mathrm{EQ}3\right)\end{array}$

[0033] wherein y_iis the i^thsampled point of the sampled waveform, and i ranges from 0 to M, where M represents the total number of samples in the sampled waveform.

[0034] In order to minimize EQ3, the technique recognizes that a best fit is achieved when the partial derivative of EQ3 equals zero, when taken with respect to each value of A_kand B_k. Because there exists a value of A_kand B_kfor each of the N frequencies, taking the partial derivative of EQ3 with respect to each of these yields a system of 2N equations: $\begin{array}{cc}\sum _{i=0}^M{y}_i\cos {\omega }_{j}i=\sum _{k=1}^N\left(A_k\sum _{i=0}^M\cos {\omega }_{k}i\cos {\omega }_{j}i+B_k\sum _{i=0}^M\sin {\omega }_{k}i\cos {\omega }_{j}i\right)& \left(\mathrm{EQ}4\right)\\ \sum _{i=0}^M{y}_i\sin {\omega }_{j}i=\sum _{k=1}^N\left(A_k\sum _{i=0}^M\cos {\omega }_{k}i\sin {\omega }_{j}i+B_k\sum _{i=0}^M\sin {\omega }_{k}i\sin {\omega }_{j}i\right),& \left(\mathrm{EQ}5\right)\end{array}$

[0035] where EQ4 and EQ5 are each repeated N times as the index of “j” ranges from 1 to N.

[0036] Defining the following coefficients simplifies the description: $\begin{array}{cc}\mathrm{Let}{\mathrm{cc}}_{\mathrm{kj}}=\sum _{i=0}^M\cos {\omega }_{k}i\cos {\omega }_{j}i& \left(\mathrm{EQ}6\right)\\ \mathrm{Let}{\mathrm{sc}}_{\mathrm{kj}}=\sum _{i=0}^M\sin {\omega }_{k}i\cos {\omega }_{j}i& \left(\mathrm{EQ}7\right)\\ \mathrm{Let}{\mathrm{cs}}_{\mathrm{kj}}=\sum _{i=0}^M\cos {\omega }_{k}i\sin {\omega }_{j}i& \left(\mathrm{EQ}8\right)\\ \mathrm{Let}{\mathrm{ss}}_{\mathrm{kj}}=\sum _{i=0}^M\sin {\omega }_{k}i\sin {\omega }_{j}i& \left(\mathrm{EQ}9\right)\end{array}$

[0037] These coefficients can be mathematically simplified by recognizing that ${\mathrm{cc}}_{\mathrm{kj}}=\left(\sum _{i=0}^M\cos \left({\omega }_k+{\omega }_j\right)i+\sum _{i=0}^M\cos \left({\omega }_k-{\omega }_j\right)i\right)/2$ ${\mathrm{sc}}_{\mathrm{kj}}={\mathrm{cs}}_{\mathrm{kj}}=\left(\sum _{i=0}^M\sin \left({\omega }_k+{\omega }_j\right)i+\sum _{i=0}^M\sin \left({\omega }_k-{\omega }_j\right)i\right)/2$ ${\mathrm{ss}}_{\mathrm{kj}}=\left(\sum _{i=0}^M\cos \left({\omega }_k+{\omega }_j\right)i-\sum _{i=0}^M\cos \left({\omega }_k-{\omega }_j\right)i\right)/2,$

[0038] and that $\sum _{i=0}^M\cos \alpha i=\left(\cos \alpha \left(N-1\right)-\cos \alpha N-\cos \alpha +1\right)/2\left(1-\cos \alpha \right)\mathrm{and}$ $\sum _{i=0}^M\sin \alpha i=\left(\sin \alpha \left(N-1\right)-\sin \alpha N+\sin \alpha \right)/2\left(1-\cos \alpha \right),$

[0039] wherein α is any arbitrary value of ω.

[0040] Rewriting EQ4 and EQ5 using the coefficients defined in EQ6-EQ9 yields the following matrix: $\begin{array}{cc}\stackrel{V}{\left(\begin{array}{c}\sum _{i=0}^M{y}_i\cos {\omega }_{1}i\\ \sum _{i=0}^M{y}_i\sin {\omega }_{1}i\\ \sum _{i=0}^M{y}_i\cos {\omega }_{2}i\\ \sum _{i=0}^M{y}_i\sin {\omega }_{2}i\\ \sum _{i=0}^M{y}_i\cos {\omega }_{3}i\\ \sum _{i=0}^M{y}_i\sin {\omega }_{3}i\\ \cdots \\ \sum _{i=0}^M{y}_i\cos {\omega }_{N}i\\ \sum _{i=0}^M{y}_i\cos {\omega }_{N}i\end{array}\right)}=\stackrel{X}{\left(\begin{array}{ccccccc}{\mathrm{cc}}_{11}& {\mathrm{sc}}_{11}& {\mathrm{cc}}_{21}& {\mathrm{sc}}_{21}& \dots & {\mathrm{cc}}_{\mathrm{N1}}& {\mathrm{sc}}_{\mathrm{N1}}\\ {\mathrm{cs}}_{11}& {\mathrm{ss}}_{11}& {\mathrm{cs}}_{21}& {\mathrm{ss}}_{21}& \dots & {\mathrm{cs}}_{\mathrm{N1}}& {\mathrm{ss}}_{\mathrm{N1}}\\ {\mathrm{cc}}_{12}& {\mathrm{sc}}_{12}& {\mathrm{cc}}_{22}& {\mathrm{sc}}_{22}& \dots & {\mathrm{cc}}_{\mathrm{N2}}& {\mathrm{sc}}_{\mathrm{N2}}\\ {\mathrm{cs}}_{12}& {\mathrm{ss}}_{12}& {\mathrm{cs}}_{22}& {\mathrm{ss}}_{22}& \dots & {\mathrm{cs}}_{\mathrm{N2}}& {\mathrm{ss}}_{\mathrm{N2}}\\ {\mathrm{cc}}_{13}& {\mathrm{sc}}_{13}& {\mathrm{cc}}_{23}& {\mathrm{sc}}_{23}& \dots & {\mathrm{cc}}_{\mathrm{N3}}& {\mathrm{sc}}_{\mathrm{N3}}\\ {\mathrm{cs}}_{13}& {\mathrm{ss}}_{13}& {\mathrm{cs}}_{23}& {\mathrm{ss}}_{23}& \dots & {\mathrm{cs}}_{\mathrm{N3}}& {\mathrm{ss}}_{\mathrm{N3}}\\ & & & \dots & & & \\ {\mathrm{cc}}_{1N}& {\mathrm{sc}}_{1N}& {\mathrm{cc}}_{2N}& {\mathrm{sc}}_{2N}& \dots & {\mathrm{cc}}_{\mathrm{NN}}& {\mathrm{sc}}_{\mathrm{NN}}\\ {\mathrm{cs}}_{1N}& {\mathrm{ss}}_{1N}& {\mathrm{cs}}_{2N}& {\mathrm{ss}}_{2N}& \dots & {\mathrm{cs}}_{\mathrm{NN}}& {\mathrm{ss}}_{\mathrm{NN}}\end{array}\right)\stackrel{\mathrm{ab}}{\left(\begin{array}{c}{A}_1\\ B_1\\ A_2\\ B_2\\ A_3\\ B_3\\ \dots \\ A_N\\ B_N\end{array}\right)}}& \left(\mathrm{EQ}10\right)\end{array}$

[0041] One can solve EQ10 for each value of A_kand B_kin the ab vector by determining the inverse of the matrix X and multiplying it by the vector V on the left-hand side of EQ10.

[0042] Once A_kand B_kare known for each value of k from 1 to N, the amplitude of each k^thfrequency from the list of frequencies can be determined by computing {square root over (A²k+B²k)}. The value of each phase can be determined by computing the 2-argument ArcTangent of B_kand A_k.

[0043] The computational burdens imposed by solving EQ10 can be somewhat reduced by imposing a constraint that 2N=M+1 (the number of frequencies is half the number of samples). With this constraint in place, the vector V from EQ10 can be rewritten as follows: $\begin{array}{cc}\stackrel{V}{\left(\begin{array}{c}\sum _{i=0}^M{y}_i\cos {\omega }_{1}i\\ \sum _{i=0}^M{y}_i\sin {\omega }_{1}i\\ \sum _{i=0}^M{y}_i\cos {\omega }_{2}i\\ \sum _{i=0}^M{y}_i\sin {\omega }_{2}i\\ \sum _{i=0}^M{y}_i\cos {\omega }_{3}i\\ \sum _{i=0}^M{y}_i\sin {\omega }_{3}i\\ \cdots \\ \sum _{i=0}^M{y}_i\cos {\omega }_{N}i\\ \sum _{i=0}^M{y}_i\cos {\omega }_{N}i\end{array}\right)}=\stackrel{C}{\left(\begin{array}{cccccc}\cos 0{\omega }_1& \cos 1{\omega }_1& \cos 2{\omega }_1& & \dots & \cos M{\omega }_1\\ \sin 0{\omega }_1& \sin 1{\omega }_1& \sin 2{\omega }_1& & \dots & \sin M{\omega }_1\\ \cos 0{\omega }_2& \cos 1{\omega }_2& \cos 2{\omega }_2& & \dots & \cos M{\omega }_2\\ \sin 0{\omega }_2& \sin 1{\omega }_2& \sin 2{\omega }_2& & \dots & \sin M{\omega }_2\\ \cos 0{\omega }_3& \cos 1{\omega }_3& \cos 2{\omega }_3& & \dots & \cos M{\omega }_3\\ \sin 0{\omega }_3& \sin 1{\omega }_3& \sin 2{\omega }_3& & \dots & \sin M{\omega }_3\\ & & & \dots & & \\ \cos 0{\omega }_N& \cos 1{\omega }_N& \cos 2{\omega }_N& & \dots & \cos M{\omega }_N\\ \sin 0{\omega }_N& \sin 1{\omega }_N& \sin 2{\omega }_N& & \dots & \sin M{\omega }_N\end{array}\right)}\stackrel{y}{\left(\begin{array}{c}{y}_0\\ y_1\\ y_2\\ y_3\\ y_4\\ y_5\\ \cdots \\ y_{M-1}\\ y_M\end{array}\right)}& \left(\mathrm{EQ}11\right)\end{array}$

[0044] Since imposing the above constraint forces C and X⁻¹to be square matrices of the same rank, EQ10 and EQ11 can be combined to form

[0045] ab=(X⁻¹C)y  (EQ12)

[0046] Thus, it is possible to solve for ab without having to compute V.

[0047] The following method is suggested for solving for ab:

[0048] First, construct C quickly using the following recurrence relations:

[0049] cos(a+1)ω_k=2cos ω_kcos aω_k−cos(a−1)ω_k

[0050] sin(a+1)ω_k=2cos ω_ksin aω_k−sin(a−1)ω_k

[0051] Next, compute X⁻¹via L-U decomposition

[0052] Apply X⁻¹to C

[0053] Compute ab by multiplying X−1C by y

[0054] Once X⁻¹C has been constructed, computing ab requires roughly N²multiply-accumulates.

[0055] Preferably, the technique described herein is implemented as a function in a software library. The function preferably receives an input array that stores the list of frequencies and a pointer to the sampled data. The function preferably returns an array containing the values of A_kand B_k, from which amplitude and phase can be computed. Alternatively, the function returns amplitude and phase directly. The software library preferably resides on an automatic test system where it is accessible to test programs running on the test system.

[0056] Example

[0057] FIGS. 4A-4C illustrate a simulated prediction of the performance of the best-fit technique according to the invention versus other techniques for reducing leakage. A listing of the test code used to produce the data for these charts is provided at the end of the specification. Each of the three graphs of FIGS. 4A-4C compares the power spectrum of a single tone obtained under four different conditions:

[0058] 1. a Fast Fourier Transform (FFT) of uncorrected data (i.e., acquired with a rectangular window, labeled “uncorrected”)

[0059] 2. an FFT of data shaped by a Hanning window (“windowed”),

[0060] 3. an FFT of resampled data (“resampled”), and

[0061] 4. the best-fit technique described herein (“least squares”).

[0062] The horizontal axis of each graph corresponds to frequency—specifically, frequency bins 0-63. To allow for direct comparisons with methods that employ an FFT, the best fit technique was run with N=64, where each of the N frequencies corresponds to an FFT bin. The vertical axis corresponds to amplitude in db.

[0063] Each graph shows a single tone that is slightly askew from the center of the 9^thfrequency bin, i.e., which is not coherently sampled. In particular ω_k=2π(k−1)(1+ε), where ε=10⁻⁶for FIG. 4A, 10⁻⁹for FIG. 4B, and 10⁻¹²for FIG. 4C. Each of these graphs reveals the strength of the best-fit technique. As compared with the other techniques, the best-fit technique maintains an exceedingly narrow peak without elevated areas (“skirts”) surrounding the peak.

[0064] Advantages

[0065] The disclosed technique offers many advantages over the conventional DFT, particularly in the context of automatic test equipment. Using the technique described, highly accurate spectral analysis can be performed using relatively inexpensive tester electronics. The sample clock need not be coherent with the frequencies to be measured, and spectral leakage is substantially eliminated. As compared with a conventional DFT, in which frequencies are sorted into bins that have finite widths, the technique disclosed does not employ frequency bins but rather discrete frequencies. Thus, the techniques disclosed enable exceedingly closely spaced frequencies to be resolved in a way that would be impossible using a conventional DFT.

[0066] The technique is also scalable, because its computation time varies as a function of N—the number of frequencies to be analyzed. Thus, the technique can be performed relatively quickly if only a small number of frequencies is being analyzed. In addition, because the technique employs a best-fit algorithm, it can be used to determine the amplitude and phase of frequency components that do not fully complete a cycle within the sample window. The requirement that frequencies be prescribed in advance does not generally operate as a disadvantage in automatic test equipment, where frequencies produced by a DUT are largely known in advance, and where the tester generates the stimuli for driving the DUT.

[0067] Although X⁻¹C requires significant time to compute, it need not be recomputed each time a waveform is analyzed. As long as the frequency list and the number of samples M remains constant, newly acquired sampled can be analyzed by retrieving a stored copy of X⁻¹C and multiplying it by y. Libraries can be provided that include a variety of different X⁻¹C combinations for different frequencies and numbers of samples. Users can select from among the combinations to provide rapid analysis of waveforms.

[0068] Alternatives

[0069] Having described one embodiment, numerous alternative embodiments or variations can be made.

[0070] For example, as described above, the technique for analyzing sampled waveforms is used in the context of automatic test equipment. However, the technique can be more generally applied to any sampled data for which analysis of frequency information is desired. Although the technique is illustrated with respect to the particular test scenario of FIG. 1, it is not limited to any particular test scenario.

[0071] The embodiment disclosed herein specifies the use of linear least squares for obtaining a best fit between the waveform model (EQ2) and the actual sampled data. However, other best-fit techniques can be used, such as Cauchy-Lorentz distributions and techniques that attempt to minimize the absolute value of difference between the model and the sampled data. Therefore, the invention is not limited to the use of least-squares.

[0072] In addition, the description provides that software be used for manipulating matrices and performing the requisite computations. Alternatively, specialized hardware circuits or processors could be provided to more efficiently perform these functions.

[0073] As described above, the technique requires that the list of N frequencies be known in advance. This requirement can be avoided, however, by performing a DFT on the sampled data and examining the results. This variation is illustrated in FIG. 3. At step 310, a DFT is performed. At step 312 peaks are identified in the power spectrum that results from the DFT. The frequencies corresponding to the peaks are then be appended to the list of frequencies for performing a more accurate analysis. Steps 316 and 318 proceed as already described: by assuming that the sampled data conforms to a model and obtaining a best fit between the samples and the model. This technique need not operate in a vacuum. It can also take attendant circumstances into account, such as the frequencies of stimuli applied to the DUT, harmonics of those frequencies, and known characteristics of the DUT.

[0074] The technique described above assumes that waveforms are sampled at a uniform rate. However, according to an alternative embodiment, waveforms can be sampled at non-uniform rates. In particular, replacing the discrete index “i” in the equations and matrices above with the term “t_i” (i.e., the actual sample times) allows for arbitrarily non-uniform sampling. With non-uniform sampling, the computational simplifications that follow EQ9 may not be used; however, the change is transparent to the remainder of the technique described.

[0075] Each of these alternatives and variations, as well as others, has been contemplated by the inventors and is intended to fall within the scope of the instant invention. It should be understood, therefore, that the foregoing description is by way of example, and the invention should be limited only by the spirit and scope of the appended claims.

[0076] Computer Listing

[0077] A software listing of the test code used to produce the data for these charts is provided below:

[Inline image omitted. See PDF.]



[Inline image omitted. See PDF.]



[Inline image omitted. See PDF.]



[Inline image omitted. See PDF.]



[Inline image omitted. See PDF.]