Analysis of the Precision of Bispectral Estimation Methods

Author(s)
Wenbing Wu

ABSTRACT

In order to evaluate the accuracy of bispectral estimation method, signals of cosine function were adopted. Because the cosine signals’ three order moment spectrum and three order cumulant spectrum are zero, Non zero part of the bispectrum estimated by no matter which method is the estimation error. Through the comparison of three kinds of estimation methods: the direct method, indirect method and AR parameter method, errors of various estimation methods were obtained, then changing the values of different parameters in all methods, and observing the bispectral error values changes with the parameters, so as to provide a basis for the various bispectrum estimation method and the selection.

In order to evaluate the accuracy of bispectral estimation method, signals of cosine function were adopted. Because the cosine signals’ three order moment spectrum and three order cumulant spectrum are zero, Non zero part of the bispectrum estimated by no matter which method is the estimation error. Through the comparison of three kinds of estimation methods: the direct method, indirect method and AR parameter method, errors of various estimation methods were obtained, then changing the values of different parameters in all methods, and observing the bispectral error values changes with the parameters, so as to provide a basis for the various bispectrum estimation method and the selection.

1. Preface

Higher-order cumulant can suppress the effects of gaussian background noise automatically, because of this, high order cumulants have received more and more attention and have become a very useful tool in signal processing. Double spectrum obtained from third order cumulant contains information asymmetric and nonlinear signal, can be used to describe the nonlinear phase coupling, quadratic phase coupling, in particular, has been widely applied in the fault diagnosis.

The contemporary history of the spectrum estimation began with the breakthrough of Tuckey in 1949. Blackman-tuckey gave a method to obtain the power spectrum estimation from the sample data sequence using Wiener correlation method—BT method.

The occurrence of FFT (fast fourier transform) in 1965 produced the periodogram; the fatal weakness of BT and periodogram is the limitation of frequency resolution. In order to overcome this shortcoming, Burg was inspired by the linear prediction method in the study of earthquake application in 1967, and derived the maximum entropy spectrum estimation method. E.p. arzen formally proposed the AR spectrum estimation method in 1968; Since then, many high resolution spectral estimation methods have been developed in the past decade, which is called the modern spectral estimation method, and the BT method and the periodogram are called the traditional spectral estimation method.

Either way, each spectrum estimation technique can be considered a model method. Specifically, based on the prior knowledge of the process, a model of approximate actual process is established. Secondly, the model parameters of the hypothesis are estimated by observation data or autocorrelation function, finally to make a spectrum estimation [1] [2] .

The commonly used models include periodogram and BT (sine harmonic summation model), AR (Auto Regression Model), MA (the moving average model), ARMA (the Auto Regression Integrated with Moving), Prony, maximum likelihood, etc. The variation of performance between various estimates is caused by the difference between the hypothetical model and the actual process. Although different models can produce similar results, some models may need fewer parameters.

It is general to estimate the high order cumulant spectrum or high order moment spectra using finite length data. There are two main methods:

1) Conventional (Fourier) method;

2) Parameter methods: AR, MA, ARMA and other methods.

This paper discusses the conventional method and the double spectrum estimation effect of AR parameter method. Conventional methods can be divided into the following two categories:

1) Indirect method: based on the definition of Formula (6) or Formula (7) to estimate;

2) Direct method: based on the definition of.

Formula (8) to estimate.

Due to the inevitable error of various estimation methods, this paper will analyze the error of these methods by the double spectrum estimation of cosine signal [3] [4] .

2. High Order Cumulant

If {x(n)} is the zero mean k order stationary random process, the k order cumulant ${c}_{kx}\left({\tau}_{1},{\tau}_{2},\cdots ,{\tau}_{k-1}\right)$ of the process is defined as the k order joint cumulant of random variable $\left\{x\left(n\right),x\left(n+{\tau}_{1}\right),\cdots ,x\left(n+{\tau}_{k-1}\right)\right\}$ , namely

${c}_{kx}\left({\tau}_{1},{\tau}_{2},\cdots ,{\tau}_{k-1}\right)=cum\left\{x\left(n\right),x\left(n+{\tau}_{1}\right),x\left(n+{\tau}_{2}\right),\cdots ,x\left(n+{\tau}_{k-1}\right)\right\}$ (1)

