##

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:

- 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);

- Compute the ACF features
z=A'*x;

- 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 .

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

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

- 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