Characterization and simulation of gunﬁre with wavelets

Gunﬁre is used as an example to show how the wavelet transform can be used to characterize and simulate nonstationary random events when an ensemble of events is available. The structural response to nearby ﬁring of a high-ﬁring rate gun has been characterized in several ways as a nonstationary random process. The current paper will explore a method to describe the nonstationary random process using a wavelet transform. The gunﬁre record is broken up into a sequence of transient waveforms each representing the response to the ﬁring of a single round. A wavelet transform is performed on each of these records. The gunﬁre is simulated by generating realizations of records of a single-round ﬁring by computing an inverse wavelet transform from Gaussian random coefﬁ-cients with the same mean and standard deviation as those estimated from the previously analyzed gunﬁre record. The individual records are assembled into a realization of many rounds ﬁring. A second-order correction of the probability density function is accomplished with a zero memory nonlinear function. The method is straightforward, easy to implement, and produces a simulated record much like the measured gunﬁre record.


Introduction
The simulation of the response of nearby external stores to the firing of high-fire rate guns has been described using several methods by two authors [2][3][4][5]8]. The first attempts, which were not completely satisfactory, used a complex periodic waveform to simulate the gunfire. Merritt used short-time Fourier spectra to simulate the gunfire, and Smallwood used the product model for nonstationary random vibration together with band-limited temporal moments to describe the waveform. Both of these methods result in very good simulation but are somewhat complicated to implement. This paper will illustrate a simpler approach using wavelets.
Gunfire belongs to a class of environments characterized by waveforms that are a mixture of a deterministic component and a nonstationary random component. The large random component complicates the characterization of these waveforms. To separate the deterministic part, requires averaging to reduce the effects of the random contribution. If the deterministic part is periodic, averaging can be performed in the time domain over an integer number of cycles. If the deterministic part is not periodic, an ensemble average must be used. Gunfire is between these two extremes; the firing rate of the guns is almost periodic. Both Merritt and Smallwood estimated the deterministic part by breaking apart the gunfire record into records of one-round firing and performed an ensemble average of the resulting records. The estimate of the deterministic part is then subtracted from the records, leaving an estimate of a zero mean nonstationary random residual.
Characterizing the residual is difficult because of the short length of the records. The uncertainty theorem limits the product of the time resolution and the frequency resolution of the characterization [1,6]. If we attempt to analyze a single record with a resolution near the theoretical limit, we cannot achieve more than a single statistical degree of freedom for each of the coefficients. We can calculate the Fourier coefficients and resolve the frequencies with a resolution given by the frequency spacing of the coefficients, but we lose any time resolution. The time resolution for all components is the length of the data record. We lose any information about the nonstationary character of the random waveform. We can improve the time resolution in several ways. Merritt breaks each record into smaller time slices and performs a spectral estimate in each time slice with the resulting loss of frequency resolution. Smallwood filtered the records using rather wide bandpass filters and used temporal moments to char- The wavelet transform offers a third alternative to attack the problem of time and frequency resolution.