The k order moment ${m}_{kx}\left({\tau}_{1},{\tau}_{2},\cdots ,{\tau}_{k-1}\right)$ of the process is defined as the k order joint moment of the random variable $\left\{x\left(n\right),x\left(n+{\tau}_{1}\right),\cdots ,x\left(n+{\tau}_{k-1}\right)\right\}$ , namely

${m}_{kx}\left({\tau}_{1},{\tau}_{2},\cdots ,{\tau}_{k-1}\right)=mom\left\{x\left(n\right),x\left(n+{\tau}_{1}\right),\cdots ,x\left(n+{\tau}_{k-1}\right)\right\}$ (2)

Here, mom() represents the joint moment, the third order cumulant is:

${c}_{3x}\left({\tau}_{1},{\tau}_{2}\right)=E\left\{x\left(n\right)x\left(n+{\tau}_{1}\right)x\left(n+{\tau}_{2}\right)\right\}$ (3)

The k-order cumulant spectrum is defined as k − 1 dimension Fourier transform of k order cumulant.

Namely

${C}_{kx}\left({\omega}_{1},\cdots ,{\omega}_{k-1}\right)={\displaystyle \underset{{\tau}_{1}=-\infty}{\overset{\infty}{\sum}}\cdots}{\displaystyle \underset{{\tau}_{k-1}=-\infty}{\overset{\infty}{\sum}}{c}_{kx}\left({\tau}_{1},\cdots ,{\tau}_{k-1}\right)}\mathrm{exp}\left[-j{\displaystyle \underset{i=1}{\overset{k-1}{\sum}}{\omega}_{i}{\tau}_{i}}\right]$ (4)

k order moment spectrum is defined as k − 1 dimension Fourier transform of k order moment.

Namely

${M}_{kx}\left({\omega}_{1},\cdots ,{\omega}_{k-1}\right)={\displaystyle \underset{{\tau}_{1}=-\infty}{\overset{\infty}{\sum}}\cdots}{\displaystyle \underset{{\tau}_{k-1}=-\infty}{\overset{\infty}{\sum}}{m}_{kx}\left({\tau}_{1},\cdots ,{\tau}_{k-1}\right)}\mathrm{exp}\left[-j{\displaystyle \underset{i=1}{\overset{k-1}{\sum}}{\omega}_{i}{\tau}_{i}}\right]$ (5)

The most commonly used high order spectrum is the third-order spectrum (also known as the double spectrum), and the estimation formula as follows:

${B}_{c}\left({\omega}_{1},{\omega}_{2}\right)={\displaystyle \underset{{\tau}_{1}=-\infty}{\overset{\infty}{\sum}}{\displaystyle \underset{{\tau}_{2}=-\infty}{\overset{\infty}{\sum}}{c}_{3x}\left({\tau}_{1},{\tau}_{2}\right){\text{e}}^{-j\left({\omega}_{1}{\tau}_{1}+{\omega}_{2}{\tau}_{2}\right)}}}$ (6)

${B}_{m}\left({\omega}_{1},{\omega}_{2}\right)={\displaystyle \underset{{\tau}_{1}=-\infty}{\overset{\infty}{\sum}}{\displaystyle \underset{{\tau}_{2}=-\infty}{\overset{\infty}{\sum}}{m}_{3x}\left({\tau}_{1},{\tau}_{2}\right){\text{e}}^{-j\left({\omega}_{1}{\tau}_{1}+{\omega}_{2}{\tau}_{2}\right)}}}$ (7)

For periodic power signal, the high order moment spectrum is generally estimated, order

$X\left(k\right)={\displaystyle \underset{n=0}{\overset{N-1}{\sum}}x\left(n\right){\text{e}}^{-i\frac{\text{2\pi}}{N}nk}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}0\le k\le N-1$

The Formula (7) can be used to calculate the double spectrum by the following formula:

${B}_{x}\left({k}_{1},{k}_{2}\right)=\frac{1}{N}X\left({k}_{1}\right)X\left({k}_{2}\right){X}^{\ast}\left({k}_{1}+{k}_{2}\right)$ (8)

3. Double Spectral Estimated Methods

1) AR parameter method

Assuming y_{1}(t) is the signal of the actual output signal y(t) of the system, the random vibration signal of the system output is caused by the non-gaussian white noise a(t), whose mean value is equal to zero, and establishes the AR model:

${y}_{1}\left(t\right)+{\displaystyle \underset{i=1}{\overset{p}{\sum}}{\psi}_{i}{y}_{1}\left(t-i\right)}=a\left(t\right),\left(t=1,2,\cdots ,N\right)$

