Auto-regressive moving average parameter estimation for 1/f process under colored Gaussian noise background

Abstract

Current algorithms for estimating auto-regressive moving average parameters of transistor 1/f process are usually under noiseless background. Transistor noises are measured by a non-destructive cross-spectrum measurement technique, with transistor noise first passing through dual-channel ultra-low noise amplifiers, then inputting the weak signals into data acquisition card. The data acquisition card collects the voltage signals and outputs the amplified noise for further analysis. According to our studies, the output transistor 1/f noise can be characterized more accurately as non-Gaussian α-stable distribution rather than Gaussian distribution. We define and consistently estimate the samples normalized cross-correlations of linear SαS processes, and propose a samples normalized cross-correlations-based auto-regressive moving average parameter estimation method effective in noisy environments. Simulation results of auto-regressive moving average parameter estimation exhibit good performance.

Introduction

A problem attracting much attention in the past decades belongs to parameter estimation of a model from noisy observations which is built by assuming that its output was generated by exciting a linear system with unobservable independent and identically distributed (i.i.d.) non-Gaussian noises. Typical of such models is auto-regressive moving average ^{1 ,2} (ARMA) model. Existing research on this subject is mainly based on the Gaussian noise assumption. Gaussian assumption brings the conveniences of theoretical analysis, for the existence of second and higher order statistics. Unfortunately, for certain realistic conditions, the Gaussian assumption is proved to be inadequate. Many types of non-Gaussian noises in real world manifest impulsive characteristics in the time domain while exhibiting heavy-tail behavior reflected in the probability density. ^{3 ,4} According to our research, in the case of Zener diode, its output 1/f noise can be characterized more accurately as non-Gaussian α-stable process ⁵ rather than Gaussian process. Under unavoidable additive noise environments, parameter estimation of a model from noisy observations which is built by assuming that its output was generated by exciting a linear system with unobservable i.i.d. non-Gaussian α-stable noise. With the additive colored Gaussian noise (CGN), ARMA parameter estimation for 1/f noise conforming to α-stable distribution becomes a more challenging problem.

The low frequency 1/f noises are manifested as voltage fluctuation at both ends of the transistor. Traditional method of accelerating electrical aging technique has caused a certain degree of internal damage to transistors and affects the reliability. Therefore, the noise analysis method is proposed to measure noise of the transistor under regular working conditions. This method possesses advantages of high efficiency and non-destructive to the test samples; thus, it becomes the universal noise measurement method.

The main difficulty to ARMA parameter estimation for α-stable distribution lies in the infinite variance. Fortunately, the fractional lower order moments (FLOM) ⁶ exist in α-stable process; thus, available parameter estimation methods of α-stable process are basically based on FLOM and its derived fractional lower order statistics ⁷ (FLOS), e.g. fractional lower order covariance (FLOC) and covariation.

Although FLOS-based methods can effectively solve the ARMA parameter estimation problem of α-stable process, their inherent non-linear characteristics make it difficult for suppression of additive CGN for α-stable process. To solve such problem, we propose a cross-spectrum method effective in noisy environments based on samples normalized cross-correlations (SNCC) function, and comprehensive studies will be carried out on the definition of SNCC function, the ARMA model parameter estimation method of α-stable noise and its performance. The simulation results confirm the effectiveness of the proposed method.

In the next section, the observed data model of linear SαS process is given. Then, we propose a new function of SNCC, the definition of SNCC, its operation properties, and SNCC-based Yule–Walker equation of ARMA model. Next, an SNCC-based parameter estimation algorithm of ARMA model is developed, and auto-regressive (AR) and moving average (MA) coefficients estimation methods are given. Then, experiments are carried out. Finally, some conclusions are drawn.

Linear SαS process in the observed data

The observed data are the measured weak voltage and additive system interference noise. The former is reflected as the level of the AC voltage at both ends of a diode, including 1/f noise and white noise. In the light of our research, taking Zener diodes for example, the experimental results indicated that diode 1/f noise in time domain is impulsive. By the examinations on infinite variance test, the consistency of the model parameter estimation results, the statistical characteristics of 1/f noise signal, and the combination of hypothesis test and real-data experiment, we can draw the conclusion that 1/f noise can be modeled as α-stable process. ⁵

Since α-stable process can be expressed by an ARMA process driven by $S\alpha S$ , the single channel observed time series is modeled as the output of an ARMA system excited by unobservable $S\alpha S$ input, and corrupted by additive CGN. We describe ARMA ( $p,q$ ) process of an $S\alpha S$ process $X\left(n\right)$ ⁸

$$X\left(n\right)=-{\displaystyle \sum _{k=1}^{p}a\left(k\right)}X\left(n-k\right)+{\displaystyle \sum _{j=0}^{q}b\left(j\right)}U\left(n-j\right)$$

(1)

given only by the observed data $Y\left(n\right)$ , presenting 1/f spectrum characteristics. We assume that input $U\left(n\right)$ is i.i.d. $S\alpha S$ with characteristic exponent $\alpha \left(1\prec \alpha \prec 2\right)$ , where A $\prec$ B means that A is less than B. $a\left(k\right)$ is the AR coefficients, $b\left(j\right)$ is the MA coefficients, $X\left(n\right)$ is the noiseless system output, $p$ is the order of the AR model, $q$ is the order of the MA model. The transfer function of this system is

$$H\left(z\right)=\frac{B\left(z\right)}{A\left(z\right)}=\frac{1+b\left(1\right)z^{-1}+\cdots +b\left(q\right)z^{-q}}{1+a\left(1\right)z^{-1}+\cdots +a\left(p\right)z^{-p}}$$

(2)

To measure voltage signal precisely, non-destructive noise analysis method is adopted, which is to filter the voltage signals through ultra-low noise amplifier, and then the weak noise is amplified and input into data acquisition card (DAQ). The DAQ collects the voltage signal and outputs the amplified signal for further analysis. Amplifier will introduce background noise into the output signal, which obviously leads to decline of spectrum estimation accuracy. In order to solve this problem, a method of cross-spectrum measurement technique was proposed by Dai ⁹ to effectively eliminate the influence of noise produced by the amplifiers on the spectrum estimation accuracy, as depicted in Figure 1.

