Identification of a Time-Varying Curve in Spectrogram

In this study a process of acquiring a trend of significant spectral coefficients of Photonic Doppler Velocimetry (PDV) data is proposed. The novelty of the paper is the design of a methodology which will allow to find a specific curve describing data of aluminium metal plate acceleration by detonation products of brisant high explosive obtained using PDV in time on the basis of frequency response. The paper combine short time Fourier Transform (STFT), timefrequency varying autoregressive process (TFAR) to specify the description of detonation products from both time and frequency perspectives. We also investigate the identification of a curve describing such behavior of frequency response on time in processed spectrogram.


Introduction
Many engineering applications are focused mostly on analysing the frequency character of a signal and then its time perspective.Because the time character of a signal contains a lot of interesting and important information, it is good to analyse this component especially in connection with the frequency behaviour of the signal.In such a case, the appropriate methodological instrument is the time-frequency (TF) approach predominating in the last decade in many fields of science.It is a useful instrument in natural sciences [1][2][3], engineering [4], [5], biology, medicine [6], or social and economic sciences.
The TF representation of a signal can be obtained via several approaches.The conventional method is the Short Time Fourier Transformation (STFT) which has found usage in a wide range of areas in engineering applications.Similarly, a periodogram or its modification such as the multiple window method using Slepian sequences [7] is also used.We can also perform an estimation via the Wigner-Ville distribution [8], the time-frequency varying Autoregressive Process (TFAR) [9], wavelet analysis (CWT) [10] or alternatively the Modified empirical mode decomposition method [11].
While the periodogram belongs to the group of classic estimators for stationary signals, multiple windows or STFT can be a useful instrument for non-stationary signals.According to Xu et al. [7] multiple window TF distribution was developed to estimate a time-varying spectrum for random non-stationary signals with low bias and variance.Huang et al. [12] uses the Fourier transform and its windowed form and wavelet transform for fringe pattern processing.Adaptive short time Fourier transform proposed by Zhong and Huang [13] is also applicable for TF representation of non-stationary signals.As Jiang and Mahadevan [14] wrote, the advantage of the wavelet analysis is that it can capture the features of a non-stationary signal due to the simultaneous TF decomposition of inputs.The TFAR process is a simplification of the general Auto Regressive Moving Average (ARMA) model.Via estimation of coefficients from the autoregressive process (AR) we can estimate spectral representation on its time-varying form.As Wear et al. [15] wrote, AR methods performed better at short record lengths than the traditional discrete Fourier Transform.In case of multicomponent type of signal Orović and Stanković [4] recommend a class of highly concentrated TF distribution for efficient estimation of frequency.Multicomponent signals were also in focus of Choi and Wiliams [17] working with Cohen's class of TF distribution.They introduced exponential distribution and its windowed form and examines their properties.Boashas and Ristic [16] used specific class of polynomial TF distributions that deals effectively with multicomponent signals and introduced higher-order spectra.
Among the advantages of Fourier transformation and its derivatives we can include low computational complexity, a wide range of software and hardware implementations, and by selecting optimal parameters this method provides sufficient results.The AR process used for estimating signal spectrum representation provides fair results especially in the case of very short signals when STFT tends to fail.For longer signals it provides good results [9].In such cases the variance of insignificant cyclical components which usually takes the character of noise commonly have lower levels than in the case of STFT.This advantage can be useful when we investigate thresholding such as in [18].The TFAR process provides a more complex view in comparison with the simple spectrum estimate in only the frequency domain.It has time and frequency resolution corresponding to the size of the window and the size of the window overlap which must be selected.In such a way it is similar to the STFT method.Unfortunately, the disadvantage of this method is its accuracy which strongly depends on the selection of optimal lag order.Another disadvantage is that there are not many existing implementations; most of them are only on a software level.Considering the methods mentioned above and with respect to the aim of the paper described below we investigated the time-varying autoregressive process.The main reason is based on the fact given in Proakis et al. [9] that AR has good frequency resolution.As shown in Klejmova [19] this fact was confirmed with preliminary analyses on simulated as well as real recorded data from metal plate acceleration by detonation products.
In the context of the stated goal arose the task how to describe the curve shape of the frequency response of explosive material in time.This problem is usually solved using time series modelling mostly done in the time domain.This is widely used in econometrics or social science, but in engineering we can also find its usage.Conventional time domain approaches are regression analyses with deterministic, stochastic or combination trends.The widely used approach in engineering is regression analyses with a linear or polynomial model.Li et al. [20] use a linear model in modelling the electrical conductivity of soil in three dimensions.Fredette et al. [21] uses a cubic polynomial regression model for fitting measured pressure dependent volumes of pumping chambers on a production bushing.Harding et al. [22] uses a generalized additive model in an analysis of long-term trends in climate effects research.Focusing on stochastic modelling, Beveridge-Nelson decomposition provide a definition of the permanent component is often calculated using an ARMA model.Such a model is designed to capture the auto-covariance structure of the first difference of an input series [23].This method assumes that the permanent component following a random walk and the transitory component is stationary with an unconditional mean of zero.We can also apply simplification of ARMA in the form of an AR process of higher lag order.As a consequence, resultant estimate is shorter of number of observations given by the lag with respect to the input time series.In some application areas, filtering techniques operating in the frequency domain are used but the result is usually displayed in time domain.We can use high pass or band pass filters [24].An alternative approach is non-parametric kernel estimation [25] which does not assume knowledge of the data distribution.
In our study we use recorded data of aluminum metal plate acceleration by the detonation products of a brisant high explosive.The data was obtained using a Photonic Doppler Velocimetry (PDV) [26][27][28].Therefore, the data for trend modelling are spectral coefficients which represents material properties at a very rapid load.On the basis of preliminary analysis it seems to be inappropriate to use of deterministic models, because they are not so flexible to capture rapid load.Especially the description of data trends in edge areas can be corrupted by such inflexibility.In the case of a stochastic model, especially in the case of an AR model, the quality depends on lag order causing a loss of data.Because the data character is fitted better with a higher lag order (random walk is insufficient), data loss in the edge area is not welcomed.Applying filtering methods is possible, but we also have here a limitation in the form of parameter selection or pre-filtering data, data loss or less ability to capture a trend in edge boundaries.Such limitations and assumptions can bring forth inaccurate estimations.Therefore, we decided to focus on non-parametric kernel estimation which is quite flexible and does not assume knowledge of distribution.Note, that identifying such curves, i.e. a spectrogram coefficient trend, allows the determination of material properties at a very rapid load.