${\psi}_{i}\left(i=1,2,\cdots ,p\right)$ is the autoregressive coefficient, p is the order number of the autoregressive model.AR bispectrum expression:

${B}^{AR}\left({\omega}_{1},{\omega}_{2}\right)={\gamma}_{a.3}H\left({\omega}_{1}\right)H\left({\omega}_{2}\right){H}^{*}\left({\omega}_{1}+{\omega}_{2}\right)$

$H\left(\omega \right)$ is the deliver function,

$H\left(\omega \right)=\frac{1}{1+{\displaystyle \underset{i=1}{\overset{p}{\sum}}{\psi}_{i}{\text{e}}^{-ji\omega}}}$

For y_{1}(t), the parameter method of AR model is used to estimate the model coefficient
$\psi =\left[\beta \right]$ and AR double spectrum amplitude expression

$\left|{B}^{AR}\left({\omega}_{1},{\omega}_{2}\right)\right|=\frac{\left|{\gamma}_{\alpha .3}\right|}{\left|1+{\displaystyle \underset{i=1}{\overset{p}{\sum}}{\beta}_{i}{\text{e}}^{-ji\left({\omega}_{1}+{\omega}_{2}\right)}}\right|{\displaystyle \underset{k=1}{\overset{2}{\prod}}\left|1+{\displaystyle \underset{i=1}{\overset{p}{\sum}}{\beta}_{i}{\text{e}}^{-ji{\omega}_{k}}}\right|}}$

To estimate the model coefficient, the model must be ordered first. Singular value decomposition (SVD) is a model order method, the normalized Frobenius norm method is a kind of order method based on singular value decomposition. The normalized Frobenius norm is defined as:

$\mu \left(p\right)=\frac{{\left[{\lambda}_{1}^{2}+{\lambda}_{2}^{2}+\cdots +{\lambda}_{p}^{2}\right]}^{1/2}}{{\left[{\lambda}_{1}^{2}+{\lambda}_{2}^{2}+\cdots +{\lambda}_{n}^{2}\right]}^{1/2}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}p=1,2,\cdots ,n$ (9)

the denominator in the formula is Frobenius norm.

Obviously, $\mu \left(p\right)\le 1$ .

That $\alpha \approx 1$ is a threshold value is close to 1, when $\mu \left(p\right)\ge \alpha $ , it thinks the first p singular value dominated, which choose to satisfy this condition of minimum p value as matrix A (AR) order of effective rank estimation results. Typically, the minimum p value corresponding to $\alpha =0.995$ , namely, $\mu \left(p\right)\ge 0.995$ is the optimal order of the AR model.

2) The direct method

For deterministic signals, the step of the direct method of double spectrum estimation is:

Step 1. Divide the data into K segments. Each segment contains M sample points, i.e., N = K * M, and the date of every piece is averaged out. If it’s a deterministic signal, then only one record (N = M) is used. If you want to get the length that the FFT algorithm fits, a zero can be added to each piece of data.

Step 2. Suppose $\left\{{x}^{i}\left(q\right)\right\},q=0,1,\cdots ,M-1$ is the i data, and the DFT coefficient is calculated.

${X}^{i}\left(\lambda \right)={\displaystyle \underset{n=0}{\overset{N-1}{\sum}}{x}^{i}\left(q\right){\text{e}}^{-i\frac{\text{2\pi}}{M}q\lambda}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}0\le \lambda \le M-1,\text{\hspace{0.17em}}1\le i\le K$

Step 3. The final K segment data is averaged to get a double spectrum of the given data.

${B}_{x}\left({\lambda}_{1},{\lambda}_{2}\right)=\frac{1}{K}X\left({\lambda}_{1}\right)X\left({\lambda}_{2}\right){X}^{\ast}\left({\lambda}_{1}+{\lambda}_{2}\right)$ (10)

3) The indirect method

For deterministic signals, the two-spectral estimation indirect method is:

Step 1. Divide the data into K segments, and each segment has M sample points, i.e., N = K ∙ M;

Step 2. the date of every piece is averaged out;

Step 3. Suppose $\left\{{x}^{i}\left(q\right)\right\},q=0,1,\cdots ,M-1$ is the i data, and the third order moment is calculated

${M}_{3}\left({\tau}_{1},{\tau}_{2}\right)=\frac{1}{M}{X}^{i}\left(q\right){X}^{i}\left(q+{\tau}_{1}\right){X}^{\ast}\left(q+{\tau}_{2}\right)$ (11)