Figure 1.

Noise measurement platform.

From Figure 1, we can see that the dual-channel output diode signals, which are passing through different amplifiers, can be expressed by α-stable noise and additive white Gaussian noise (WGN) through the same linear system, and both signals are corrupted by different additive CGN introduced, respectively, by the two amplifiers. The two different CGN of amplifiers are independent, because the amplifiers are powered by independent electric pools in this cross-spectrum measurement technique. This measurement method can improve the accuracy of spectrum estimation by effectively eliminating the influence of zero-mean CGN and/or WGN on it.

Applying the cross-spectrum measurement technique, by passing through dual amplifiers, yields dual-channel observed data. We address the problem of estimating ARMA coefficients of processes $Y\left(n\right)$ and $Z\left(n\right)$

$$\begin{array}{l}Y\left(n\right)=X\left(n\right)+N_1\left(n\right)=-{\displaystyle \sum _{k=1}^{p}a\left(k\right)X\left(n-k\right)}\\ +{\displaystyle \sum _{j=0}^{q}b\left(j\right)}U\left(n-j\right)+N_1\left(n\right)\end{array}$$

(3)

$$\begin{array}{l}Z\left(n\right)=X\left(n\right)+N_2\left(n\right)=-{\displaystyle \sum _{k=1}^{p}a\left(k\right)X\left(n-k\right)}\\ +{\displaystyle \sum _{j=0}^{q}b\left(j\right)}U\left(n-j\right)+N_2\left(n\right)\end{array}$$

(4)

ARMA model assumptions are as follows:

[AS1] The order ( $p,q$ ) of ARMA model is determined.

[AS2] The system is causal, linear time invariant, exponentially stable and minimum phase, i.e. the roots of the $A\left(z\right)$ polynomial lie within the unit circle; in the complex Z domain, the transfer function $H\left(z\right)$ of the model is conforming to pole-zero cancellations, which mean all poles lie inside of the unit circle. $A\left(z\right)$ and $B\left(z\right)$ have no common factors.

[AS3] $a\left(0\right)=b\left(0\right)=1$ , $a\left(p\right)\ne 0$ , $b\left(q\right)\ne 0$ .

[AS4] The input sequence $\left\{U\left(n\right)\right\}$ is considered to be unobservable, zero-mean, i.i.d. non-Gaussian, and α-stable process with characteristic exponent $\alpha \left(1\prec \alpha \prec 2\right)$ .

[AS5] The 1/f signal $\left\{X\left(n\right)\right\}$ is considered to be unobservable, zero-mean, non-Gaussian and α-stable process, with characteristic exponent $\alpha \left(1\prec \alpha \prec 2\right)$ ,which is the same as $\left\{U\left(n\right)\right\}$ .

[AS6] The additive noise $N_1\left(n\right)$ and $N_2\left(n\right)$ are unobservable, zero-mean colored Gaussian processes, respectively, with variance ${\sigma }_1^2$ and ${\sigma }_2^2$ , and they are, respectively, independent of the input $\left\{U\left(n\right)\right\}$ .

[AS7] The noisy outputs $Y\left(n\right)$ and $Z\left(n\right)$ are considered to be observable, zero-mean, non-Gaussian, and α-stable processes, with characteristic exponent ${\alpha }_Y\left(1\prec {\alpha }_Y\prec 2\right)$ and ${\alpha }_Z\left(1\prec {\alpha }_Z\prec 2\right)$ , which is different from $\left\{X\left(n\right)\right\}$ and $\left\{U\left(n\right)\right\}$ .

SNCC and SNCC-based spectrum estimation in SαS process

Definitions of SNCC

From the above, we can conclude that the accuracy of spectrum estimation can be improved by a method utilizing cross-spectrum measurement technique, while the method can achieve such goal by means of SNCC function.

The definition of SNCC is similar to that of samples normalized self-correlations (SNSC). According to literature, ^{10 ,11} SNSC $R_{\mathit{xx}}$ of $X\left(n\right)$ converges in probability to a certain constants as $n\to \infty$ :

$$\begin{array}{c}{R}_{\mathit{xx}}=\frac{\frac{1}{L}{\displaystyle {\sum }_{n=1}^L\left[X\left(n\right)X\left(n+m\right)\right]}}{\frac{1}{L}{\displaystyle {\sum }_{n=1}^L\left[X^2\left(n\right)\right]}}=\frac{{\displaystyle {\sum }_{n=1}^L\left[X\left(n\right)X\left(n+m\right)\right]}}{{\displaystyle {\sum }_{n=1}^L\left[X^2\left(n\right)\right]}}\end{array}$$

(5)

where $L$ is samples number and $\left|m\right|$ is far smaller than $L$ .

The following identities can be used to establish the existence and boundedness of SNSC

$${\displaystyle \sum _{n=1}^L{\left[X\left(n\right)\pm X\left(n+m\right)\right]}^2}\ge 0$$

(6)

By expanding equation (6)

$$\begin{array}{l}{\displaystyle \sum _{n=1}^L{\left[X\left(n\right)\pm X\left(n+m\right)\right]}^2}\\ ={\displaystyle \sum _{n=1}^L\left[X^2\left(n\right)\right]}+{\displaystyle \sum _{n=1}^L\left[X^2\left(n+m\right)\right]}\pm 2{\displaystyle \sum _{n=1}^L\left[X\left(n\right)X\left(n+m\right)\right]}\end{array}$$

(7)

In virtue of $X\left(n\right)$ satisfying α-stable distribution, therefore it has a stable ergodicity, ¹² when $L\to \infty$ , it has

$${\displaystyle \sum _{n=1}^L\left[X^2\left(n\right)\right]}={\displaystyle \sum _{n=1}^L\left[X^2\left(n+m\right)\right]}$$

(8)

In order to ensure the existence of $E\left[X\left(n\right)\right]$ , the bounds for $\alpha$ is: $1\prec \alpha \prec 2$ . Substituting equation (8) in equation (6)