Wavelets
Wavelets are a new method (which became popular about ten years ago) for analyzing data. Unfortunately, space will permit only a brief description of wavelets. Strang [10] contains a good introduction in an article. A book by Strang and Nguyen [11] contains a thorough discussion of wavelets. Wavelets are a powerful new way to look at transient data. Basically a function is described as a combination of basis functions, w jk (t), the wavelets Equation (1) is an inverse transform: If you know the coefficients, b jk , you can reconstruct the original waveform. The forward transform exists and is typically implemented with a pair of filters as will be described later. A typical wavelet is compressed in time and shifted: The wavelet is normalized with the factor 2 j/2 on the interval [k∆t, (k + N )∆t] and is zero outside the interval, where ∆t is the sampling interval. The basis wavelet (when j and k are zero) w 00 (t) = w(t) often has a quite complicated shape but is completely described by a pair of filters: a high pass filter and a low pass filter. A number of useful wavelets have been defined, and one of the problems of a user is to pick the most appropriate wavelet for a particular application. The wavelets typically have compact support. Simply this means that the wavelets are nonzero only over a short finite period of time. N is typically small.
In practice the continuous form of the wavelet, Eqs (1) and (2), is seldom used; we almost always are dealing with finite waveforms sampled at discrete times. The wavelet coefficients C are given by Z is the number of levels and N is the normalized length of the wavelet.
An advantage of wavelets is that they provide a description of a waveform that is localized in both time and frequency. Although users of wavelets seldom refer to frequency, they use the term scale instead. Some of the wavelets are members of a family of wavelets identified by the order number. In general, the higherorder wavelets provide better frequency localization at the expense of time localization. High frequencies are localized in time more precisely than low frequencies.
The product of the rms bandwidth and the rms duration of a wavelet is almost constant. The wavelet chosen for this paper was the Daubechies wavelet of order 20, abbreviated as the db20 wavelet in the MATLAB Wavelet Toolbox. The basic wavelet is plotted as Fig. 1. This wavelet provides an uncertainty product of rms duration and rms bandwidth of about 5 times the theoretical minimum. It is one of the better wavelets in this regard.
Each level of decomposition separates the signal into two bands, with a high-pass filter and a low pass filter. The output of the high-pass filter is called the "detail", and the output of the low-pass filter is called the "approximation". The approximation and the detail are decimated by two, and the process is repeated for the approximation for the next level. The detail is not processed again in this paper because we will be satis- fied with an octave band characterization of the signal content. Remarkably the desired wavelet coefficients are the decimated detail. At first glance it would appear the decimation of the detail would cause serious aliasing problems, but this is not the case. The frequencies are aliased, but the two filters are carefully constructed such that the information lost in the decimation of the detail is retained in the approximation. Leakage, which is aliased in the approximation, is retained in the detail.
Exact reconstruction is possible from the approximation and details. This is one of the advantages of wavelet decomposition as contrasted with decomposition using filter banks. If more traditional filter banks are used for decomposition the decomposed waveforms cannot be added to reconstruct the original waveform. Exact reconstruction is not usually possible.
Since each level of approximation has about half the bandwidth of the previous level, the upper frequency limit of the n-th level of approximation is about f s /2 n+1 , where f s is the sampling frequency. Thus the level chosen provides frequency localization. In this case the approximation contains the low frequencies, and the detail contains the high frequencies. We also would like localization in time to resolve the nonstationary character of the waveforms. The wavelet transform inherently provides the localization in time. Each wavelet coefficient is a sample in time of a filtered and decimated waveform.
If we take the wavelet transform of a set of independent samples, the resulting wavelet coefficients are asymptotically independent. We can see this because the wavelet transform is invertible, and the number of wavelet coefficients is essentially the same as the number of samples in the original waveform. We have neither lost nor added information with the transform. The wavelet transform is a linear transform with no loss of information.
This implies that each coefficient has about one statistical degree of freedom. To improve our estimates, we must perform an average. If an ensemble is avail- able, an ensemble average is the obvious choice. If a sufficient number of realizations are not available, averaging wavelet coefficients in the time domain will improve the estimates at the expense of time resolution. At some point the advantage of the wavelet transform over a short-term Fourier transform representation or the output of a filter bank for the characterization is lost.
In the analogous Fourier transform, the real and imaginary parts of the DFT each have one degree of freedom. The difference is that for the Fourier transform, all frequency components have the same time resolution. If the input sequence consists of independent numbers, we need both the real and imaginary numbers of the DFT to reconstruct the original waveform. When we find the magnitude squared of a spec-tral line, we have two degrees of freedom. This leads to the Chi squared distribution, with two degrees of freedom for the spectral estimation. An ensemble average is used to increase the degrees of freedom.
If we perform a simple ensemble average of the wavelet coefficients, the result will be near zero for a zero mean process. Noting that the wavelet transform can be represented as a linear transformation where {b} j are the wavelet coefficients at level j, H j is the wavelet linear transformation at level j, and {x} is the sampled waveform. Take the expected value Since the transformation has an inverse, the inverse transform of the expected value of the coefficients is the same as the expected value of the original waveform. This implies that the deterministic part of the waveform can be estimated in an equivalent manner by taking the mean of the ensemble or by taking the inverse transform of the mean of the wavelet coefficients. If x is a nonzero mean process, the mean, or the reconstruction from the mean of the wavelet coefficients, is an estimate of the deterministic part of the waveform. If x is a zero mean process, the mean of the wavelet coefficients b is also zero.

