Abstract:
A randomization unit for use in a vibrational control system for generating a series of randomized drive terms where the magnitude of each randomized drive term is the convolution of a psuedo-random sequence of bistable signals and a drive data sequence characterized by a drive power spectral density (PSD). The randomization unit does not require multipliers to perform the convolution operation.

Description:
This is a continuation of application Ser. No. 07/283,445 filed Dec. 15, 1988. 
    
    
     BACKGROUND OF THE INVENTION 
     It is often required to test a particular piece of equipment to determine whether it will properly function in an intended environment. One important factor is the vibrational energy generated in the intended environment. For example, the piece of equipment may be utilized in a space craft and therefore must survive the intense vibrational energy generated during lift off. 
     The actual vibrational energy experienced is the response to many independent vibrational sources and thus the distribution function of the amplitude over time is represented by a Gaussian distribution. Being Gaussian, this vibrational energy is completely specified by the amplitude distribution and the spectral power density function (PSD) (the first and second order statistics). In practice this reference PSD is often determined by recording the actual environmental vibrational energy and determining the PSD. 
     A well known system for subjecting a device under test (DUT) to vibrational energy is illustrated in FIG. 1. The DUT 10 is mounted on a shaker table 12 which receives a drive signal, y(t), from a driver 14. It would appear, at first glance, that the drive signal should be characterized by the reference PSD. However, the drive signal is applied to the shaker/DUT system which modifies the actual PSD of the energy applied to the DUT because of their combined transfer function. Accordingly, the drive PSD must be modified so that the PSD characterizing the actual motion of the sample is equal to the reference PSD. 
     A feedback system for generating a drive signal for subjecting the DUT to vibrational energy characterized by the reference PSD is depicted in the remainder of FIG. 1. A transducer 16 is coupled to the DUT 10 to convert the vibrations of the DUT 10 into an electrical signal. The output of the transducer is coupled to an analog to digital unit (ADC) 18 subsystem that may include signal conditioning such as filters, which has its output coupled to PSD estimation unit (S(ω)) 20. An equalization unit 22 has one input coupled to the output of the PSD estimator (S(ω)) 20 and the second input connected to a reference spectrum, S(ω), and its output coupled to the input of a frequency domain randomization unit 24. An inverse Fourier transform unit (IFT) 26 has its input coupled to the output of the frequency domain randomization unit 24 and its output coupled to a window function unit 27. The output of the window function unit 27 is coupled to the input of a time domain randomization unit 28. The output of the time domain randomization unit 28 is coupled to the input of the driver 14 which includes a DAC, power amplifier, and may include additional signal conditioning. In practice, a given system may omit the time domain randomization unit 28, the frequency domain randomization unit 24, the window function unit 27, or some combination thereof. 
     The basic purpose of the Random Vibration Control System is to have the device under test (DUT) experience a prescribed vibration as sensed at one or more points on the DUT by an appropriate transducer (accelerometer, typically). Said prescribed vibration is usually required to be Gaussian (normal) which implies that the first and second order statistics are necessary and sufficient for complete specification. The first order statistic is simply the amplitude probability density function ##EQU1## where σ is the RMS value of the random process. The second order statistic describes the frequency content or spectrum of the normally distributed signal which controls the correlation (or covariance) between successive values of the random signal. Either, the PSD, S(ω), or its Fourier transform, the auto correlation function, φ(t), are usually used to specify this statistic because these two functions constitute a Fourier transform pair 
     
         S(ω)←→φ(t) 
    
     Because digital signal processing techniques are used in these control systems, uniformly sampled time domain and frequency domain functions are used in the control loop: analog to digital converters and discrete Fourier transform usually providing the means for obtaining sampled functions. 
     It is important to recognize that an essential feature of these normal random signals is that the discrete Fourier transform X(k) of a sample segment of x(n), consisting of N samples duration given by ##EQU2## is a complex function of the frequency index, k, such that ##EQU3## Note that if the set of values {x(n)} are randomly distributed, the real part P x  (k) and the imaginary part, Q x  (k), being weighted sums of these randomly distributed samples must also be random. Further, if x(n) is normally distributed, P x  (k) and Q x  (k) will be, for certain, normally distributed. Thus, the Fourier transform of a normally distributed signal is a complex random function with normally distributed real and imaginary parts. P x  (k) and Q x  (k) are, also, uncorrelated because the weights in cos 2πkn/N and sin 2πkn/N used to obtain each part are orthogonal (uncorrelated). 
     If the complex function, X(k), is expressed in polar coordinate form, ##EQU4## then θ(k) can be shown to be uniformly distributed on the interval 0≦θ≦2π, while |X(k)| is Rayleigh distributed. 
     In order to establish control of the vibrations, the spectrum of the vibrations must be estimated and compared with the desired spectrum so that the spectral content of the drive signal may be appropriately tailored to produce the desired result. Conceptually, it is a straightforward problem: estimate the transfer function by estimating the transfer function between the input to the DAC 14 and the output of the ADC 18 from which the necessary drive signal may be determined. In practice, problems may arise because the vibrations are random, the system may be noisy and may also be nonlinear. 
     Because the sensed response, z(n) is random, the Fourier transform, in accordance with the previous discussion, is also random. This means that each of the individual estimated spectrums or periodograms will be a random function. Specifically, the spectral estimator of the response signal that we desire is ##EQU5## τ is the sampling interval and z m  (n) is the m th  response data frame. Hence ##EQU6## where P zm   2  (k) and Q zm   2  (k) are each squared Gaussian variables that are, therefore, chi-squared distributed with one degree of freedom each. The sum of the terms ##EQU7## known as the &#34;raw&#34; spectral estimates or periodograms, is also chi-squared distributed but with two degrees of freedom, and has a relative variance of unity with respect to its true mean value. The purpose for averaging over an infinite ensemble of such terms is to estimate S z  (k) by reducing these variances. In practice a finite weighted sum of &#34;raw&#34; spectral estimates, S zm  (k), are used as a basis of control. For example 
     
         S.sub.z,m+1 (k)=αS.sub.z,m (k)+(1-α) S.sub.z,m+1 (k) 
    
     where 
     S z ,m (k) is the m th  spectral estimate; ##EQU8## α is a &#34;discount&#34; averaging factor (0&lt;α&lt;1). Thus, this spectral estimate will have errors (variance) because of the finite number of degrees of freedom used in the estimate. The variance decreases inversely with m. Additional variance may be introduced by such non-linear PSD estimating techniques as windowing of the sampled response data, z(n), for purposes of reducing leakage in the spectral estimate. 
     Because the system may also be noisy and non-linear a stochastic approximation algorithm is used that introduces less than the total correction implied by the latest average response spectrum. Thus, a typical correction for the drive Fourier transform magnitude may take on the form ##EQU9## where 0&lt;ε&lt;1/2 is the stochastic weighting factor. This procedure results in a deterministic magnitude function whereas what we desire is 
     
         X.sub.m+1 (k)=|X.sub.m+1 (k)|e.sup.j.spsp.θ.sup. m+1.spsp.(k) = P.sub.x,m+1 (k)+jQ.sub.x,m+1 (k) 
    
     whee |X m+1  (k)| is Rayleigh distributed, θ m+1  (k) is uniformly distributed or equivalently P x ,m+1 and Q x ,m+1 are independently normally distributed in order to produce a proper Gaussian distributed drive function, x m+1  (n) after taking the inverse Fourier transform of X m+1  (k). 
     Prior art has established that an acceptable approximation to the normal drive signal may be achieved by only randomizing X(k) by introducing a random phase function θ(k). See U.S. Pat. No. 3,710,082. Because the central limit theorem asserts that the sum of a large number of random variables tends to be Gaussian, almost independent of their particular distributions, the inverse Fourier transform, x m+1  (n), will also be approximately Gaussian. This procedure has been dubbed &#34;frequency domain&#34; randomization even though only the phase is randomized and not the amplitude also. 
     A further variant was also introduced by Sloane et al., see U.S. Pat. No. 3,848,115 that has been popularly referred to as time domain randomization that is based on randomizing phase only by random circular rotation of x(n). This procedure also only randomizes phase, and not both phase and amplitude. 
     The time domain function, x m  (n), is a finite length array of N data points that could be concatenated with other such arrays to produce a continuous stream, x(m,n) such that ##EQU10## where 0≦n≦N-1. However, it was found in Sloane et al. that the uniform weighting implied by construction of x(m,n) introduced unacceptable spectral leakage that limited the ability to control systems exhibiting large sharp resonances or being able to control when the desired reference spectrum, S(k), had large dynamic range characteristics. As a result, non-uniform windowing of x m  (n) was introduced in order to control this spectrum. 
     FIG. 2A and 2B depict systems utilizing phase domain and time domain randomization, respectively, followed by windowing and overlap and add. 
     From a practical point of view, non-uniform windowing by itself could not be arbitrarily applied to each segment, x m  (n), because concatenation of such a sequence, ##EQU11## where w(n) is the window function, would produce a non-stationary output sequence whose instantaneous power was a periodic function proportional to w 2  (n) in violation of the usual requirement that, at steady-state or at convergence of control to the proper spectrum, the drive function must also be a stationary normal process. This led to the introduction of overlapping in addition to windowing. A further implicit constraint was thus established on the choice of window functions because, for example, a 50% overlap of windowed arrays required that  where c is a constant. For the 50% overlap case, this restricted the choice of window to the half-sine window ##EQU12## Unfortunately, this constraint on the window function restricted the choice to one of the least effective of all window functions. Thus, we conclude that the state of the art is deficient on several counts: 
     1. Dynamic range is limited because of the windowing constraint. 
     2. True phase and frequency randomization are lacking. 
     3. Even under the best circumstances, the second order statistics are non-stationary because the overlap and windowing operations result in a correlation function that is a function of two times, φ(t 1 ,t 2 ), rather than the stationary function which is only dependent on the time difference T=t 1  -t 2 , or φ(t 1  -t 2 )=φ(T). 
     SUMMARY OF THE INVENTION 
     The objective of this invention is to provide true randomization of both amplitude and phase and optimal leakage control in a practical and efficient manner by either time domain or frequency domain operations that also preserve the quality of the second order statistics. Specifically, the method involves digital filtering of a random sequence by a windowed impulse response function that is determined from a comparison of response and control spectrums. Two methods, one based on time-domain convolution and the second on the equivalent frequency-domain operation, are disclosed that lead to comparably equivalent results. 
     According to one aspect of the invention, a time domain digital filter is implemented so that each term in the drive data sequence is gated by a corresponding bistable signal in the randomized sequence. The value of the gated signal is zero if the corresponding bistable signal is in a first state and the value is equal to the corresponding term in the drive data sequence if the corresponding bistable signal is in a second state. The gated signals are then summed to form convolution term. These convolution terms are utilized as the randomized drive data sequence. 
     The result of the gating and summing operations is to generate a linear combination of the terms in the drive data sequence, with the coefficient of each term equal to zero or one depending on the state of the corresponding bistable signal. Accordingly, the convolution is performed without utilizing multipliers. According to a further aspect of this embodiment, the drive and randomized sequences are stored in shift registers with the parallel outputs of the shift register coupled to the inputs of gating elements for gating the terms in the drive data sequence as described above. 
     According to another aspect of the invention, the bistable signals in the randomized control sequence are utilized to control an adder/subtracter unit (ASU) having a first input port for receiving a corresponding term in the drive data sequence and second input port coupled to an output port of the ASU. A linear combination is formed with each term either added or subtracted to the combination as determined by the state of the corresponding bistable signal. 
     Thus, a linear combination of the terms of the drive data sequence are formed with the coefficient of each term being -1 or +1 as determined by the state of the corresponding bistable signal. Accordingly, the convolution is performed without using multipliers. 
     Additional features and advantages of the invention will become apparent in view of the drawings and the following detailed description. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a schematic diagram of a typical prior art vibration control system; 
     FIGS. 2A and 2B are schematic diagrams of prior art systems utilizing windowing and overlap and add; 
     FIG. 3 is a functional block diagram of the invention. 
     FIG. 4 is a block diagram of a randomization system; 
     FIG. 5 is a block diagram of the invention; 
     FIG. 6 is a schematic diagram of a unipolar embodiment of the present invention; 
     FIG. 7 is a detailed schematic diagram of a gating element utilized in the embodiment of FIG. 4; 
     FIG. 8 is a schematic diagram of a bipolar embodiment of the present invention; 
     FIG. 9 is a timing diagram illustrating the operation of the embodiment of FIG. 6; and 
     FIG. 10 is a schematic diagram of a vibration control system using a frequency domain digital filter. 
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     Referring now to the drawings, where like reference numerals identify identical or corresponding parts, FIG. 3 is a functional block diagram of an embodiment of the present invention. Referring to FIG. 1, a vibration control system having many elements common to the to the prior art system described with reference to FIG. 1 is depicted. 
     In FIG. 3, a window function generator 8 is coupled to a first input of a windowing unit 10. The windowing unit 10 has a second input coupled to the output of the IFT and receives a time domain digital drive data sequence, x(n), characterized by the reference PSD. The output of the windowing unit 10 is coupled to first input of a special digital filter 12. The second input of digital filter 12 is coupled to the output of a psuedo-random number generator (PRNG) 13. 
     The operation of the embodiment of FIG. 3 will now be described. The window selection function generated by the windowing function generator 8 is not limited as in the prior art. Efficient windows with low spectral leakage, such as Dolph-Tchebychev or Kaiser-Bessel windows, may be selected. The windowed drive data sequence functions as a digital representation of the impulse response of the digital filter 12. Accordingly, the result of filtering the output of the PRNG 13 is to generate a randomized drive data sequence characterized by the reference PSD, having low spectral leakage, and a Gaussian amplitude distribution. The average power of the randomized drive data sequence is a constant, regardless of the form of the window function. Specific time and frequency domain implementations of the digital filter 12 are described below. 
     FIG. 5 is a block diagram of an embodiment of the invention. Referring to FIG. 5, a unit 40 for generating a random, or psuedo-random, sequence of bistable signals has an output port coupled to the serial input port of a first N-stage, serial in, parallel out shift register 42. The design of the afore-mentioned random signal unit is generally well-known and described, for example, in the above-reference patents. 
     The first shift register also has N parallel output ports for transferring the signals stored in respective shift register stages 44 and a serial input port. The serial input port transfers a signal to the first stage. 
     A second N-stage multibit, serial in, parallel out shift register memory 46 includes a multibit serial input port, for transferring signal to its first stage, coupled to the y(n) output port of unit 10, and N parallel multibit output ports each coupled to one multibit stage 47, each multibit stage, 47, holding one output word of unit 10. 
     A timing unit 48 generates a CLK signal for generating a series of CLK signal pulses that define a sequence of clock cycles. A linear combiner 50 receives the bistable signals stored in the first shift register 42 and the terms of the impulse response sequence stored in the second shift register 46. 
     In operation, the linear combiner utilizes the bistable signals, x(n-k), stored in the first shift register 42 as control signals to determine the contribution of each term, y&#39;(k), in the drive data sequence to a term y(n) in the randomized drive data sequence. Two embodiments will be described. The first, described with reference to FIG. 6, utilizes the bistable signals as gating signals to control a set of gating elements in a parallel gating network. This embodiment generates a linear combination with the y&#39;(k) terms multiplied by unipolar coefficients. The second, described with reference to FIG. 6, utilizes the random bistable signals to control an adder/subtracter unit (ASU) and to generate a linear combination with the y&#39;(k) terms multiplied by bipolar coefficients. 
     FIG. 6 illustrates an embodiment for forming the linear combination of eq. 6 with coefficient term x(n-k) equal to 0 or 1 (unipolar coefficient). Referring now to FIG. 6, a parallel gating network 70 includes a set of N gating elements 72. Each gating element has two input ports coupled to a single bit input from each element 44 of shift register 42 and to a p-bit input from each corresponding element 47 of the shift register memory 46. The stages of each shift register are identified by the index k, and each shift register coupled to the serial input port identified by the index k=N-1. 
     Each gating element 72 also has a multibit output port coupled to one of the input ports of adder 74. The adder 74 also has an input port coupled to the output port of an offset correction unit 76. The output port of the adder 74 functions as the output port of the digital filter. 
     The operation of the system of FIG. 6 will now be described. At time corresponding to the index n, the contents of shift register 42 contain the bistable values x(n-N+1) through x(n) corresponding to x(n-k), 0≦k≦N-1 where k is the address index. Simultaneously, the contents of the shift register memory 70 contains the last (m th ) N-word array y&#39; m  (n) from unit 10 stored in locations 0≦k≦N-1 with k=n. The k th  stored binary value of memory 42 and the each bit of the k th  stored word of memory 70 logically interact through gate 72(k) by the Boolean &#34;and&#34; operation to produce an output equal to y&#39; m  (k) if x(n-k) is true (x(n-k)=1) and zero if x(n-k) is false (x(n-k)=0). thus, this is equivalent to forming the product ##EQU13## The output of the N gating elements 72 are digitally summed by the adder/accumulator 74 to form ##EQU14## This expression represents the discrete convolution of x(n) with y m  (n). Because both x(n) and Y m  (n) are time varying, either sequence (array) may be considered as the convolution kernel. However, because the spectral update changes are generally slower than the correlation time of the PRNG 40, the above expression takes on the meaning of a time varying filter with an impulse response of y&#39; m  (n) operating on a stationary random sequence x(n). Significantly, no multiplier is required to generate the product that is summed in the above expression for y m  (n). 
     The generation of the filter ouput sequence requires that a new sequence of N words representing x(n) be available for forming the next output word y m  (n+1). It is sufficient that the next noise generator sequence be 
     
         [x(n+1), x(n), x(n-1) . . . x(n-N+2)] 
    
     a shifted N element array corresponding to the old array less then oldest element x(n-N+1) but with x(n+1) occupying the newest data position, k=0. 
     FIG. 7 is a detailed schematic diagram of a gating element 72. Each gating element 72 includes p AND gates 80 having a first input port for receiving the signal x(n-k), a second input port coupled to a corresponding conductor of the first p-bit databus 78, and an output port coupled to a corresponding conductor of the second p-bit databus 79. 
     In operation, each AND gate 80 forms the Boolean product of the one bit unipolar signal x(n-k) and a single bit of y&#39;(k). The complete set of p Boolean products forms the arithmetic product of x(n-k) and y&#39;(k). 
     Because x(n) and y&#39;(n) are statically independent, the mean or expected value at the output of the filter is ##EQU15## where E{w(n)} is the average value of w(n), while usually E{x m  (n)}=0 because the test is run at a net displacement of zero. Thus, ##EQU16## This implies that the means value (DC-term) out of the filter is zero as it should be for a Gaussian zero means drive signal even though the unipolar PRNG sequence has a net average value of 1/2 because the filter impulse response function corresponds to a bandpass filter with zero gain at zero frequency. 
     FIG. 8 illustrates an embodiment for forming the linear combination of eq. 6 with the coefficient term x(n-k) equal to +1 or -1 (bipolar coefficient). Referring to FIG. 6, the stages of the first and second shift registers are identified by the index k as described above with reference to FIG. 4. An adder/subtracter unit (ASU) 90 has first and second input ports, an add/subtractor control port, and an enable control port, and an output port. The control port is coupled to the serial output port of the first shift register 42 by a dataline and the first input port is coupled to the serial output port of the second shift register 44 by a first p-bit databus. 
     The accumulator register 92 has an input port coupled to the output port of the ASU 90, a clear input for receiving a CLR signal, and an output port connected to a selector switch S3 95 that is controlled by the signal 102 which switches the output of the accumulator register 92 to the output buffer 94 each time a new output drive sample is accumulated in register 92. 
     The serial input port of the first shift register 42 is coupled to the output of the bistable noise generator 40 or to its own serial output port via selector switch S1. The serial input port of the second shift register 44 receives the drive data sequence, y&#39;(n), or its own serial output via selector switch S2. The serial output port of the second SR 44 is also coupled to the first input port of ASU 90 by a first p-bit databus. 
     The output port of the accumulator register 92 is coupled to the second input port of the ASU 90 by a second multibit databus or to the filter output port via selector switch S3 95. 
     FIG. 9 is a timing diagram for the system of FIG. 8. Referring to FIG. 9, a CLK signal 100 is generated and the ENABLE, 102, SELECT/STORE 104, and SHIFT 106 are derived therefrom. The timing control signals have a periodicity of N+1 cycles of the CLK signal 100. The CLK pulses are identified by the index number n. When the ENABLE 102 is high, the ASU 90 is inoperative. The shift 106 signal increments the two shift registers 42 and 44. When 102 is true, the contents of the accumulator register is stored in the output buffer 94. The selector switch S1 96, is normally connected to its own output. When control signal 104 is true, however, the latest output from the PRNG 40 is stored in shift register 42, once during each cycle of N+1 clock cycles. 
     The operation of the bipolar system of FIG. 8 will now be described with reference to FIGS. 8 and 9. It is assumed that N terms of the random bistable sequence, x(n), are stored in the first shift register 42 and that N terms of the drive data sequence, y&#39;(n), are stored in the second shift register 44. 
     The first input of the ASU 90 is coupled to the accumulator register 92 via S3 95, while the second input port is coupled to the second shift register 44 output. The selector switch S2, allows a means for updating the filter kernel, y&#39; m  (n), stored in shift register 44 by means of the S2 97 control line labelled LOAD FILTER KERNEL. 
     When the control signal 102 is true, the selector S3 95 connects the accumulator register 92 to the output buffer 94. Usually, in practice, the impulse response of the filter y&#39; m  (n) is not updated as frequently as the contents of shift register 42. 
     After N accumulations, n=N and control signal 102 switches S3 95, connecting the accumulator register R92 to the write port of the output buffer memory 94 which is simultaneously enabled to write by control signal 102. Also, 102 disables the ASU 90 allowing the oldest stored value in SR 42 to be shifted out and lost while the latest output of PRNG 40 is switched into SR 42 via S1 96. Control signal 104 is true during the next half clock cycle thus clearing the accumulator register R92 in preparation for the next N accumulations. 
     From n=0 to n=N-1, the contents of shift register 42 and 44 are incremented by the shift 106 signal which is staggered by a half cycle with respect to the CLK 100 signal. The accumulator register 92 contents are modified by the contents of shift registers 42 and 44 once each CLK 100 period: the output of shift register 44 being added to the contents of the accumulator when the ADD/SUBT signal is true and subtracted when false. Thus, this has the effect of multiplying, term by term, the contents of shift registers 42 and 44 and accumulating the resulting products. In this way, the contents of the accumulator register 92 at the end of N shift cycles is ##EQU17## which corresponds to the desired convolution. Accordingly, a new value of the convolution of x(n) and y&#39;(n) is directly formed in N clock cycles. 
     From the above it is apparent that one term, y(n), in the randomized drive data sequence will be generated every N+1 clock cycles. Further, although the effect of each random signal is to multiply the corresponding y&#39;(n) term by +1 or -1, no multipliers are required. Further, because the x(n) terms are used as control signals to the ASU 90, these signals may be either bipolar or unipolar, although the effect is multiplication by a bipolar signal. No offset term is required because the time average of the bipolar coefficients is zero. 
     Both methods, i.e., the method described earlier using a unipolar (0,1) PRNG 100 signal or the latter method using bipolar (±1) interpretations of the PRNG 100 output signals are functionally equivalent except for a scale factor: the unipolar method being equivalent to an add or not-add instruction results in an RMS scale factor of √2/2 relative to the bipolar add/subtract method because the product being accumulated is, on average, zero half the time. 
     We have described a method and apparatus for an automatic digital random vibration control system capable of producing at equilibrium the required normal (Gaussian) signal with appropriate stationary first and second order statistics. We will now demonstrate that the method described for introducing windowing of the drive signal is superior to the current state of the art. 
     The output y&#39; m  (n) of the digital filter 12 of FIG. 3, when expressed as a convolution may be represented as 
     
         y(n)=y.sub.m &#39;(n)*x(n)                                     Eq. 7 
    
     where y m  (n) is the impulse response of the filter, x(n) the binary noise sequence, and the asterisk, *, represents the convolution operation. 
     The discrete Fourier transformation of Eq. 7 is the product 
     
         Y(k)=Y.sub.m &#39;(k)·X(k)                            Eq. 8 
    
     where the Y(k), Y m  (k) and X(k) are the Fourier transforms of y(n), y m  (n), and x(n) respectively. Because 
     
         y.sub.m &#39;(n)=ω(n)·x.sub.m (n)               Eq.9 
    
     where ω(n) is the window function and x m  (n) is the inverse Fourier transform of ##EQU18## where S(k) is the updated drive spectrum and θ(k) is the standardized phase function. Hence, the Fourier transform of Eq. 9 is 
     
         Y.sub.m &#39;(k)=W(k)*X.sub.m (k) 
    
     or 
     
         Y.sub.m &#39;(k)=σ.sub.i W(i) X.sub.m (k-i)              Eq. 10 
    
     Substituting Eq. 10 into Eq. 8 yields 
     
         Y(k)=X(k)·ΣW(i)X.sub.m (k-i)                Eq. 11 
    
     Thus, the drive spectrum is proportional to 
     
         E{|Y.sub.m (k)|.sup.2 }=E{|X(k) ΣW(i)X.sub.m (k-i)|.sup.2 }                Eq. 12 
    
     
         Because X(k) and ΣW(k)X.sub.m (k-i) are uncorrelated, ##EQU19## Because for the bipolar case E{|X(k)|.sup.2 }=1 and because the binary noise sequence has a white spectrum and the system is at equilibrium, ##EQU20## so that Eq. 13 becomes ##EQU21## 
    
     This expression is proportional to the expression that would be obtained by widowing and overlapping. 
     To demonstrate this assertion, consider a single windowed drive vector x m  (n)w(n). Its contribution to the spectrum would be proportional to ##EQU22## Because x m  (n) has a white spectrum ##EQU23## 
     We conclude both methods introduce windowing into the drive spectrum and thus control leakage which limits dynamic range. However, the conventional overlap and add method is limited in its choice of window functions, typically settling for a Hanning window while the method described is free to use any window function because stationarity of the mean square value is maintained. 
     A frequency domain equivalent to the time domain filter may be implemented by using a well known Fourier transform relationship corresponding to the time domain operation of filtering (convolution). Specifically, the convolution of two time domain functions, say h(t) and x(t), corresponds to the multiplication of the corresponding Fourier transformations, H(ω) and X(ω): ##EQU24## This Fourier transform pair implies that frequency domain operations may be used for purposes of filtering. 
     Thus, the continuous process of time domain filtering may be approximated by a piece-wise frequency domain operation in which the data sequence, x(t), is partitioned into concatenated sections of finite duration, T, such that ##EQU25## 
     Unfortunately, the use of discrete Fourier Transforms (DFTs) such as the well known fast Fourier transform or FFT introduces a further complication. Because DFTs are sampled versions of the Fourier transform they implicitly correspond to periodic versions of the original time domain data. This necessitates the appending of zeros in order to avoid &#34;circular&#34; convolution as described by T. G. Stockham, Jr. in &#34;Digital Processing of Signals&#34;, Gold and Rader, McGraw-Hill, 1969, pp. 203-213. 
     A preferred embodiment of this technique is shown in the following block diagram of a Vibration Control System, FIG. 10, wherein the time domain filter has been replaced by the frequency domain filter shown within the dotted outline. The data to be filtered is supplied by the psuedo-random noise generator (PRNG). The filter characteristic is determined by the impulse response, h(n), from the array multiply operation, w(n)·x(n). If this sequence is N words long, N zeros are appended to avoid circular convolution, before the Fourier transform operation. Note that the FFT used must now operate on 2N real data points to produce H(k). Similarly, the PRNG output is segmented to consist of N data values and N zeros from which Y(K) is obtained via a 2N real FFT. The product X(K)·H(k) is formed and its inverse Fourier transform (IFT) is taken to produce a time domain sequence of 2N duration. The first N words are added to the previously obtained last N words (overlap and add operation) while the last N words are temporarily stored for the next overlap and add operation. The output is fed to the system DAC as previously described. 
     Many other variants of this frequency domain filtering operation such as &#34;select-saving&#34; are possible and will suggest themselves to those versed in the art of digital signal processing. (See H. D. Helms, &#34;Fast Fourier Transform Method of Computing Difference Equations and Simulating Filters&#34;, IEEE Trans. Audio Electroacoust., Vol. AU-15, pp. 85-90, June 1967). (Also see Otnes and Enochson, &#34;Applied Time Series Analysis&#34;, Vol. 1, John Wiley and Sons, 1978, pp. 280-306). 
     Several embodiments of the invention have been described. Modifications and substitutions will now be apparent to persons of ordinary skill in the art. In particular, the functions of the hard-wired embodiments of FIGS. 6 and 8 may implemented by a programmed digital computer or microprocessor. The randomized bistable sequence could be stored or generated in the computer and the drive data sequences could be stored in the computer memory and the processor could be programmed to supply the correct terms to the gating elements of FIG. 4 or to the control port and first input port of ASU 90. Additionally, the processor could control the generation of the required control signals. Other hard-wired embodiments are also within the scope of the invention. Accordingly, the invention is not intended to be limited except as provided by the appended claims