$$\begin{array}{l}{\displaystyle \sum _{n=1}^L{\left[X\left(n\right)\pm X\left(n+m\right)\right]}^2}\\ =2{\displaystyle \sum _{n=1}^L\left[X^2\left(n\right)\right]}\pm 2{\displaystyle \sum _{n=1}^L\left[X\left(n\right)X\left(n+m\right)\right]}\ge 0\end{array}$$

(9)

and ${\sum }_{n=1}^L\left[X^2\left(n\right)\right]\succ 0$ , then $\pm {\sum }_{n=1}^L\left[X\left(n\right)X\left(n+m\right)\right]\ge -{\sum }_{n=1}^L\left[X^2\left(n\right)\right]$ . So

$$\left|\frac{{\displaystyle {\sum }_{n=1}^L\left[X\left(n\right)X\left(n+m\right)\right]}}{{\displaystyle {\sum }_{n=1}^L\left[X^2\left(n\right)\right]}}\right|\le 1$$

(10)

when $L\to \infty$ , it has

$$\begin{array}{l}\underset{N\to \infty }{\lim }{R}_{\mathit{XX}}=\frac{\underset{L\to \infty }{\lim }\frac{1}{L}{\displaystyle {\sum }_{n=1}^L\left[X\left(n\right)X\left(n+m\right)\right]}}{\underset{L\to \infty }{\lim }\frac{1}{L}{\displaystyle {\sum }_{n=1}^L\left[X^2\left(n\right)\right]}}\\ =\frac{E\left[X\left(n\right)X\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]}\end{array}$$

(11)

Then

$$\left|\underset{L\to \infty }{\lim }{R}_{\mathit{XX}}\right|=\left|\frac{E\left[X\left(n\right)X\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]}\right|\le 1$$

(12)

since $X\left(n\right)$ , $N_1\left(n\right)$ , and $N_2\left(n\right)$ are independent, $N_1\left(n\right)$ and $N_2\left(n\right)$ are zero-mean CGN, according to equations (3) and (4). Let

$$\begin{array}{l}{g}_{\mathit{XX}}=\frac{E\left[X\left(n\right)X\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]},\\ g_{\mathit{YZ}}\left(m\right)=\frac{E\left[Y\left(n\right)Z\left(n+m\right)\right]}{E\left[Y\left(n\right)Z\left(n\right)\right]},\\ g_{\mathit{ZY}}\left(m\right)=\frac{E\left[Z\left(n\right)Y\left(n+m\right)\right]}{E\left[Z\left(n\right)Y\left(n\right)\right]},\end{array}$$

Then

$$\begin{array}{l}{g}_{\mathit{YZ}}\left(m\right)\\ =\frac{E\left[Y\left(n\right)Z\left(n+m\right)\right]}{E\left[Y\left(n\right)Z\left(n\right)\right]}\\ =\frac{E\left\{\left[X\left(n\right)+N_1\left(n\right)\right]\left[X\left(n+m\right)+N_2\left(n+m\right)\right]\right\}}{E\left\{\left[X\left(n\right)+N_1\left(n\right)\right]\left[X\left(n\right)+N_2\left(n\right)\right]\right\}}\\ =\frac{\left\{\begin{array}{c}E\left[X\left(n\right)X\left(n+m\right)\right]+E\left[X\left(n\right)N_2\left(n+m\right)\right]\\ +E\left[X\left(n+m\right)N_1\left(n\right)\right]+E\left[N_1\left(n\right)N_2\left(n+m\right)\right]\end{array}\right\}}{\left\{\begin{array}{c}E\left[X^2\left(n\right)\right]+E\left[X\left(n\right)N_2\left(n\right)\right]\\ +E\left[X\left(n\right)N_1\left(n\right)\right]+E\left[N_1\left(n\right)N_2\left(n\right)\right]\end{array}\right\}}\\ {\to }_p\frac{E\left[X\left(n\right)X\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]}\end{array}$$

(13)

Then

$$g_{\mathit{YZ}}\left(m\right)=\frac{E\left[X\left(n\right)X\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]}=g_{\mathit{XX}}\left(m\right)$$

(14)

So $g_{\mathit{YZ}}\left(m\right)$ of observed data can be used as a substitute for $g_{\mathit{XX}}\left(m\right)$ of noiseless data. Similarly, $g_{\mathit{ZY}}\left(m\right)$ of observed data can also be used as a substitute for $g_{\mathit{XX}}\left(m\right)$ .

Properties of SNSC

Following are some significant properties of SNSC:

If $X$ , $X_1$ , and $X_2$ are jointly $S\alpha S$ , then SNSC of $g_{\mathit{XX}}\left(m\right)$ is linear in $X$

$$\begin{array}{l}{g}_{\mathit{XX}}\left(m\right)=\frac{E\left\{X\left(n\right)\left[A{X}_1\left(n+m\right)+B{X}_2\left(n+m\right)\right]\right\}}{E\left[X^2\left(n\right)\right]}\\ =A\frac{E\left[X\left(n\right)X_1\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]}+B\frac{E\left[X\left(n\right)X_2\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]}\end{array}$$

for arbitrary real constants $A$ and $B$ .

2. If $X$ and $U$ are independent and jointly $S\alpha S$ , according to Kanter and Steiger, ¹¹ then

$$g_{\mathit{XU}}\left(m\right)=\frac{E\left[X\left(n\right)U\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]}=0$$

3. Let $X\left(n+m\right)={\displaystyle {\sum }_{k=1}^{j}a\left(k\right)}X\left(n+m-k\right)$ , then

$$\begin{array}{l}{g}_{\mathit{XX}}\left(m\right)=\frac{E\left[X\left(n\right)X\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]}\\ =\frac{E\left\{X\left(n\right)\left[{\displaystyle {\sum }_{k=1}^{j}a\left(k\right)X\left(n+m-k\right)}\right]\right\}}{E\left[X^2\left(n\right)\right]}\\ ={\displaystyle \sum _{k=1}^{j}a\left(k\right)}\frac{E\left[X\left(n\right)X\left(n+m-k\right)\right]}{E\left[X^2\left(n\right)\right]}\\ ={\displaystyle \sum _{k=1}^{j}a\left(k\right)}{g}_{\mathit{XX}}\left(m-k\right)\end{array}$$

