autocorrelation & family

[Francesco Potorti`]

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.

Comments welcome.


# 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)
    x = postpad (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