Step 4. Double spectrum estimation is performed according to Equation (7).

4. Two Spectra of Cosine Signal

Supposecosine signal $x\left(n\right)=\mathrm{cos}\omega t$ , By third-order cumulants calculation formula

${c}_{3x}\left({\tau}_{1},{\tau}_{2}\right)=E\left\{x\left(n\right)x\left(n+{\tau}_{1}\right)x\left(n+{\tau}_{2}\right)\right\}$ (12)

can be known, cosine signal of third-order cumulants is 0, so whatever according to type (7) or type (8), the double spectrum are calculated for 0. See literature (5) and (6).

5. Double Spectrum Estimation Error Analysis of Cosine Signal

It can be seen from the above analysis that the cosine signal is zero for both its third-order cumulant or the third-order moment spectrum. Based on the above three estimation methods. Due to the known cosine signal bispectrum is zero, so no matter what kind of method to estimate the cosine signal power spectrum, as long as there is not zero, it can be regarded as the estimate method of error.

First, the AR parameter method is used to estimate.

This article takes the cosine signal y(t) = cos(2 * PI * f_{0} * t), respectively in f_{0} 1 HZ, 5 HZ and 50 HZ, the sampling frequency is 1000 HZ, this kind of selection is mainly in random, and considering the interval of frequencies simultaneously, the sampling time is 1 second, adopt the AR parameters method, direct method and indirect method for double spectrum estimation, respectively is shown in Figures 1-3. Can be seen from the above image, when f_{0} to 1 HZ, estimate the double spectrum distribution minimum non-zero part, but the value is bigger. Under the condition of the sampling frequency must, therefore, with the increase of signal frequency, whatever estimation method, the double spectrum error estimate its distribution range is increased, but the value is the trend of decrease. In keeping the signal frequency of 1 HZ again under the condition of not changing, sampling frequency, respectively, using 100 HZ and 100 HZ and 2000 HZ, three methods to estimate the bispectrum are shown in Figures 4-6. You can see from the picture, under the condition of the signal frequency is fixed, when the sampling frequency from 100 HZ to 100 HZ, parameter method and direct method of the estimation results are present a bigger change, estimation precision is improved more, and had less indirect rule to improve. When the

Figure 1. The double spectrum of sampling frequency 1000 HZ and AR parameter estimation. (a) Signal frequency 1 HZ; (b) Signal frequency 5 HZ; (c) Signal frequency 50 HZ.

Figure 2. The double spectrum of sampling frequency 1000 HZ, direct method estimation. (a) Signal frequency 1 HZ; (b) Signal frequency 5 HZ; (c) Signal frequency 50 HZ.

Figure 3. Double spectrum of sampling frequency 1000 HZ, indirect method estimation. (a) Signal frequency 1 HZ; (b) Signal frequency 5 HZ; (c) Signal frequency 50 HZ.

Figure 4. The double spectrum of signal frequency 1 HZ and parameter estimation. (a) Sampling frequency 100 HZ; (b) Sampling frequency 500 HZ; (c) Sampling frequency 2000 HZ.

Figure 5. Double spectrum of signal frequency 1 HZ and direct method estimation. (a) Sampling frequency 100 HZ; (b) Sampling frequency 500 HZ; (c) Sampling frequency 2000 HZ.

Figure 6. The double spectrum of signal frequency 1 HZ and direct method estimation. (a) Sampling frequency 100 HZ; (b) Sampling frequency 500 HZ; (c) Sampling frequency 2000 HZ.

change from 500 HZ to 2000 HZ, the results of the three methods changed very little. Figure 7 and Figure 8 show the results of the number of steps in the direct and indirect methods by changing the number of steps in the direct and indirect methods by changing the frequency of the signal at the fixed level of 1 HZ and the sampling frequency fixed bit 1000 HZ. In the figure, the accuracy of the two estimation methods has been improved greatly with the number of data segments from 64 to 16. Comparing Figure 2(a) and Figure 3(a), both represent the signal frequency of 1 HZ, the result of the data block is equal to the number 8 section, do not Figure 7(a) and Figure 8(a), its data segment is equal to the number 4, the two ways to estimate the result has little difference. In the AR method, then by changing the type (9) normalized Frobenius norm values, it is concluded that the bispectrum as follows: Figure 9(a) and Figure 9(b) of the image is in the signal frequency of 50 HZ, under the condition of the sampling frequency is 1000 HZ. Compared with that Figure 1(c), Figure 9(a) and,