SNCC-based Yule–Walker equation of ARMA model

The SNCC functions of processes $Y\left(n\right)$ and $Z\left(n\right)$ is

$$\begin{array}{l}{g}_{\mathit{YZ}}\left(m\right)=g_{\mathit{ZY}}\left(m\right)=g_{\mathit{XX}}\left(m\right)=\frac{E\left[X\left(n\right)X\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]}\\ =\frac{E\left\{X\left(n\right)\left[{\displaystyle {\sum }_{k=1}^p-\left[a\left(k\right)X\left(n+m-k\right)\right]}+{\displaystyle {\sum }_{j=0}^{q}b\left(j\right)}U\left(n+m-j\right)\right]\right\}}{E\left[X^2\left(n\right)\right]}\\ =\frac{-a\left(1\right)E\left[X\left(n+m-1\right)X\left(n\right)\right]-\cdots -a\left(m\right)E\left[X^2\left(n\right)\right]}{E\left[X^2\left(n\right)\right]}\\ +\frac{\left\{\begin{array}{l}-a\left(m+1\right)E\left[X\left(n-1\right)X\left(n\right)\right]-\cdots -a\left(p\right)E\left[X\left(n+m-p\right)X\left(n\right)\right]\\ +E\left\{X\left(n\right)\left[{\displaystyle {\sum }_{j=0}^{q}b\left(j\right)}U\left(n+m-j\right)\right]\right\}\end{array}\right\}}{E\left[X^2\left(n\right)\right]}\\ =-a\left(1\right)g_{\mathit{XX}}\left(m-1\right)-\cdots -a\left(m\right)-a\left(m+1\right)\left[g_{\mathit{XX}}\left(1\right)-\frac{X\left(N\right)X\left(N+1\right)}{E\left[X^2\left(n\right)\right]}\right]-\cdots \\ -a\left(p\right)\left[g_{\mathit{XX}}\left(p-m\right)-\frac{\left\{\begin{array}{l}X\left(N\right)X\left(N+p-m\right)-\cdots \\ -X\left(N-p+m+1\right)X\left(N+1\right)\end{array}\right\}}{E\left[X^2\left(n\right)\right]}\right]+\frac{E\left\{X\left(n\right)\left[{\displaystyle {\sum }_{j=0}^{q}b\left(j\right)}U\left(n+m-j\right)\right]\right\}}{E\left[X^2\left(n\right)\right]}\end{array}$$

(15)

$$\begin{array}{l}\frac{E\left[X\left(n\right)\left[{\displaystyle \sum _{j=0}^{q}b\left(j\right)U\left(n+m-j\right)}\right]\right]}{E\left[X^2\left(n\right)\right]}\\ =\frac{b\left(0\right)E\left[X\left(n\right)U\left(n+m\right)\right]+b\left(1\right)E\left[X\left(n\right)U\left(n+m-1\right)\right]+\cdots +b\left(m\right)E\left[X\left(n\right)U\left(n\right)\right]}{E\left[X^2\left(n\right)\right]}\\ +\cdots \frac{b\left(m+1\right)E\left[X\left(n\right)U\left(n-1\right)\right]+\cdots +b\left(q\right)E\left[X\left(n\right)U\left(n+m-q\right)\right]}{E\left[X^2\left(n\right)\right]}\\ =b\left(0\right)g_{\mathit{XU}}\left(m\right)+b\left(1\right)g_{\mathit{XU}}\left(m-1\right)+\cdots +b\left(m\right)g_{\mathit{XU}}\left(0\right)+b\left(m+1\right)\left[g_{\mathit{XU}}\left(1\right)-\frac{X\left(N\right)U\left(N+1\right)}{E\left[X^2\left(n\right)\right]}\right]\cdots +\\ +b\left(q\right)\left[g_{\mathit{XU}}\left(q-m\right)-\frac{X\left(N\right)U\left(N+q-m\right)-\cdots -X\left(N-q+m+1\right)U\left(N+1\right)}{E\left[X^2\left(n\right)\right]}\right]\end{array}$$

(16)

According to Shi et al., ¹²

$$\begin{array}{l}\frac{X(N)X(N+k)+\cdots +X(N-k+1)X(N+1)}{E[X^2(n)]}{\to }_{p}0\\ \left(\text{convergence}\text{in}\text{probability}\right).\end{array}$$

For a fixed integer $m$ ,

$$\begin{array}{l}{g}_{\mathit{XU}}\left(m\right)=\frac{E\left[X\left(n\right)U\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]}{\to }_{p}0\\ \left(\text{convergence}\text{in}\text{probability}\right).\end{array}$$

For

$$\begin{array}{l}1\le j\le q-m,\frac{\left\{\begin{array}{l}X(N)U(N+j)+\cdots \\ +X(N-j+1)U(N+1)\end{array}\right\}}{E[X^2(n)]}{\to }_{p}0,\\ \left(\text{convergence}\text{in}\text{probability}\right).\end{array}$$

So

$$\begin{array}{c}\frac{E\left[X(n)\left[{\displaystyle {\sum }_{j=0}^{q}b(j)U(n+m-j)}\right]\right]}{E[X^2(n)]}{\to }_{p}0,b(j)\ne 0,\\ \left(\text{convergence}\text{in}\text{probability}\right).\end{array}$$

Then

$$\begin{array}{l}-g_{\mathit{XX}}\left(m\right)=a\left(1\right)g_{\mathit{XX}}\left(m-1\right)\\ +\cdots +a\left(m\right)+a\left(m+1\right)g_{\mathit{XX}}\left(1\right)\\ +\cdots a\left(p\right)g_{\mathit{XX}}\left(p-m\right),1\le m\le p\prec q\end{array}$$

(17)

