Digital Spectral Analysis by means of the Method of Averag Modified Periodograms Using Binary-Sign Stochastic Quantization of Signals

The method of averaging modified periodograms is one of the main methods for estimating the power spectral density (PSD). The aim of this work was the development of mathematical and algorithmic support, which can increase the computational efficiency of signals digital spectral analysis by this method. The solution to this problem is based on the use of binary-sign stochastic quantization for converting the analyzed signal into a digital code. A special feature of this quantization is the use of a randomizing uniformly distributed auxiliary signal as a stochastic continuous quantization threshold (threshold function). Taking into account the theory of discrete-event modeling the result of binary-sign quantization is interpreted as a chronological sequence of instantaneous events in which its values change. In accordance with this we have a set of time samples that uniquely determine the result of binary-sign quantization in discrete-time form. Discrete-event modeling made it possible to discretize the process of calculating PSD estimates. As a result, the calculation of PSD estimates was reduced to discrete processing of the cosine and sine Fourier transforms for window functions. These Fourier transforms are calculated analytically based on the applied window functions. The obtained mathematical equations for calculating the PSD estimates practically do not require multiplication operations. The main operations of these equations are addition and subtraction. As a consequence, the time spent on digital spectral analysis of signals is reduced. Numerical experiments have shown that the developed mathematical and algorithmic support allows us to calculate the PSD estimates by the method of averaging modified periodograms with a high frequency resolution and accuracy even for a sufficiently low signal-to-noise ratio. This result is especially important for spectral analysis of broadband signals. The developed software module is a problem-oriented component that can be used as part of metrologically significant software for the operational analysis of complex signals.


Introduction
Spectral analysis is one of the most significant and applied methods for signal investigation. It is used in many areas of engineering and physical sciences. In particular, the spectral estimation of the frequency composition of signals is of practical interest in the process of non-destructive testing and functional diagnostics of technical systems used for various purposes.
In real life signals are exposed to random disturbances and noise interference which leads to their statistical uncertainty. Statistical analysis methods are used to investigate such noisy signals. In this case, the spectral analysis of signals will be associated with the need to obtain a robust estimate of the power spectral density (PSD).
PSD is a continuous function of frequency and determines the average power per unit frequency interval. Estimation of PSD allows us to get an idea of the distribution of the average signal power in the monitored frequency range and to identify within this range the dominant (resonant) frequency components present in its composition. In practice, the PSD is estimated based on the calculation of the direct Fourier transform. Currently, the discrete Fourier transform is used. This is due to the advantages of digital signal processing. First of all, such advantages are the repeatability and reproducibility of digital measuring procedures [1].
With the advent of powerful computer technology and as a result of the development of fast Fourier transform algorithms the periodogram method for estimating the PSD is widely used. According to this method, if the signal is stationary and ergodic in time, then the PSD estimate is calculated by processing single of its centered realization using the finite Fourier transform [2,3]. Centering the signal realization implies preliminary removal of the constant component, if any. However, it was shown in [3] that such a periodogram estimate of the PSD will yield statistically unstable results of spectral analysis. At the same time, an increase in the duration of the processed signal realization does not improve the quality of the PSD estimate. This only leads to an increase in the number of frequencies for which the PSD estimate can be calculated. The smoothing of the PSD estimate can be obtained by averaging the periodogram estimates calculated for individual segments of the signal realization under the assumption that it is stationary. In this case, the observed realization of the signal is divided into segments of finite duration, which form a pseudo-ensemble. Local periodogram estimates are calculated for these segments. Finally, the PSD estimate is calculated by averaging the local periodograms. If the duration of the signal realization is limited in time, then the formation of a pseudo-ensemble with partial overlapping of segments is allowed. This allows you to get more periodogram local estimates for averaging. In order to reduce spectral leakage and reduce distortion of periodogram estimates by strong harmonics, each segment is processed using windowed weighting operations. This reduces the level of side lobes in the spectrum estimate. In addition, the use of windowing functions gives less weight to the values at the ends of the segments. As a result, the overlapping segments are less correlated with each other, which make it possible to obtain a more effective variance reduction for averaged periodogram estimates. This approach to estimating the PSD is known as the method of averaging modified periodograms [2][3][4][5].
The method of averaging modified periodograms leads to a decrease in the variance of the calculated PSD estimates. It also allows you to control the level of the side lobes, which is especially important when the amplitudes of the harmonics in the signal vary greatly. However, this method is a complex signal processing procedure. The need to carry out such complex signal processing leads to the fact that classical algorithms implementing this method in digital form assume the organization of computational processes that require performing a significant number of digital multiplication operations. It is known that it is the multiplication operations that are the most computationally laborious among all arithmetic operations [6,7]. The execution of such algorithms can lead to significant time costs and reduce the efficiency of signal processing, even if the fast Fourier transform is used. The development and use of increasingly powerful computer technology cannot solve this problem completely. This is due to the constant complication of technological processes for monitoring and diagnosing complex systems and objects, which leads to the need to analyze large amounts of data sets. In accordance with the above, the aim of this work was the development of mathematical and algorithmic support, which can increase the computational efficiency of signals digital spectral analysis by the method of averaging modified periodograms. It is necessary to note, that a positive solution to this problem can be obtained