Paper Contribution and Organization
The aim of our paper is to design a methodological approach leading to the best possible identification of expected spectral peaks and their consequent description.The designed approach will then be applied on a data set obtained from real experiments.Therefore, the paper's aim is focused on: i) identification and ii) specification of vectors containing spectral peaks and their position in TF response of explosive materials.An additional aim is focused on iii) the description of the curve describing such frequency response of explosive materials on explosion in time represented by this vector of peaks.
The novelty presented in the paper is the design of a methodology which will allow to find a specific vector of spectral peaks describing the reaction of an explosive material in time on the basis of frequency response.In the case of an explosion, there is an expected behaviour of the used materials.This expectation can be described by a curve with a specific shape that allows us to determine material properties under difficult experimental conditions (very rapid load, degraded surfaces, ejecta, etc) [29].
The paper is organized as follows: after the introduction, in the section Methodology, we introduce the background of methods for time-frequency analysis.In the section Data Acquisition we describe explosive experiments and real data.In the section Application we present and summarize achieved results.The paper ends with Conclusion.

Methodology
According to the description of the paper's aim, we suggest the following approach.Firstly, we identify a rough region via STFT where a representation of explosion response is found.Secondly, we are focused on specifying the curve describing such a frequency response in the preselected region.Consequently, we suggest using AR for modelling the TF representation because this method is suitable when the frequency resolution is considered.As a supplement, we also use the wavelet transform because the advantage of wavelets is better time resolution.After that, we compare results from the TFAR approach and wavelet approach to specify starting point of the curve describing reaction of material on explosion.To support the decision and selection of the best TF point representing the curve we use test of significance of TF transform on preselected curve points.As the last step, we describe a trend in time of the obtained points (describing frequency response of explosive material) via kernel estimation.