Consistency of parameter estimation for ARMA model

Theorem 1: The natural estimate ${\widehat{g}}_{\mathit{XX}}\left(m\right)$ of $g_{\mathit{XX}}\left(m\right)=\frac{E\left[X\left(n\right)X\left(n+m\right)\right]}{E\left[X^2\left(n\right)\right]}$ is consistent.

Proof: Since ${\widehat{g}}_{\mathit{XX}}\left(m\right)=\frac{{\sum }_{n=1}^L\left[X\left(n\right)X\left(n+m\right)\right]}{{\sum }_{n=1}^L\left[X^2\left(n\right)\right]}$ is absolutely summable, the order of $\lim$ and $\Sigma \hspace{0.17em}$ can be exchanged

$$\begin{array}{l}E\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]=\frac{1}{L}{\displaystyle \sum _{n=1}^L\left\{\frac{{\displaystyle {\sum }_{n=1}^L\left[X\left(n\right)X\left(n+m\right)\right]}}{{\displaystyle {\sum }_{n=1}^L\left[X^2\left(n\right)\right]}}\right\}}\\ \underset{L\to \infty }{\to }\frac{1}{L}{\displaystyle \sum _{n=1}^L\left[g_{\mathit{XX}}\left(m\right)\right]}=g_{\mathit{XX}}\left(m\right)\end{array}$$

(18)

So ${\widehat{g}}_{\mathit{XX}}\left(m\right)$ is an unbiased estimate of $g_{\mathit{XX}}\left(m\right)$ .

$$\begin{array}{l}V\mathit{ar}\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]\\ =E\left\{{\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)-E\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]\right]}^2\right\}\\ =E\left\{{\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]}^2\right\}-{\left\{E\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]\right\}}^2\\ =\underset{L\to \infty }{\lim }\frac{1}{L}{\displaystyle \sum _{n=1}^L\left\{{\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]}^2\right\}}-{\left\{E\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]\right\}}^2\end{array}$$

(19)

According to equation (12), $\left|{\widehat{g}}_{\mathit{XX}}\left(m\right)\right|\le 1$ , so ${\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]}^2\le 1$ . Since ${\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]}^2$ is absolutely summable, the order of $\lim$ and $\Sigma \hspace{0.17em}$ can be exchanged, then

$$\begin{array}{l}V\mathit{ar}\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]\\ =\frac{1}{L}{\displaystyle \sum _{n=1}^L\underset{L\to \infty }{\lim }}\left\{{\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]}^2\right\}-{\left\{E\left[{\widehat{g}}_{\mathit{XX}}\left(m\right)\right]\right\}}^2\end{array}$$

$$\begin{array}{l}=\frac{1}{L}{\displaystyle \sum _{n=1}^L\left\{{\left[g_{\mathit{XX}}\left(m\right)\right]}^2\right\}}-{\left[g_{\mathit{XX}}\left(m\right)\right]}^2\\ ={\left[g_{\mathit{XX}}\left(m\right)\right]}^2-{\left[g_{\mathit{XX}}\left(m\right)\right]}^2=0\end{array}$$

(20)

It suffices to establish the consistency of parameter estimation for ARMA model.

A SNCC-based ARMA parameter estimation algorithm

Estimation of AR coefficients

Given the SNCC-based equation of ARMA model, we apply equation (17) to write into the form of Yule–Walker equation

