Method and apparatus for waveform reconstruction for sampled data system
BACKGROUND OF THE INVENTION 
     This invention relates to sampled data signal analysis and in particular to methods and apparatus for sampling and reconstructing the wave shape of substantially nontransient data wherein processing of samples of signals occurs at a rate substantially less than the Nyquist rate of the highest frequency components of the sampled signal but higher than the Nyquist rate of the lowest frequency components. A related, but distinguishable field is that of direct measurement of signals. In that field, measurement techniques are constrained by the absolute speed limitations of the components in the system. In the present field, techniques have been developed for making measurements wherein the signal analysis components of the measurement system operate at frequencies or speeds substantially less than the frequency components of the measured signal. Such techniques are primarily employed for measurements of nontransient or slowly-changing signal components of high-frequency (e.g. microwave) steady-state signals. 
     In connection with the need to measure slowly-changing signal components in high-frequency, steady-state signals, such as signals above 1 GHz, there is a need to examine and record fast rise-time (or fall-time) signals. In the past, it has been necessary to provide an input signal trigger. Such a scheme may require circuitry capable of operation at speeds comparable to or better than the signal to be analyzed. Since this is not always practical, because of limitations imposed by the available analysis tools, an alternative is clearly needed to the traditional hardware-triggered sampling data capture system. 
     More specifically, what is desired is an ability to record for analysis fast rise-time signals with a conventional signal digitizer (which operates at a substantially lower speed than the signal to be analyzed). For example, a 20 million sample per second digitizer sampling a signal with 20 GHz components implies analysis by circuit components operating at 20 MHz upon a 1/1000th harmonic with a best resolution of 50 nanoseconds between each sample. In theory, a signal with a 50 nanosecond-range rise time waveform could be recorded for analysis. This is an extremely difficult task with a conventional digitizer employing state-of-the-art components and circuitry. In reality, a digitizer with a resolution of 50 nanoseconds can only resolve for identification purposes a signal having frequency components with a rise time of greater than about 100 nanoseconds. At best, it is possible to determine that a signal has a rise time of 50 nanoseconds, but it is not possible to determine important characteristics, including exactly when the event started or the shape of the rising signal between samples. 
     Random event signals are particularly difficult to analyze. Depending on the point of random sampling, one might be able to record the occurrence of an event at a mid-point between a series of samples representing stable states of a signal. If so, one would know only that an input signal experienced a transition during an interval of somewhat less than 100 nanoseconds because the signal would require a minimum period of two samples to move between stable states. Little more would be known about the signal. 
     It is nevertheless desirable to be able to examine in detail an event represented by an input signal whose resolution is on the order of only twenty picoseconds. 
     It is believed that there is no prior art directly relevant to the present invention. However, in the past, an instrument manufactured by Hewlett-Packard Company under the model designation 54100 has used random repetitive sampling whereby triggered samples of a repetitive waveform are taken every 25 nanoseconds, and a hardware trigger interpolator is employed to determine where each sample has occurred with respect to the trigger, i.e., whether it occurred before the trigger or after the trigger and by how much before or after the trigger it occurred. Based on where it occurred with respect to the trigger, a dot or sample point representing the position of the sample in a magnitude versus time relationship is directed to an output device such as a sampling oscilloscope screen or a printer, wherein the dot is stored and/or displayed as an element of a reconstructed or synthesized waveform. Due to the inherent shortcomings on the accuracy of triggering, there is a limitation on the accuracy of such a synthesized waveform at high frequencies. 
     SUMMARY OF THE INVENTION 
     A method according to the invention comprises the steps of registering a time record of a large plurality of samples of a repetitive signal, determining the frequency of the repetitive signal very precisely, sorting based on frequency and the sequence of the samples, and presenting the time record on a reconstructed common time base (one period) in order to represent the repetitive signal with a much finer resolution than would be possible for a sampling technique based on sampling a single period. The invention does not employ a circuit for triggering each sampling operation. 
     The invention herein described does not require triggering based on a hardware circuit. Rather, the invention relies on software reconstruction of a repetitive waveform which has been randomly sampled and then analyzed to calculate the precise frequency of the fundamental component and then from that the beginning of each cycle and then the phase of all other sampled points. This technique is limited to the analysis of a single repetitive waveform. Its reliance on so-called software triggering means eliminates the need for a high-speed triggering circuit in order to build an effective signal sampling instrument. 
     According to a specific embodiment of the invention, a string of samples of a repetitive, input signal with high frequency components is registered relative to an arbitrary time base for a recording time which extends over many periods of the repetitive signal. The samples are captured without triggering with relatively low time resolution to determine an approximate waveform from the low resolution samples, then digital signal processing techniques in the form of a fast Fourier transform (FFT) are applied to a reconstructed time record of the input signal to obtain an accurate frequency for the signal, and finally the sampled waveform is reconstructed by overlaying samples with reference to a common time or phase reference. 
     In accordance with a specific embodiment of the invention, an FFT is employed to determine the fundamental frequency of the signal component precisely. Other methods can be employed as alternatives, including for example a discrete Fourier transform (DFT) of zero crossings, a chirp transform or other time to frequency transforms. Further processing, such as bin interpolation as explained herein, is employed to improve frequency resolution. For example, if the shape of the sampling window in the time domain is known and a particular result in the frequency domain is produced, then it is possible to use interpolation to find the frequency of the input to sub-bin resolution. Thereby it is possible to overlay samples taken over multiple periods to resolve rise times that are extremely fast, e.g., as short as ten nanoseconds, and to take advantage of the entire bandwidth of the sampling system. 
     According to the invention, therefore, the function of the traditional trigger circuit is synthesized so that random repetitive sampling of any signal can be accomplished without reference to a triggering signal. Using sample values representing hundreds of samples of signal edges, it is possible to determine the shape of the waveform as well as its frequency to an accuracy in excess of known trigger-based sampling techniques. 
     The invention will be better understood by reference to the following detailed description taken in conjunction with the accompanying drawings. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a functional block diagram of the method according to the invention. 
     FIG. 2 is a block diagram of an apparatus according to the invention. 
     FIG. 3A is a diagram illustrating a frequency correction function useful in bin interpolation. 
     FIG. 3B is a diagram illustrating a magnitude correction function useful in bin interpolation. 
     FIG. 4 is a block diagram of a specific embodiment of a frequency correction and time reordering circuit in accordance with the invention. 
    
    
     DESCRIPTION OF SPECIFIC EMBODIMENTS 
     Referring to FIG. 1, there is shown a flow chart or functional block diagram of a method in accordance with the invention. A source code listing of selected elements of the method according to the invention is contained in Appendix A attached, in which the word &#34;bin&#34; is replaced by the word &#34;bucket.&#34; The source code listing is internally documented to explain the functions performed. 
     In FIG. 1, the following steps are repeated for each selected frame to be displayed (a frame being typically one cycle, which could be duplicated on a time axis for display purposes). First, an unknown repetitive input signal is sampled over many periods in order to make a time record of magnitude data (Step A). For example, 512, 1024 or 2048 analog values are registered with respect to an arbitrary time base (the number of values being selected for convenience of standard FFT processing techniques). The analog values are accurately recorded. However, because the sampling is essentially random with respect to the period of the sampled waveform, the recorded time of occurrence does not provide an accurate indication of the phase position of each of the recorded values relative to the actual signal. For example, the first sample may be at time T=50 nanoseconds, the second sample may be at time T=100 nanoseconds, the third sample might be at time T=150 nanoseconds, and so forth. However, no reference is made as to the position of the samples with respect to a signal phase, since the exact timing of the samples with respect to the actual signal is only within the uncertainty between sample times since sampling occurs at a rate substantially less often than the highest frequency components of the sampled signal. In this example, the uncertainty of the signal is plus or minus 50 nanoseconds. 
     The next step is to determine the frequency of the input signal (Step B). A very rough estimate can be determined by dividing the time of the last signal edge sensed by the number of edges detected. This form of computation is insufficiently accurate to reconstruct a sampled waveform of the input signal. In the inventive method, it is necessary that the frequency be known precisely. Otherwise, the phase position of each of the samples will be inaccurately recorded. The procedure involves the setting up of an approximate representation of the waveform, thereafter performing a time to frequency transform, such as a fast Fourier transform or discrete Fourier transform, on the results to determine the frequency of the fundamental component, then using a bin interpolation procedure (Appendix A) to make the final frequency adjustment to interpolate between two sample points in the frequency output determined by the time to frequency transform. 
     The next step is to determine the phase relationship of each sample or in other words to determine the position of each sample relative to all other samples as if they were taken during the same cycle in order to place each sample&#39;s amplitude in the appropriate horizontal location (Step C). One such procedure is outlined in Appendix A wherein the raw data is reread and processed in view of the determination of the fundamental frequency and the relative position in phase of each sample. 
     Having determined the phase position of each of the samples, the samples are then reordered relative to a single frame or display period (Step D). The display period may be one complete cycle of 360° C. phase or any fraction thereof. 
     Finally, the reordered samples are presented in an output device such as a CRT (Step E). If a conventional x-y display is employed, the phase positions drive the x axis and the sampled analog values drive the Y axis. 
     This loop is repeated for each varying frame of the sampling system. 
     Referring to FIG. 2, there is shown a block diagram of a representative circuit 10 in accordance with the invention. An appropriate probe 12 receives signals which are applied to a sample and hold circuit 15 which samples at a rate determined by a clock 50. The clock 50 synchronizes all subsystems of the circuit 10. The output of the sample and hold circuit 15 is applied to an analog to digital converter 20 which loads each sampled value in an array (column) format into a first memory 25. The first memory 25 thus contains simply a table listing of the analog values in the order received with the approximate time of occurrence (within the uncertainty of the sampling system). 
     The values in the first memory 25 are then provided to means for computing the fundamental frequency 30. This frequency computing means 30 is explained hereinafter. The output of the frequency computing means 30 is the sample number &#34;n&#34; and the accurate time t i  for the i th  sample. It is supplied together with the analog values from the first memory 25 to a second memory 35. 
     The second memory 35 is for storing the sample values resorted in phase order, the phase having been accurately determined by the frequency computing means 30. The output of the second memory 35 is supplied to digital to analog converter means 40. 
     Digital to analog converter means 40 comprises a first digital to analog converter (DAC) 40A and a second digital to analog converter (DAC) 40B. The first DAC 40A reconstructs the analog values for the magnitude A i  and drives the y axis input of a CRT 45 or like output device. The second DAC 40B receives the time inputs t i  &#39; and provides the corresponding analog output signals to drive the x axis input of the CRT 45. The result is a display of one frame of a repetitive waveform reconstructed from samples of many waveforms. 
     There are many specific structures of devices which would operate according to the invention. For example, the first memory 25 and the second memory 35 may be the same physical memory device. In other forms, each memory 25, 35 may be a shift register means or other pipeline device in which data is received at an input port and provided thereafter at an output in a continuous fashion. It is also possible that the magnitude values be stored in analog fashion as in charged coupled devices. 
     If it is necessary to account for delays inherent in the processing of the frequency computing means 30, a delay device 55 may be provided which inserts for the magnitude data A i  an appropriate processing time delay between the first memory 25 and the second memory 35. This time delay assures synchronism between data in pipelined data paths. 
     Referring now to FIG. 4, there is shown one specific example of a frequency computing means 30 according to the invention. The frequency computing means 30 comprises for example a fast Fourier transform processor FFT 130 coupled to a bin interpolator 132 which in turn is coupled to an index calculator 134, which in turn is coupled to a modulo calculator 136, which in turn is coupled to a time index sorter 138. 
     Data is received at the input of the FFT 130 from the first memory 25 in the form of the array values for magnitude A i  and time t i  in the order of sampling. These values are processed in segments suitable for fast Fourier transform processing, i.e., 512, 1024 or 2048 sample points. The FFT 130 provides as its output the bin number and the value N to a bin interpolator 132, which represent the frequency of the fundamental of the sampled signal within the uncertainty of a sample bin. The bucket interpolator 132 processes the result to generate a more highly resolved value for bin number and N as its output, as explained hereinafter and in Appendix A. The bin interpolator 132 in turn provides the results to the index calculator 134. The index calculator 134 is a functional unit which divides the bin number by the value N and multiplies it by a single cycle of 360° to obtain a value K. The function is thus: 
     
         K=(BIN #)×360°/N 
    
     The index calculator 134 provides the value K to a modulo calculator 136. The modulo calculator recomputes the time sequence by the function: 
     
         t&#39;.sub.i =(i×K) mod 360°. 
    
     By this function, the value t&#39; i  is determined for each magnitude value. Thereafter, all values t&#39; i  are sorted into ascending order by the sorter 138. The output of the sorter 138 is then supplied to the second memory 35. 
     Bin interpolation is very helpful in accurately determining a frequency derived from sampled data. It has been found for example that a fast Fourier transform when used for spectrum analysis creates a difference in magnitude between a signal that falls directly on a sampling bin versus one that falls halfway between two bins, and therefore suggestions have been made to use a correction factor to compensate for frequency and magnitude errors. The correction factor is based on a window function whereby the magnitude in adjacent bins is estimated by examination of the signal at a bin in view of the shape of the main lobe of the chosen window function. 
     Stated mathematically: 
     
         M.sub.N =M.sub.S -W(x) 
    
     
         M.sub.N+1 =M.sub.S -W(1-x) 
    
     so x can be solved for from: 
     
         M.sub.N -M.sub.N+1 =W(1-x)-W(x) 
    
     where M S  is the magnitude of the signal to be estimated, M N  is the magnitude at bin N, M N+1  is the magnitude at bin N+1, and W(x) is the window response as a function of x. 
     The ideal window function W(x) exhibits a broad main lobe with little or no skirts. An alternative is to correct for droop within a window by examining and calculating a correction based on the two largest values in a local area within the frame of response data. From examination of these values, it is possible to correct not only the magnitude but to determine the true frequency to sub-bin precision, assuming there is only one signal within the response data. 
     An approximation of the window function for the Hanning window is of the form: W(x)=a 1  x 2  +a 2  x 4  where a 1  is 5.592 and a 2  is 0.405. Other coefficients apply to other window functions. The value of the prior and following magnitudes is a function of W(x), where x is to be solved for. 
     Because it is difficult to find an approximation for W(x) in the form above from known discrete magnitudes M, it is often easier to find an approximation for frequency and magnitude corrections as a function of delta M, i.e., M N  -M N+1 , by simple curve fitting to the window function. It has been determined for example that the curve fitting correction factors for frequency and magnitude for a Hanning window is as follows: 
     
         (delta F)=0.500904+8.89781E-2(delta m)+8.07123E-4(delta m.sup.2) 
    
     
         (delta M)=1.4191+0.487777(delta m)+4.19929E-2(delta m.sup.2) 
    
     Implementation of a specific embodiment of bin interpolation involves the following: First, creation of a function for magnitude correction and then the creation of a function for frequency correction. The input to both of these functions is the amplitude difference between the largest two points in a signal peak of a sampled signal. The output of both of the functions is the correction value to be applied to the larger of the two points. Since the functions are of fairly low order, these functions can be approximated by simple polynomials, such as given above. For frequency correction, the correction function 202 illustrated by FIG. 3A is preferably used. For magnitude correction, the correction function 204 illustrated by FIG. 3B is preferably used. FIG. 3A shows a correction function 202 which adds to the larger point of the pair of points in the direction toward the smaller point. In magnitude correction, as shown in FIG. 3B, the correction 204 provides that the larger point be increased by up to one-half bin droop. The position at the point labeled MAX is defined as the difference between the main lobe height and the height one bin away. 
     One procedure which computes frequency adjustments for two points in order to interpolate the position of the signal between them is given in Appendix A. In the embodiment of Appendix A, the window used is a Hanning window. Other window functions might be used such as the Hamming window or the ideal flattop. The interpolation equations shown in Appendix A thus are based on a simple curve fitting of several points. 
     This invention has now been explained with reference to specific embodiments. Other embodiments will be apparent to those of ordinary skill in this art. It is therefore not intended that the invention be limited except as indicated by the appended claims. ##SPC1##