Figure 7. The double spectrum of sampling frequency 1000 HZ and signal frequency 1 HZ direct method estimation. (a) Section 64 of the data; (b) Section 16 of the data; (c) Section 4 of the data.

Figure 8. The double spectrum of sampling frequency 1000 HZ and signal frequency 1 HZ direct method estimation. (a) Section 64 of the data; (b) Section 16 of the data; (c) Section 4 of the data.

Figure 9. Frobenius norm value to 0.8 and 0.999. (a) Frobenius norm value to 0.8; (b) Frobenius norm value equal to 0.999.

Figure 9(b) of Frobenius norm value equal to 0.8 and 0.999 respectively, Figure 1(c) value equal to 0.995, by contrast, you can see that in the AR method, Frobenius norm estimation precision of the numerical change can also lead to the change, when the Frobenius norm values change from 0.8 to 0.995, the estimated accuracy change is bigger, because Figure 9 compared with Figure 1(c), double spectrum distribution significantly more of a non-zero value, and when the Frobenius norm values change from 0.995 to 0.999, the change is not obvious, which indirectly proves that the Frobenius norm value to 0.995 as the rationality of the threshold. As can be seen from the above analysis results, the estimation effect of parametric method is slightly better than that of direct and indirect methods.

6. Brief Summary

Based on the spectrum and the theoretical value of zero cosine signal bispectrum estimation, three kinds of methods for analysis of the various led to the error of estimation methods, and by changing the parameters, analysis of different estimation methods in different cases, the change of the error, to the precision of estimation methods provide a judgment standard [5] [6] .

Cite this paper

Wu, W. (2019) Analysis of the Precision of Bispectral Estimation Methods.*International Journal of Modern Nonlinear Theory and Application*, **8**, 62-71. doi: 10.4236/ijmnta.2019.83005.

Wu, W. (2019) Analysis of the Precision of Bispectral Estimation Methods.

References

[1] Zhang, X.D. (2002) Modern Signal Processing. Tsinghua University Press, Beijing, 263-281.

[2] Wu, W.B. and Huang, Y.J. (2011) Based on the Fault Diagnosis of the Pressure-Reducing Valve of Bicoherent Spectrum. Computer Measurement and Control, 19, 2413-2416.

[3] Shao, R.P., Huang, X.N. and Liu, H.Y. (2008) Fault Detection and Diagnosis of Gear System Based on High Order Cumulant. Journal of Mechanical Engineering, 44, 161-168.

[4] Kocur, D. and Stanko, R. (2000) Order Bispectrum: A New Tool for Reciprocated Machine Condition Monitoring. Mechanical Systems & Signal Processing, 14, 871-890.

https://doi.org/10.1006/mssp.2000.1307

[5] Zhang, Y. and Wang, S.X. (1998) Analysis of the Method of Slice Spectrum Analysis of Nonlinear Phasecoupling. Electronic Journal, 26, 104-109.

[6] Yang, J.T., Zhou, L.H. and Chen, J.J. (1999) Higher Order Spectra and Its Application in Machinery Fault Diagnosis. Transactions of Tianjin University, 5, 190-194.

[1] Zhang, X.D. (2002) Modern Signal Processing. Tsinghua University Press, Beijing, 263-281.

[2] Wu, W.B. and Huang, Y.J. (2011) Based on the Fault Diagnosis of the Pressure-Reducing Valve of Bicoherent Spectrum. Computer Measurement and Control, 19, 2413-2416.

[3] Shao, R.P., Huang, X.N. and Liu, H.Y. (2008) Fault Detection and Diagnosis of Gear System Based on High Order Cumulant. Journal of Mechanical Engineering, 44, 161-168.

[4] Kocur, D. and Stanko, R. (2000) Order Bispectrum: A New Tool for Reciprocated Machine Condition Monitoring. Mechanical Systems & Signal Processing, 14, 871-890.

https://doi.org/10.1006/mssp.2000.1307

[5] Zhang, Y. and Wang, S.X. (1998) Analysis of the Method of Slice Spectrum Analysis of Nonlinear Phasecoupling. Electronic Journal, 26, 104-109.

[6] Yang, J.T., Zhou, L.H. and Chen, J.J. (1999) Higher Order Spectra and Its Application in Machinery Fault Diagnosis. Transactions of Tianjin University, 5, 190-194.