Fourier Transform (FT)
One of the most common methods used for spectrum estimation is the Fourier transform (FT) and its modifications.If an input signal s(n) is a discrete time series, then the Discrete Fourier transform (DFT) is used.A slight modification of this method is called Short Time Fourier Transform (STFT).Then the Fourier transform is calculated using a sliding observation window.The individual spectrum estimations are then sorted in time and can be plotted in a 2D graph.It can be defined as [9] where w(n − m) is the window function.

Autoregressive (AR) Process
Another approach in comparison with the methods mentioned above is the time-frequency varying AR process (TFAR).Comparing with non-parametric FT, TFAR uses parametric approach and creates a model generating an input signal [9].The analyzed signal s(n) is then regarded as the output of a linear filter influenced by white noise w with variance σ 2 w .The autoregressive process can be described by the AR(p) model given by the equation where a i , i = 1, . . .p are the parameters of the autoregressive model of the order p, c is a constant and w n is white noise.
The output spectrum can be described as [9] S( f where H e j2π f T is a linear time variant filter.Thus the spectrum estimation, when we use the AR(p) process, is done according to formula [9] where a i are estimates of the AR(p) parameters, p is the lag order.Several methods for estimating AR(p) model parameters can be used.The most common are the Burg method, Yule-Walker method, unconstrained least-squares method or sequential estimation methods.For the AR process, appropriate selection of lag order p is an important aspect.Selecting a low level order leads to excessive smoothing of the spectrum.Furthermore, if the level of p is selected too high, a non-significant spectral coefficient can arise as a high peak.
For optimal selection, several criteria can be used [9].

Nadaraya-Watson (NW) Kernel Estimator
As written in Introduction, several approaches for trend modelling can be used.With respect to our application and data character we selected non-parametric kernel approach which is quite flexible and does not assume knowledge of the data distribution.Thus the dependency of value Y on value x can be described by the following regression function where m is an unknown function, [30].In our case Y represents the frequency value of the significant spectral coefficient obtained from a spectrogram and x is transformed equidistantly in [0; 1] corresponding to time values in Y .
In our study we are going to investigate the Nadaraya-Watson (NW) estimator [25], [30].Let us assume that If N i=1 K h (x i − x) = 0, then we can define m(x, h) = 0.An optimum bandwidth h estimation can be done via the cross-validation method [25], [31], [32].
For quality evaluation of the estimated trend we can use the mean square error (MSE) or coefficient of determination (R 2 ) calculated using theoretical and estimated spectrum representation for each process.This can be done according to the formula [9] where Y i is the value of the theoretical spectral coefficient, Ȳi is its mean and Y i is its estimate counterpart.In the case of the kernel estimate Y i = m(x, h).The lower the MSE is, the more accurate the estimation is.The closer R 2 is to one, the more accurate the estimation is.

Data Acquisition
For applying the selected methods of TF analysis, we used recorded data of aluminium metal plate acceleration by detonation products of brisant high explosive.The data was obtained using Photonic Doppler Velocimetry (PDV) [26], [27].The PDV setup used is shown in Fig. 1a.A fiber laser with wavelength λ 0 of 1550 nm was used to feed a 3port circulator.Light from the circulator travels towards the probe (collimator or bare fiber end) where it partially exits   and partially reflects back.Light that exits the probe reflects back from the surface of a measured object and reenters the probe.If the object is in motion, the frequency of the reflected light is changed by the Doppler shift f d .This frequency is then combined with the non-Doppler-shifted frequency of the laser source f 0 .The resulting signal with frequency f b equal to the Doppler-shift is captured by the detector and visualized by the oscilloscope.For a more detailed description of PDV, see [28], [29].
In the case of our measurement no amplifier was used between the detector and the oscilloscope.A detailed schematic of the measured object can be seen in Fig. 1b [26].To recalculate the Doppler-shift to the velocity of the moving target v the following formula can be used v = f b λ 0 /2, where λ 0 is the wavelength of the non-Doppler-shifted light (the laser wavelength), f b is the Doppler-shift and v is the resulting velocity.
Two different signals (A and B; Fig. 2) obtained using PDV were selected.Both with a sampling rate of 25 Gs/s and length of 20 µs giving us 500k samples.The maximum available frequency that our data contained was limited by the bandwidth of the used oscilloscope which was approximately 4 GHz.The red arrow in Fig. 2 denotes the time when the Al flyer plate accelerated by the detonation destroyed the PDV fiber probe (level of signal rapidly falls).After this moment the optic fiber is compromised and the measured signal contains unusable noise.

Application
To make the process of describing the trend of significant spectral coefficients more manageable we divided the algorithm into several steps.After data acquisition the procedure was done according to the following steps: 1. Selection of region of interest using STFT.2. Application of STFT, TFAR and supplementary CWT. 3. Combination of AR and STFT results.4. Kernel estimate of trend of significant spectral peaks.

Results: Step 1
The first step was to select a part of the signal containing the required information.The range of the considered signal was defined as follows.The beginning of this time range was taken as the start of the shock wave in the Al plate.The end was defined by the destruction of the optic fiber (rapid fall of the signal).The exact time of destruction is easily visible even in simple time representation.Unfortunately, the exact location of the start of the shock wave is difficult to precisely identify based only on time representation.However, we can perform a preliminary TF analysis of the signal.Therefore, based on this TF analysis we can select an appropriate sample range.For this purpose the STFT method is sufficient (low computational requirements, easy implementation).
The preliminary STFT analysis results can be seen in Fig. 3 where the area of interest is highlighted by a circle.Based on the detailed analysis, we selected, for signal A, samples within a range of 5.00 to 6.80 µs which corresponds to 45k samples, and for signal B a range of 5.00 to 6.88 µs which corresponds to 47k samples.To have access to the low frequency slope we avoided any filtering of data.

Results: Step 2
In the second step we repeatedly applied STFT on the preselected range described above.Establishment of parameters was motivated by the intention of keeping a balance between time and frequency resolution.Therefore, we used a window length of 1024 samples with an overlap of 900 samples.Due to wide overlap, the Hann filtering window was used to suppress the influence of earlier signal samples.Results are presented in Fig. 4a, b.
In the next step we obtained TF representation using the AR process.Its application was also done on the preselected range described above.Motivation for parameter selection was the same as for STFT, therefore, we chose a window length of 1024 samples with an overlap of 900 samples and the Hann filtering window.Considering the sample size of the signal and the size of the window we chose the Yule-Walker method [33], [34] for AR coefficients estimation which is more suitable for long signals.Selection of lag order was done separately for each window (i.e.1024 samples).To determine optimal lag order we used the Akaike information criterion [19] providing good results in the case of a similar signal.Results are presented in Fig. 4c, d.
As we can see in Fig. 4b,d the results show a rising edge.Based on the knowledge from energetic material research, this rising edge represents the initial moment of detonation.However, under some circumstances, it might not be captured.To support its existence we also provided CWT   for signal B which confirmed the rising edge.Therefore, in such a case it is not edge effect of the estimation.

Results: Step 3
To obtain the best possible TF representation we combined results from the STFT and TFAR approach.We considered two approaches, multiplication and correlation [9] of amplitude part of TF transform for each signal to noise suppression.In the case of signal A (Fig. 5a) lower signal noise to ratio (SNR) was achieved using correlation while for the signal B (Fig. 5b) by simple multiplication.Because the main focus was on the amplitude part of the spectra we multiplied only the amplitude part of STFT with TFAR.Doing this we managed to smooth the scatter of background noise.Both methods (STFT with TFAR) used the same window length and the same overlap meaning that the time and frequency resolution of both methods were compatible.
In step 3 we used the signal gained in step 2. We firstly designed thresholding according to the significant values of the spectrogram.The threshold value was established as follows: with respect to the frequency value we took the highest spectrogram value ±∆ to specify close surrounding area.According to empirical results we established ±∆ = 7.5 %.
Using this thresholding we specified a mask which contained the region of our interest (region of required signal).After applying the mask on the signal from the previous step we were able to mark the position of spectral peak belonging to the area defined by the mask.By implementing this step we gained a vector of peak positions suitable for kernel analysis.

Results: Step 4
The last step was focused on kernel estimation of the trend of significant spectra peaks.Identification of such trend allows better determination of materials properties at very rapid load than simple maxims connection.We firstly created vectors containing the position of detected peaks (step 3).The length of vectors was approximately 200 points.For both signals we used the NW estimator for a kernel of the order ν = 0; k = 2, smoothness µ = 2 and bandwidth h = 0.02 [30], [31].The value of bandwidth was optimized via the generalized cross-validation method [35].The NW estimator was selected for providing good results in case of PDV signal, for the details of kernel estimator application see [25].
In both signals, Fig. 6c and f, we can identify the same shape of the curve.Comparing them, we can see a rising edge of the curve at the start of signal B (Fig. 6f).This shape of curve path corresponds with expectations of material behavior during explosion.It is missing in signal A and is probably caused and influenced by the explosive event.Therefore, the start of signal A was identified after this part.In both signal cases we can identify a three-level decrease with osculation followed by its increase.The general tendency has a similar parabolic shape in both signals.
A detailed analysis of Fig. 6c and f reveals the following facts.Signal A (Fig. 6c): the graphical representation reveals a short time decrease followed by a short time increase and stagnation of the signal at a higher level.The second and third part copies this tendency, but both parts are moved to the lower level having a stair shape.The last part takes a parabolic shape with the second half of the signals.
In the case of signal B (Fig. 6f) we can identify three similar parts of decrease.In comparison with signal A, the first part is rather stagnation with a very slow decrease followed by a sharper decline divided by structural breaks.Generally, the dynamics of the curve is smaller.The final part of the curve for signal B, taking a parabolic shape, is similar to the curve for signal A (Fig. 6c).

Comparison of Curve Shape Modelling
For modelling the trend, we used vectors containing position of detected peaks for both signals described above.To support decision beneficial to kernel estimate we investigated two parametric and one non-parametric model.Namely the deterministic polynomial model (Fig. 6b,e), stochastic AR model (Fig. 6a,d) and non-parametric NW kernel estimate (Fig. 6c,f).In Tab. 1 you can see the evaluation of model fit via MSE and R 2 with optimized parameters.
The comparison of estimated curve shapes (Fig. 6) and measured quality values (Tab. 1) reveal, that for signal A the best approach is the NW estimate while for signal B it seems to be the AR process.However, AR appears problematic when focusing on capturing the rising edge because in such a case it was not captured.This is caused by a loss of data in the initial part of the signal -area of omission depends on lag order of the AR process corresponding to the signal character.Another problem with AR is in the case of missing points (non equidistant) approximation the resultant curve looks under-smoothed.Making the sentence nonsense (missing points as adjective -describing the noun approximation.If I understand correctly, perhaps: Another problem with AR is in the case of missing points (non equidistant); the approximation of the resultant curve looks under-smoothed.We investigated lag order p = 1 − 20; optimal orders were selected using information criterion (Tab.1).For this reason, we prefer the NW estimate (Fig. 6c and f) as the best even when it has a little bit worse measured quality values for signal B, because the rising edge was clearly captured.The worst curve fit estimation was achieved via the polynomial model (Fig. 6b and e).We investigated order q = 2 − 30.As the results show, in such an approach, the estimated fit is over-smoothed.Also, its ability to describe a rising edge is worse compared to the NW estimate.The acquired measurements were the worst among all.To sum up, the kernel estimate provides the best curve shape estimate for such data.

Conclusion
The aim of the presented paper was to design a methodological approach leading to the best possible expected curve identification of significant spectral coefficients for timefrequency representation of Photonic Doppler Velocimetry.We proposed the following steps of approach.Firstly, we recommended FFT to preselect an area of interest.Consequently, we estimated TF representation via FFT and the time-varying AR process.By multiplication or correlation of both transformations we acquired better frequency resolution.Additionally, to support the results, we proposed applying CWT.Then, after thresholding, we established a mask to obtain significant spectral coefficients of processed timefrequency representation of the detonation product.Finally, we recommended identifying the curve describing the behavior of frequency response in the processed spectrogram by using a non-parametric kernel estimator.

( a )
Spectrogram of signal A using STFT.(b) Spectrogram of signal B using STFT.(c) Spectrogram of signal A using AR.(d) Spectrogram of signal B using AR.

( a )
Combined spectrogram of signal A. (b) Combined spectrogram of signal B.
Modelling of the curve fit (signal A: a-c; signal B: d-f).Note: the velocity was calculated from frequency using v = f λ 0 /2, where λ 0 is the the laser wavelength.Evaluation of the model fit via MSE and R 2 .