Application of CLT to Auto-Correlation Function Analysis

We repeat the feature extraction of Section 5.2.2, but this time use the CLT. Once again, because the FFT is used, this approach falls under feature extraction methods for circularly stationary processes (See section 10.1).

The main problem in CLT is finding a positive-valued mean for the floating reference hypothesis, $\bar{{\bf x}}$, for which the feature mean is equal to the current feature value. Let ${\bf z}^*$ be the current ACF feature value. We need to find a positive-valued $\bar{{\bf x}}$ such that ${\bf A}^\prime \bar{{\bf x}}={\bf z}^*.$ A convenient approach comes from AR analysis, where it is known that the ACF of the AR spectrum is equal to the original ACF [31]. Thus, the AR spectral estimate serves as $\bar{{\bf x}}$.

Using the Levinson-Durbin recursion [31], we may transform the ACF ${\bf z}$ into the AR coefficients $\{a_1,a_2 \ldots a_P,\sigma^2\}$, and $a_0=1$. The corresponding AR spectrum is written

$$\bar{x}_k = {N_t \sigma^2 \over \left\vert \sum_{i=0}^P\; a_i \exp \left\{ -\frac{j2\pi i k}{N_t}\right\} \right\vert^2}, \;\;\; 1\leq k \leq N.$$

Under ${\bf H_0\mbox{\small $({\bf z})$}}$, the elements of ${\bf x}$ are independent and Chi-squared with 1 or 2 degrees of freedom with mean $\bar{k}_k$ (See Section 17.1.2). Bins $k=1,N$ are distributed according to

$$p_0(x_k,\bar{x}_k) = \frac{1}{\bar{x}_k \sqrt{2\pi}}\; (x_k/\bar{x}_k )^{-1/2} \; \exp\left\{ -\frac{x_k}{2\bar{x}_k }\right\},$$

while bins 2 through $N-1$ are exponentially distributed according to

$$p_1(x_k,y^r) = \frac{1}{\bar{x}_k}\; \exp\left\{-\frac{x_k}{\bar{x}_k}\right\}.$$

Putting it all together,

$$\log p({\bf x}\vert H_0$$$({\bf z})$$$) = \log p_0(x_1,\bar{x}_1) + {\displaystyle \sum_{k=2}^{N-1}} \log p_1(x_k,\bar{x}_k) + \log p_0(x_N,\bar{x}_{N}).$$

In summary, we take the following steps:

1. Create matrix ${\bf A}$ for ACF.

        % Nt = FFT size
        % x is dimension n by 1, where n=Nt/2+1;
        % A is dimension n by P+1
        % k is vector of degrees-of-freedom
        k = [1;2*ones(n-2,1); 1];
        A = cos([0:Nt/2]' *(2*pi/Nt)*[0:P]).*repmat(k,1,P+1)*(1/Nt^2);
   

2. Compute the ACF features

        z=A'*x;
   

3. Using the Levinson algorithm, compute the AR coefficients corresponding to the model-$P$ AR model, where $P+1=D$ the length of venctor ${\bf z}$.

        [a,e]=levinson(z);
   

   The output variable e is the AR innovation variance parameter $\sigma^2$.

4. Compute the AR spectral estimate.

     af = fft(a(:),Nt);
     xz = Nt * e  ./ abs(af(1:n)).^2;
   

5. Compute $p({\bf x}\vert H_0$$({\bf z})$$)$:

   lpxHz = expon_dist(x([2:n-1],:),xr([2:n-1],:)) + ...
       chisq1_dist(x([1 n],:),xr([1 n],:));
   

6. Verify the reproduction of the features

      format long
      [A'*xr z]
   

We then apply the results in Section 5.2.4. The variance of the elements of ${\bf x}$ under $H_0$$({\bf z})$ is $q_i$:

   q = 2./k .* xz.^2;

The covariance of ${\bf z}$ under $H_0$$({\bf z})$ is ${\bf\Sigma}_z$:

   Sz = A'*diag(q)*A;

Then, $\log p({\bf z}\vert H_0$$({\bf z})$$)$ is

    ldetSr=log(det(Sz));
    lpzHz =  -(P+1)/2*log(2*pi) - .5*ldetSz;

See software/module_acf_clt.m for additional details. Use software/module_acf_synth.m for inversion (re-synthesis). The module can be tested using software/module_acf_test.m with TYPE=1. We compare this method with other approaches in Section 10.4.10.