Characterizing gunfire using wavelets
The characterization begins just as Merritt and Smallwood did in the previous papers. The original gunfire waveform, supplied by Ron Merritt, is shown as the top figure in Fig. 2. The bottom half of the figure will be discussed later.

The gunfire record is broken into records, each
representing the firing of a single round. 2. The start time for each record is determined by performing a cross correlation between all of the records and one record arbitrarily chosen as the reference. Time is shifted for each record to maximize the cross correlation at time equal to zero. The variation in the time between rounds firing is characterized with a mean and standard deviation. 3. The mean record is computed by averaging all the records. The mean is subtracted from the records leaving a zero mean residual. The original record for the 25th round is plotted as the top  This is called the residual record. Notice that the mean correlates well with the original for about the first 5 ms. After 5 ms the high frequency correlation is lost. This is characteristic of a narrow band random signal. The residual looks like a broad band random that looks to the eye to be approximately stationary over the duration of the record. We will see later that the residual is not stationary over the duration of the record. The characterization to this point has followed the method of the previous paper [8]. 4. The wavelet coefficients are computed for each residual record. A Daubechies order 20 (db20) wavelet, with 6 levels, was used. The mean and standard deviation for each wavelet coefficient are computed over the ensemble of 50 records. It was confirmed that the mean was very near zero. A higher-order wavelet, in general, will give a little better frequency localization, without as much leakage, than a lower-order wavelet. The number of levels was chosen to be the minimum number needed for an adequate decomposition.
Computing the mean and standard deviation of the original waveform before the mean was removed would generate an equivalent characterization. In this case a reconstruction, which included the mean, would reconstruct the mean record.
A plot of the standard deviation of the coefficients is shown as Fig. 4. To facilitate understanding, the co- efficients have been zero shifted, and a multiple of 10 was added to each level. The Level 1 detail is quite small, and the values were multiplied by 10 for plotting. The coefficients have been plotted by level with the horizontal axis as an approximate time axis. The approximate bandwidth of each level of coefficients is included in the legend. The standard deviation of the coefficients also bears a resemblance to the variation in mean square from Smallwood [8, Fig. 4]. The largest values in the Level 1 detail are edge effects. The next three levels (D2, D3, and D4) show a definite nonstationary character. The lowest frequency levels (D5, D6, and A6) have a more uniform distribution in time.
The probability density of the coefficients normalized to a unity standard deviation was estimated using Kernel density estimator [7] as Fig. 5. Also plotted is a Gaussian distribution with a unity standard deviation. We can see the distribution of the wavelet coefficients is nearly Gaussian. We can check for correlation of the coefficients by estimating the autocorrelation of the coefficients at each level. The results are plotted as Fig. 6. The correlations are near zero for lags removed from zero except Level 2. This implies that most of the coefficients are nearly uncorrelated. We could probably reduce the correlation of coefficients by using wavelet packets and resolving the data with a finer frequency resolution with a loss of time resolution, but this will not be done in this paper.
Characterizing the gunfire using wavelets has several advantages over the previous methods. The methods of Merritt use short time Fourier transforms. The gunfire is broken into time segments and each segment is characterized with a random realization. He has problems in the transitions between time segments. There are also problems with picking the appropriate time segments. The time resolution is independent of frequency, all frequencies are resolved with the same time resolution. This requires the time segments be chosen compatible with the lowest frequency of interest. The band-limited method [8] faces similar problems in the frequency domain. Using this method the frequency bands must be chosen, which are not always an obvious choice. The filtering of the waveforms into bands does not allow exact reconstruction. If the bands are chosen too narrow, the filtering spreads the waveforms in time giving the appearance of duration's which are longer than reality. Using band-limited methods leads to problems with transitions between the frequency bands. Wavelets overcome this problem providing a natural way to provide both time and frequency localization. Wavelets are a very effective compromise between the conflicting demands of time and frequency localization.

