[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
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.
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
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- autocorrelation & family,
Francesco Potorti` <=