$$\begin{array}{l}\left[\begin{array}{ccccc}1& g_{\mathit{XX}}\left(1\right)& g_{\mathit{XX}}\left(2\right)& \cdots & g_{\mathit{XX}}\left(p-1\right)\\ g_{\mathit{XX}}\left(1\right)& 1& g_{\mathit{XX}}\left(1\right)& \cdots & g_{\mathit{XX}}\left(p-2\right)\\ g_{\mathit{XX}}\left(2\right)& g_{\mathit{XX}}\left(1\right)& 1& \cdots & g_{\mathit{XX}}\left(p-3\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ g_{\mathit{XX}}\left(p-1\right)& g_{\mathit{XX}}\left(p-2\right)& g_{\mathit{XX}}\left(p-3\right)& \cdots & 1\end{array}\right]\\ \times \left[\begin{array}{c}a\left(1\right)\\ a\left(2\right)\\ a\left(3\right)\\ ⋮\\ a\left(p\right)\end{array}\right]=-\left[\begin{array}{c}{g}_{\mathit{XX}}\left(1\right)\\ g_{\mathit{XX}}\left(2\right)\\ g_{\mathit{XX}}\left(3\right)\\ ⋮\\ g_{\mathit{XX}}\left(p\right)\end{array}\right]\end{array}$$

(21)

Let

$$G=\left[\begin{array}{ccccc}1& g_{\mathit{XX}}\left(1\right)& g_{\mathit{XX}}\left(2\right)& \cdots & g_{\mathit{XX}}\left(p-1\right)\\ g_{\mathit{XX}}\left(1\right)& 1& g_{\mathit{XX}}\left(1\right)& \cdots & g_{\mathit{XX}}\left(p-2\right)\\ g_{\mathit{XX}}\left(2\right)& g_{\mathit{XX}}\left(1\right)& 1& \cdots & g_{\mathit{XX}}\left(p-3\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ g_{\mathit{XX}}\left(p-1\right)& g_{\mathit{XX}}\left(p-2\right)& g_{\mathit{XX}}\left(p-3\right)& \cdots & 1\end{array}\right],$$

$$\begin{array}{l}A={\left[\begin{array}{ccccc}a\left(1\right)& a\left(2\right)& a\left(3\right)& \cdots & a\left(p\right)\end{array}\right]}^T,\\ g=-{\left[\begin{array}{ccccc}{g}_{\mathit{XX}}\left(1\right)& g_{\mathit{XX}}\left(2\right)& g_{\mathit{XX}}\left(3\right)& \cdots & g_{\mathit{XX}}\left(p\right)\end{array}\right]}^T,\end{array}$$

then $\mathit{GA}=g$ .

The coefficients $a\left(1\right),\dots ,a\left(p\right)$ are unique if and only if $g_{\mathit{XX}}\left(1\right),\dots ,g_{\mathit{XX}}\left(p\right)$ are linearly independent elements in the space of integrable random variables.

Theorem 2: $G$ is a Full rank matrix.

Proof: Let

$$G^M=\left[\begin{array}{ccccc}1& g_{\mathit{XX}}\left(1\right)& g_{\mathit{XX}}\left(2\right)& \cdots & g_{\mathit{XX}}\left(p-1\right)\\ g_{\mathit{XX}}\left(1\right)& 1& g_{\mathit{XX}}\left(1\right)& \cdots & g_{\mathit{XX}}\left(p-2\right)\\ g_{\mathit{XX}}\left(2\right)& g_{\mathit{XX}}\left(1\right)& 1& \cdots & g_{\mathit{XX}}\left(p-3\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ g_{\mathit{XX}}\left(p-1\right)& g_{\mathit{XX}}\left(p-2\right)& g_{\mathit{XX}}\left(p-3\right)& \cdots & 1\\ g_{\mathit{XX}}\left(p\right)& g_{\mathit{XX}}\left(p-1\right)& g_{\mathit{XX}}\left(p-2\right)& \cdots & g_{\mathit{XX}}\left(1\right)\\ g_{\mathit{XX}}\left(p+1\right)& g_{\mathit{XX}}\left(p\right)& g_{\mathit{XX}}\left(p-1\right)& \cdots & g_{\mathit{XX}}\left(2\right)\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ g_{\mathit{XX}}\left(M\right)& g_{\mathit{XX}}\left(M-1\right)& g_{\mathit{XX}}\left(M-2\right)& \cdots & g_{\mathit{XX}}\left(M-p+1\right)\end{array}\right]$$

(22)

where $M$ is an arbitrary positive integer greater than $p$ , let $g_m{\hspace{0.17em}}^M$ be the $m$ th line row vector of the matrix $G^M$ , it is known by equation (22)

$$g_{p+1}{\hspace{0.17em}}^M=\left[\begin{array}{cccc}{g}_p{\hspace{0.17em}}^M\cdot A& g_{p-1}{\hspace{0.17em}}^M\cdot A& \cdots & g_1{\hspace{0.17em}}^M\cdot A\end{array}\right]$$

where $g_{p+1}{\hspace{0.17em}}^M$ is the linear combination of vectors from pth to the first row vector of the matrix $G^M$ .

$$g_{p+2}{\hspace{0.17em}}^M=\left[\begin{array}{cccc}{g}_{p+1}{\hspace{0.17em}}^M\cdot A& g_p{\hspace{0.17em}}^M\cdot A& \cdots & g_2{\hspace{0.17em}}^M\cdot A\end{array}\right]$$

where $g_{p+2}{\hspace{0.17em}}^M$ is the linear combination of vectors from ( $p+1$ )th to the second row vector of the matrix $G^M$ .

It can be deduced that $g_M{\hspace{0.17em}}^M=\left[\begin{array}{cccc}{g}_{M-1}{\hspace{0.17em}}^M\cdot A& g_{M-2}{\hspace{0.17em}}^M\cdot A& \cdots & g_{M-p+1}{\hspace{0.17em}}^M\cdot A\end{array}\right]$ .

Then, the non-singularity of the matrix $G^M$ depends only on the former $p$ row vectors.

Through the reverse proving, suppose the former $p$ row vectors $g_1{\hspace{0.17em}}^M,g_2{\hspace{0.17em}}^M,\dots ,g_p{\hspace{0.17em}}^M$ of matrix $G^M$ have a significant linear correlation, since $g_{p+1}{\hspace{0.17em}}^M$ is linearly dependent of $g_1{\hspace{0.17em}}^M,g_2{\hspace{0.17em}}^M,\dots ,g_p{\hspace{0.17em}}^M$ , $g_{p+2}{\hspace{0.17em}}^M$ is linearly dependent of $g_1{\hspace{0.17em}}^M,g_2{\hspace{0.17em}}^M,\dots ,g_p{\hspace{0.17em}}^M$ , in a similar fashion, matrix $G^M$ is only linear independent of the front ( $p-1$ ) row vectors, so there exists such a $\left(p-1\right)\times 1$ dimension parameter vector:

$$\begin{array}{l}{g}_{\mathit{XX}}\left(m\right)=-c\left(1\right)g_{\mathit{XX}}\left(m-1\right)-\cdots -c\left(m\right)\\ -c\left(m+1\right)g_{\mathit{XX}}\left(1\right)-\cdots \\ -c\left(p-1\right)g_{\mathit{XX}}\left(p-m-1\right),1\le m\le p.\end{array}$$

It obviously leads to the conclusion that $a\left(p\right)=0$ , which contradicts the AR model assumption [AS3]. So it is well explained that the former $p$ row vectors $g_1{\hspace{0.17em}}^M,g_2{\hspace{0.17em}}^M,\dots ,g_p{\hspace{0.17em}}^M$ of the matrix $G^M$ must be linearly independent. As a result of matrix $G$ being symmetric, the former $p$ row vectors are linearly independent of each other. In summary, matrix $G$ is full rank with its rank $p$ .

Estimation of MA coefficients

Taking advantage of the AR model parameter $A_1={\left[\begin{array}{cccc}1& {\tilde{a}}_1& \cdots & {\tilde{a}}_p\end{array}\right]}^T$ to constitute $\tilde{A}\left(z\right)$ , $\tilde{A}\left(z\right)=1+{\sum }_{k=1}^p\tilde{a}\left(k\right)z^{-k}$ , cascading $\tilde{A}\left(z\right)$ with the original ARMA model, then: $\tilde{X}\left(n\right)={\sum }_{k=0}^{p}a\left(k\right)X\left(n-k\right)={\sum }_{j=0}^{q}b\left(j\right)U\left(n-j\right)$ , $\tilde{X}\left(n\right)$ is conversed to pure MA system

$${\displaystyle \sum _{j=0}^{q}b\left(j\right)}a\left(n-j\right)=\delta \left(n\right),n=0,\dots ,p,b\left(0\right)=1$$

(23)

Then

$$\tilde{a}\left(k\right)+{\displaystyle \sum _{j=1}^{q}b\left(j\right)}\tilde{a}\left(k-j\right)=e\left(k\right)$$

(24)

where $e\left(k\right)$ represents error sequence, $\tilde{a}\left(0\right)=0$ , with $\left\{\tilde{a}\left(k\right)\right\},1\le k\le p$ as observation sequence, and $\left\{b\left(j\right)\right\},1\le j\le q$ as coefficients of AR( $q$ ) model.

Since $\left\{\tilde{a}\left(k\right)\right\}$ is a finite method sequence, estimation of $\left\{b\left(j\right)\right\}$ is obtained by means of Yule–Walker equation.

The Yule–Walker equation of AR( $q$ ) model is as follows

$$\left[\begin{array}{cccc}r\left(0\right)& r\left(-1\right)& \cdots & r\left(1-q\right)\\ r\left(1\right)& r\left(0\right)& \cdots & r\left(2-q\right)\\ ⋮& ⋮& \ddots & ⋮\\ r\left(q-1\right)& r\left(q-1\right)& \cdots & r\left(0\right)\end{array}\right]\left[\begin{array}{c}b\left(1\right)\\ b\left(2\right)\\ ⋮\\ b\left(q\right)\end{array}\right]=-\left[\begin{array}{c}r\left(1\right)\\ r\left(2\right)\\ ⋮\\ r\left(q\right)\end{array}\right]$$

(25)

Using the estimated value $\tilde{r}\left(j\right)$ of the auto-correlation sequence $\left\{\tilde{a}\left(k\right)\right\}$ to take the place of $r\left(j\right)$ in the matrix

$$\tilde{r}\left(j\right)=\{\begin{array}{l}\frac{1}{p+1}{\displaystyle \sum _{k=0}^{p-m}\tilde{a}\left(k\right)\tilde{a}\left(k+j\right),}0\le j\le q\hfill \\ \tilde{r}\left(-j\right),1-q\le j\le -1\hfill \end{array}$$

(26)

the estimated value of the MA parameter sequence is obtained: $\left[\begin{array}{cccc}1& \tilde{b}\left(1\right)& \cdots & \tilde{b}\left(q\right)\end{array}\right]$ .

Singular value decomposition-Total Least Squares (SVD-TLS) algorithm for ARMA parameter estimation

For the AR parameter estimation part, we select an extended order $p_e$ to construct Yule–Walker equation of matrix

$$\begin{array}{l}\left[\begin{array}{ccccc}{g}_{\mathit{XX}}\left(1\right)& 1& g_{\mathit{XX}}\left(1\right)& \cdots & g_{\mathit{XX}}\left(p_e-1\right)\\ g_{\mathit{XX}}\left(2\right)& g_{\mathit{XX}}\left(1\right)& 1& \cdots & g_{\mathit{XX}}\left(p_e-2\right)\\ g_{\mathit{XX}}\left(3\right)& g_{\mathit{XX}}\left(2\right)& g_{\mathit{XX}}\left(1\right)& \cdots & g_{\mathit{XX}}\left(p_e-3\right)\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ g_{\mathit{XX}}\left(l\right)& g_{\mathit{XX}}\left(l-1\right)& g_{\mathit{XX}}\left(l-2\right)& \cdots & g_{\mathit{XX}}\left(l-p_e+1\right)\end{array}\right]\\ \times \left[\begin{array}{c}1\\ a\left(1\right)\\ a\left(2\right)\\ ⋮\\ a\left(p_e\right)\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ ⋮\\ 0\end{array}\right]\end{array}$$

(27)

where $p_e\ge p$ , $l\succ p+1$ .

Writing equation (27) in brief

$$G_e{A}_e=0$$

(28)

SVD decomposition of matrix $G_e$ can be obtained

$$G_e=U\Lambda V^{*}$$

(29)

$U$ and $V$ are, respectively, the left and right singular vectors of $G_{\mathit{ge}}$ , that is

$$G_e{G}_e^{*}{u}_i={\sigma }_i{u}_i,G_e{G}_e^{*}{v}_i={\sigma }_i{v}_i$$

$\Lambda =\mathit{diag}\left[{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_{p_e+1}\right]$ ranges according to the descending order of ${\sigma }_i$ , that is ${\sigma }_1\succ {\sigma }_2\succ \cdots \succ {\sigma }_{p_{\text{e}}+1}$ .

Computing $G_e{\hspace{0.17em}}^{\left(p\right)}={\sum }_{i=1}^p{\sum }_{k=1}^{p_e-p+1}{\sigma }_i{u}_i{\hspace{0.17em}}^k{v}_i{{\hspace{0.17em}}^k}^{*}$ , where $u_{\mathit{i}}^k$ represents an $\left(p+1\right)\times 1$ dimensional vector consisting of $\left(p+1\right)$ elements of the $i$ column vector in a matrix $U$ , and

$$\begin{array}{l}{u}_{\mathit{i}}^k=\left[u_i\left(k\right),u_i\left(k+1\right),\dots ,u_i\left(k+p\right)\right],\\ 1\le k\le p_e-p+1,1\le i\le p\end{array}$$

Similarly,

$$\begin{array}{l}{v}_{\mathit{i}}^k=\left[v_i\left(k\right),v_i\left(k+1\right),\dots ,v_i\left(k+p\right)\right],\\ 1\le k\le p_e-p+1,1\le i\le p\end{array}$$

Solving linear equations $G_e^{\left(p\right)}{A}_1=a{e}_1$ , where $A_1={\left[\begin{array}{cccc}1& {\tilde{a}}_1& \cdots & {\tilde{a}}_p\end{array}\right]}^T$ , $e_1={\left[\begin{array}{cccc}1& 0& \cdots & 0\end{array}\right]}^T$ , $A_1$ is the estimated parameter vector of the AR model.

The MA parameter estimation is analogous to that of AR. Selecting an extended order $q_e$ to construct Yule–Walker equation of matrix with extended order

$$\left[\begin{array}{cccc}r\left(0\right)& r\left(-1\right)& \cdots & r\left(1-q_e\right)\\ r\left(1\right)& r\left(0\right)& \cdots & r\left(2-q_e\right)\\ ⋮& ⋮& \ddots & ⋮\\ r\left(s\right)& r\left(s-1\right)& \cdots & r\left(s-q_e+1\right)\end{array}\right]\left[\begin{array}{c}b\left(1\right)\\ b\left(2\right)\\ ⋮\\ b\left(q_e\right)\end{array}\right]=-\left[\begin{array}{c}r\left(1\right)\\ r\left(2\right)\\ ⋮\\ r\left(q_e\right)\end{array}\right]$$

(30)

where $q_e\ge q$ , $s\succ q+1$ .

Then, MA parameter vector $B_1=\left[\begin{array}{cccc}1& {\tilde{b}}_1& \cdots & {\tilde{b}}_q\end{array}\right]$ can be estimated.

Experiments

Simulation and experiments

Analogous to dual-channel output data, the composition of the data includes the signal (1/f signal conforming to alpha-stable distribution) and the additive noise (CGN). The 1/f signal is generated according to Jeremy Usdin. ¹³ We generate a 1/f process, adding two different colored Gaussian processes respectively.

According to “Linear SαS process in the observed data” section, the characteristic parameter of $U\left(n\right)$ is set to $\alpha \left(1\prec \alpha \prec 2\right)$ , the dispersion coefficient $\gamma$ is set to 1, symmetric parameter $\beta$ and location parameter $\delta$ are both set to 0. The characteristic parameter of generated signal $X\left(n\right)$ is equal to that of $U\left(n\right)$ .

Mixed signal-to-noise ratio

With non-existence of finite variance, the variance of conventional signal-to-noise ratio (SNR) stands no point; based on that, we use mixed SNR ^{6 ,14,15} to define SNR under α-stable distribution

$$R_{\mathit{MSNR}}=10\log 10\left(\gamma /{\sigma }_{\mathit{WN}}^2\right)=m$$

(31)

where $\gamma$ is the dispersion coefficient of signal $X\left(n\right)$ , ${\sigma }_{\mathit{WN}}^2$ is the variance of WGN, $R_{\mathit{MSNR}}$ is set to $m$ dB.

Error analysis of AR coefficients and MA coefficients

The method of Wang and Zhu ¹ estimating ARMA parameters using FLOC-based method under noise-free environment is chosen to compete with our methods of SNCC and SNSC estimating ARMA parameters under noisy environment.

From Tables 1 to 4, the performance of our method of estimating ARMA parameters based on SNCC is evaluated and compared with that of SNSC-based and FLOC-based ¹ methods. It is observed that our SNCC-based method exhibits good performance.

Table 1.

Mean (std) of AR coefficients over 100 trails, with 10,000 samples, alpha = 1.8, $R_{\mathit{MSNR}}$  = 20.

	a(1)	a(2)
True value	–0.5000	–0.1250
SNCC	–0.4982 (0.0214)	–0.1271 (0.0128)
SNSC	–0.4886 (0.0245)	–0.1333 (0.0126)
FLOC	–0.4573 (0.0467)	–0.1453 (0.0625)

SNCC: samples normalized cross-correlations; SNSC: samples normalized self-correlations; FLOC: fractional lower order covariance.

Table 2.

Mean (std) of MA coefficients over 100 trails, with 10,000 samples, alpha = 1.8, $R_{\mathit{MSNR}}$  = 20.

	b(1)	b(2)	b(3)
True value	0.10000	0.0550	0.0386
SNCC	0.0912 (0.0362)	0.0491 (0.0589)	0.0217 (0.0578)
SNSC	0.0763 (0.0418)	0.0352 (0.0618)	0.0186 (0.0715)
FLOC	0.0532 (0.3566)	0.1244 (0.4564)	0.0473 (0.5633)

SNCC: samples normalized cross-correlations; SNSC: samples normalized self-correlations; FLOC: fractional lower order covariance.

Table 3.

Mean (std) of AR coefficients over 100 trails, with 10,000 samples, alpha = 1.5, $R_{\mathit{MSNR}}$ =10.

	a(1)	a(2)	a(3)
True value	–0.5000	–0.1250	–0.0625
SNCC	–0.4972 (0.0144)	–0.1263 (0.0056)	–0.0623 (0.0073)
SNSC	–0.4854 (0.0186)	–0.1311 (0.0098)	–0.0569 (0.1109)
FLOC	–0.4678 (0.0356)	–0.1567 (0.1233)	–0.0489 (0.2345)

SNCC: samples normalized cross-correlations; SNSC: samples normalized self-correlations; FLOC: fractional lower order covariance.

Table 4.

Mean (std) of MA coefficients over 100 trails, with 10,000 samples, alpha = 1.5, $R_{\mathit{MSNR}}$  = 10.

	b(1)	b(2)	b(3)
True value	0.8645	0.3617	0.1081
SNCC	0.8845 (0.0867)	0.3122 (0.1567)	0.0865 (0.2355)
SNSC	0.646 (0.2356)	0.4367 (0.2465)	0.0663 (0.4577)
FLOC	1.1467 (0.7643)	0.8674 (0.7654)	1.3568 (1.8664)

SNCC: samples normalized cross-correlations; SNSC: samples normalized self-correlations; FLOC: fractional lower order covariance.

Conclusions

Our studies have showed that 1/f noises of transistors belong to α-stable distribution, and during the noise measurement, additive CGN will be introduced inevitably, so the problem of ARMA parameter estimation is transformed into α-stable noise ARMA parameter estimation problem in noisy environments. Due to the infinite variance of alpha-stable process, we proposed SNCC function to estimate ARMA parameters in CGN environments. The proposed SNCC function facilitates the estimation of the model parameters and also helps to suppress the effects of additive CGN. The experimental results have shown that the estimated spectrum by our method is well behaved.