Simulation of gunfire using wavelets
The characterization developed in the previous section can now be used in a simple manner to generate realizations of simulated records.

A realization of the residual is formed by tak-
ing inverse wavelet transform of a set of random wavelet coefficients with the same standard deviation of the coefficients as found from the characterization. A zero mean Gaussian distribution was used. The coefficients were assumed to be independent. As discussed above, this is only approximately true. The effects if this assumption will be discussed later. It may have been possible to improve the simulation by introducing some correlation between some of the coefficients, but this was not done for this paper. The residual is then added to the mean waveform found from the characterization. All realizations were 512 samples in length.
An alternate equivalent procedure could have been used. If the mean is not removed from the characterization, the mean of the wavelets coefficients will reconstruct the mean waveform as was discussed in the previous section. 2. The length of each realization was modeled as in the previous paper [8]. Briefly about the first five records are slightly longer as the gun comes up to speed. Later in the realization the record lengths are about the same with a small random variation. Since the records in the realizations required a length of greater than 512 points, the realizations were padded to the desired length. The padding was from additional synthesized resid-ual records reversed in time. A time reversal was used because by the end of the record, the transient should have decayed to near stationary values, and we want to extend the records with near stationary data. No overlapping or windowing was used at the record boundaries. 3. The 50 individual synthesized records were then assembled into a realization of the simulated gunfire. The synthesized gunfire is shown as the bottom part of Fig. 2. As can be seen the synthesized gunfire and the original gunfire are very much alike. The autospectral densities (power spectra) of both synthesized and original gunfire are shown in Fig. 7. The autospectral densities were estimated using the Welch algorithm with a block size of 1024, a Kaiser-Bessel window, and a 75% overlap. The probability density function (PDF) was estimated for both wave- forms and is plotted as Fig. 8. As can be seen the power spectrum and the PDF are very similar for both waveforms. The tails of the PDF of the gunfire record are matched very well. Notice that the non-Gaussian character of the PDF is preserved in the simulation. A slight surplus of values near zero is observed. 4. If the slight deviation in the PDF is objectionable, a zero memory nonlinear (ZMNL) function can be derived to improve the PDF without significantly changing the spectrum [9]. The slight deviation in the PDF is probably caused by the independence assumption of the wavelet coefficients and leakage in the wavelet reconstruction. Both effects will move the distribution toward a Gaussian distribution. The ZMNL function was obtained from the estimated cumulative distribution function (CDF) of the original waveform and the estimated CDF of the simulation and is plotted as Fig. 9. A linear function is plot-ted for reference. The function is very nearly linear with a little deviation from linear except near the extremes. The simulation was transformed using the ZMNL function, and the resulting modified simulation is shown as Fig. 11. The estimated autospectral density and PDF are shown as Figs 12 and 10. The spectrum is very near the spectrum of the original waveform. A little leakage is seen at the highest frequencies and in some of the spectrum notches. Examples are seen in the notch slightly above 1 kHz and in the frequencies above 2 kHz. The modified PDF is almost indistinguishable from the PDF of the original waveform.
The result is a very credible simulation achieved in a very straightforward manner. The time history of the simulation has the same character as the original record. The autospectrum and PDF of the simulation are very near the autospectrum and PDF of the original.

Conclusions
Wavelets can be an effective tool for analyzing nonstationary random data. In this paper the mean and standard deviation of the wavelet coefficients are used to characterize the waveform. This provided a simple and effective way to characterize and simulate the gunfire record. The method can be easily applied to characterize other nonstationary random events when an ensemble of events is available for analysis.
Additional work needs to be done to clarify the statistical properties of the estimates. Additional work is needed to clarify the role of wavelet packets in characterizing nonstationary random signals. Entropy methods can be used to define an optimum wavelet packet for a deterministic waveform, but an optimum packet for nonstationary random waveforms has not been established.