## 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 9.1).

The main problem in CLT is finding a positive-valued mean for the floating reference hypothesis, , for which the feature mean is equal to the current feature value. Let be the current ACF feature value. We need to find a positive-valued such that A convenient approach comes from AR analysis, where it is known that the ACF of the AR spectrum is equal to the original ACF [25]. Thus, the AR spectral estimate serves as .

Using the Levinson-Durbin recursion [25], we may transform the ACF into the AR coefficients , and . The corresponding AR spectrum is written

Under , the elements of are independent and Chi-squared with 1 or 2 degrees of freedom with mean (See Section 16.1.2). Bins are distributed according to

while bins 2 through are exponentially distributed according to

Putting it all together,

In summary, we take the following steps:

1. Create matrix 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- AR model, where the length of venctor .
     [a,e]=levinson(z);

The output variable e is the AR innovation variance parameter .

4. Compute the AR spectral estimate.
  af = fft(a(:),Nt);
xz = Nt * e  ./ abs(af(1:n)).^2;


5. Compute :
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 under is :
   q = 2./k .* xz.^2;

The covariance of under is :
   Sz = A'*diag(q)*A;

Then, 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 9.4.10.

Baggenstoss 2017-05-19