Estimation of the power spectral density based on binary-sign stochastic quantization
Currently widely used digital algorithms for calculating periodogram estimates of PSD are traditionally developed on the basis of the classical approach to the discrete-time representation of signals. According to this approach, analog-todigital signal conversion is the result of performing uniform time sampling, multi-level quantization, and encoding of the digitized samples. Such analog-todigital conversion leads to the necessity of processing multi-bit samples of digital signals in the process of spectral analysis. In practice, in order to reduce time costs, it is necessary to make trade-offs between the quality of analog-to-digital conversion and the computational complexity of processing digital samples. In the extreme case, binary quantization methods are used to convert signals into a digital code [8][9][10][11][12][13][14].
Randomization (deliberate introduction of randomness) of binary quantization makes it possible to obtain a rational relationship between an extremely low number of quantization levels equal to two and the accuracy characteristics of computational algorithms in digital signal processing. In the case of using randomization, we have binary stochastic quantization. A special case of such quantization is binary-sign stochastic quantization. The use of binary-sign stochastic quantization for the statistical analysis of signals is substantiated in [15]. In the process of performing such quantization, a continuous auxiliary random signal is used as a variable quantization threshold (threshold function). The auxiliary signal takes on values within the values range of the analyzed signal and has a uniform distribution. The result of the binary-sign stochastic quantization is: where is a centered signal realization; auxiliary random signal.
The result of binary-sign stochastic quantization can only take on the value "-1" or "+1". The change of these values occurs sequentially at the moments of time . At the same time, the probability of simultaneous occurrence of two or more events at the same time is excluded. Taking this into account, a mathematical model was obtained for representing z(t) in discrete form [16]. It is developed based on the theory of discrete-event modeling. According to this model, it is enough for us to know only the initial value of the binary-sign quantization result z(t 0 ) and the set of time samples within the time interval of signal analysis. Based on this model, the calculation of the periodogram estimate of the PSD using binary-sign stochastic quantization is considered in [17,18]. However, the mathematical equations obtained in [17] do not provide for the use of window functions. In [18], the periodogram estimate of the PSD is calculated using window functions, but its calculation is carried out indirectly through the calculation of the correlation function estimate. This approach to the estimation of PSD requires the formation of initial data using two independent procedures of binary-sign stochastic quantization, which complicates the procedure of spectral analysis.
Taking into account the results obtained by the author in [16][17][18], we consider the application of binary-sign stochastic quantization for estimating PSD by the method of averaging modified periodograms. According to this method, the result of binary-sign stochastic quantization should be represented as a pseudo-ensemble of signal segments. For these segments, local periodogram estimates are calculated, which are averaged to obtain the desired PSD estimate. Figure 1 shows a pseudo-ensemble consisting of M segments for the result of a binary-sign stochastic quantization. The duration of each segment is T = L∆T, where L is the number of overlapping pieces of the segment, and ∆T is the segment shift. For the pseudo-ensemble shown in Figure 1 If we assume that t 0 = 0, then mathematically in time, the segments can be represented as follows: The periodogram estimate for the each weighted segment we will calculate with a frequency resolution f 0 = T -1 at frequencies f k = k f 0 : where w(t) is a window function; W characterizes the average power of the window function and is a normalizing factor.
Devices The averaged periodogram estimate of the PSD will be equal to: We represent (3) in the following form: where Let's introduce the functions: Then for (6) and (7) we get: Taking into account the discrete-event model developed in [16], any segment z(m, t) can be represented in a discrete-time form. As a result of performing binary-sign quantization, we have the initial value z(t 0 ) and a set of time samples , in which ts values change during the analysis time interval T A , where 1 ≤ i ≤ I-1, and . To uniquely identify a segment z(m, t) in time, it is enough to know the set of time samples on the time interval within which this segment is defined, and you also need to know its initial value z(t m-1 ). The result of binary-sign quantization z(t) sequentially takes on the values "-1" and "+1". Therefore, if its initial value z(t 0 ) is known, then the initial value z(t m-1 ) for the segment z(m, t) can be easily determined. Based on this, we introduce the notation The subscripts v(m) and r(m) are integers. Such designations of these indices show their dependence on the segment number. Then for the segment z(m, t) we will have the set of time samples: At time intervals, the boundaries of which are determined by the counts (12), the values for z(t) are equal to "-1" or "+1". Therefore, the integrals in (10) and (11) can be represented as a sum of integrals: (4)  A window w(t) function is a continuous and integrable function. Therefore, the functions h cos ( f, t) and h sin ( f, t) will also be continuous and integrable with respect to the time variable. Consequently, there are functions for H cos ( f, t) and H sin ( f, t) which the differentiability conditions are satisfied: Taking into account (15) and (16), the integrals in (13) and (14) are calculated analytically: where According to (5), when calculating the periodogram estimate of the PSD for the segment z(m, t), the coefficients estimates (17) and (18) must be squared. In this case, we get z 2 (t-1) = 1. It follows that there is no need to know the initial values z (t-1) of the segments. It is enough to have only a set of time samples in which the result of binary-sign stochastic quantization changes its values in the time interval of spectral analysis. This greatly simplifies the computational procedures for estimating the PSD. With this in mind, it is sufficient to calculate the estimates of the following form: Then we get: Equations (4)

