help-octave
[Top][All Lists]

## autocorrelation & family

 From: Francesco Potorti` Subject: autocorrelation & family Date: Wed, 03 May 1995 17:51 +0100 (MET)

```I wrote these simple functions for working with sequences of samples
of data.  They compute the positive autocovariance and autocorrelation
of the sequence, the sequence sample variance and the equivalent
number of independent observations, that is, the number by which the
sample variance is greater than the process variance.  Last, a
function for computing the one-sided power spectrum of a sequence.  I
tried windowing on the overlapping chunks used in this function, but
it strangely seems to give worse results than without windowing.

# Prevent Octave from thinking that this
# is a function file:
1;

function A = autocov(x)

# Given a vector X, computes a vector of the same shape, whose length
# is one less than the argument, representing the autocovariance of X.

A = x - mean(x);
A = fftconv(A, fliplr(flipud(A)));
A = A(length(x) : 2*length(x)-1);
A = A / length(x);
endfunction

function A = autocorr(x)

# Given a vector X, computes a vector of the same shape, whose length
# is one less than the argument, representing the autocorrelation of X.
# The autocorrelation is simply the normalised autocovariance.

A = autocov(x);
A = A / A(1);
endfunction

function t = samplestd(x)

# Given a vector X of samples, computes its sample variance.

weights = 1 - [1 : length(x)-1] / (length(x)-1);
x = reshape(x,1,length(x));
t = sqrt(2 * autocov(x) * [0.5, weights]);
endfunction

function t = corrtime(x)

# Given a vector X of samples, computes its "equivalent number of
# independent observations".

weights = 1 - [1 : length(x)-1] / (length(x)-1);
x = reshape(x,1,length(x));
t = autocorr(x) * [0.5, weights];
endfunction

function A = spectre(x, n)

# Given a vector X, compute its one-sided power spectre in N points plus
# the null frequency by repeatedly applying the fft to overlapping chunks of
# 2*N points.
# Faster if N is a power of 2.  The variance of each element in the computed
# vector is approximately equal to avg*N/length(X), where avg is the
# element's average.  Gives a vector of length N+1.

if (length(x) < 2*n)
endif
A = 0;
iterations = max (1, floor(length(x)/n)-1);
for i = 0:iterations-1
a = fft (x(i*n+1:(i+2)*n));
A = A + abs(a(1:n+1)).^2;
A = A + abs(flipud(fliplr(postpad(a(n+1:2*n), n+1)))).^2;
endfor
A = A / n^2;
endfunction

```