Numerical experiments
Based on (4) and (19)-(23), the author has developed algorithmic support and software for calculating the PSD estimate by the method of averaging modified periodograms. The software is made in the form of a specialized software module that can be used as part of a metrologically significant part of the application software for complex signal analysis [19]. When conducting spectral analysis of signals, it should be borne in mind that none of the existing window functions is universal in its purpose. The choice of a specific window function is due to the task and conditions in which the spectral analysis of signals is carried out. Analysis of frequency and metrological characteristics of window functions and recommendations for their application are presented in [20][21][22]. In our case, for the selected window function, it is necessary to have the functions H cos ( f, t) and H sin ( f, t). By definition, these functions are calculated analytically. Therefore, depending on which window functions are supposed to be used, a set of corresponding functions H cos ( f, t) and H sin ( f, t) can be pre-formed and arranged in the form of special collections or libraries of application routines. As an example, we will consider three widely used window functions in practice: rectangular (box car) window, triangular window (Bartlett's) and cosine window. Below for these window functions are the functions  H cos ( f, t) and H sin ( f, t), the values of which for t = 0 and t = T are shown in Table 1.

2) Triangular window (Bartlett's): 3) Cosine window:
Investigations of the computational properties of the developed algorithm for estimating the PSD were carried out on the basis of planning and carrying out numerical experiments. For this purpose, sets of test signal models were used. Each of these models simulated a noisy implementation of a centered signal with a given frequency spectrum structure. The structural composition of the frequency spectrum of the signal model was specified using a combination of summed harmonic components included in its composition. The frequency values of the harmonic components were set in relative units and varied in the range from zero to 0.5. Taking into account the Nyquist-Shannon theorem, they were interpreted as normalized with respect to the doubled value of the upper limit of the frequency band that the spectrum is supposed to occupy. The amplitudes of the harmonic components A n were set in the range from zero to one. Random initial phases φ n were set in the range from -π to +π using a generator of evenly distributed numbers. The sum of the harmonic components was subject to noise. This was achieved by generating white noise that had a zero mean value. The dispersion of white noise was set during the experiment. We considered the signal-to-noise ratio as the ratio of the power of each harmonic component to the noise power. This simulation approach allowed us to investigate the ability to detect harmonic components in noise.
As an example, let's consider the case when the signal model contained nine harmonic components. The numerical values of the frequencies and amplitudes of these harmonic components are presented in Table 2.   Table 1 The values of the functions H cos ( f, t ) and H sin ( f, t ) for the rectangular (box car), triangular (Bartlett's) and cosine windows at t = 0 and t = T, when fk = kf 0 and f 0 = T -1 Window   They were calculated with a resolution of 0.0005 conventional units of the normalized frequency within the entire analysis range from zero to 0.5. For each of the three window functions, the PSD estimates were calculated for one and ten segments. The overlap of the segments was half their length. For the PSD estimates calculated for one segment, we observe the presence of significant fluctuations relative to the true values of the frequency components. At the same time, it is difficult for us to identify the presence of weak harmonic components in the spectrum, for which the signal-to-noise ratio is rather low. In particular, the harmonic component with a frequency of 0.35 and amplitude of 0.05 is masked by noise. It is clearly seen that an increase in the number of processed segments leads to an improvement in the quality of PSD estimates and a decrease in the level of additive noise. All nine frequency components are present on the graphs of PSD estimates calculated for ten segments. The weaker harmonic components are clearly visible and their position in the spectrum corresponds to the tabular data. There are no false spectral lines.

Conclusion
The paper considers the development of mathematical and algorithmic support that allows increasing the computational efficiency of estimating the PSD by the method of averaging modified periodograms in discrete form. This development was carried out on the basis of the use of binarysign stochastic quantization to obtain the digital code of the analyzed signal. A discrete-event model is used to represent the result of binary-sign stochastic quantization in time. This model allowed us to reduce the calculation of PSD estimates to discrete processing of functions that are the result of cosine and sine Fourier transforms for window functions. A set of such functions can be formed beforehand, depending on the window functions used. As a result, we obtained mathematical equations for calculating PSD estimates in discrete form, which do not require performing numerous multiplication operations. These equations became the basis for the development of computational algorithms for estimating the PSD. The main computational operations of the algorithms are the operations of algebraic addition and subtraction. The practical implementation of these algorithms leads to a decrease in computational and time costs in the process of estimating the PSD.
Numerical experiments were carried out on the basis of simulation modeling. They showed that the developed approach to solving the problem posed provides the estimation of PSD by averaging modified periodograms with high frequency resolution. Reliable identification of spectral components is ensured even in additive noise when the signal-to-noise ratio is sufficiently low. This result is especially important for the analysis of broadband signals.
The practical result was the development of a problem-oriented software module for spectral analysis. This module can be used as part of metrologically significant software for the operational analysis of complex signals.