I'm attempting to code Pisarenko's method for some time-series data that I have.
I've had a look at a couple of papers, which explain in detail the theory behind it, but I'm having some trouble figuring out how to convert the math into an algorithm for either C or Python.
Is there any information on how one can do that? A pseudocode algorithm, or some hints on which direction to start etc?
UPDATE: On choosing the size of $R_{xx}$: If you suspect there are $M$ sources, then you need your autocorrelation matrix ($R_{xx}$) to be made up of M+1 unique (i.e. non-overlapping) snapshots. For this reason, the snapshots technically have to be less than len(x) / M+1. Estimating $R_{xx}$ is not trivial, and there are some pitfalls for certain signals and certain selections of N (size of autocorr matrix) if the segments or slices are non-overlapping. This is why I typically use as many segments of x as possible when estimating $R_{xx}$
First, you need an autocorrelation matrix. One of the more fool-proof ways to do that, is to add up successive outer products of "slices" of the input vector where a slice is x[idx:idx+N].
# x contains input vector , usually a time-series
N = 20 # length of one side of the square autocorr matrix
Rxx = NxN matrix of zeros
foreach unique slice of x of length N:
s = unique slice
Rxx += outer product of (s,s) #outer(s,s) in Numpy
Rxx /= nsnap # divide by number of snapshots that Rxx is composed of
Second, calculate the eigenweights and eigenvectors of $R_{xx}$. Save the "minimum" eigenvector. The minimum valued eigenweight and minimum eigenvector should correspond to the subspace of the noise.
w,v = eig(Rxx)
vmin = v[:,abs(w).argmin()]
At this point, you have $v_{min}$ and $e^H$ is simply
eh = exp( 1j*omega*arange(N) )
Where $PHD = \frac{1}{ | e^H v_{min} |^2}$ or
phd[omegaIdx] = 1.0 / ( abs(eh * vmin) ** 2 )
Do that for all omegas of interest. The non-pseudo-code is in Python/Numpy, hence a * b is an element-by-element multiplication.
This probably doesn't do exactly what you're after, but I put this Matlab code together many moons ago and think it may help in understanding what's happening.
function omegahat = pisarenko(signal)
%PISARENKO Find the Pisarenko frequency estimates of the signals
% in each column of signal.
%
% omegahat = pisarenko(signal)
%
% [1] V. F. Pisarenko, "On the Estimation of Spectra by
% Means of Non-linear Functions of the Covariance Matrix,"
% Geophys. J. Roy. astr. Soc., Vol. 28, pp. 511-531, 1972.
%
% [2] V. F. Pisarenko, ``The Retrieval of Harmonics from a
% Covariance Function," Geophys. J. Roy. astr. Soc.,
% Vol. 33, pp. 347-366, 1973.
%
% $Id: pisarenko.m,v 1.2 1999/01/09 11:10:35 PeterK Exp PeterK $
% File: pisarenko.m
%
% Copyright (C) 1999 Peter J. Kootsookos
%
% This software provided under the GNU General Public Licence, which
% is available from the website from which this file originated. For
% a copy, write to the Free Software Foundation, Inc., 675 Mass Ave,
% Cambridge, MA 02139, USA.
%Type cast it to double.
signal = double(signal);
%--------------------------------------------------------
% Error condition checks
[T,N] = size(signal);
if any([T N]==0)
error('pisarenko: zero size data not allowed.');
end;
if T==1
signal = signal(:);
[T,N] = size(signal);
end;
Rss1 = sum(signal(2:T,:).*signal(1:T-1,:));
Rss2 = sum(signal(3:T,:).*signal(1:T-2,:));
alpha = ( Rss2 + sqrt(Rss2.^2 + 8*Rss1.^2) ) ./ ( Rss1 + eps ) / 2;
omegahat = acos(alpha/2);
% Author: Peter J. Kootsookos ([email protected])
%
% Based on: P.J. Kootsookos, S.J. Searle and B.G. Quinn,
% "Frequency Estimation Algorithms," CRC for Robust and
% Adaptive Systems Internal Report, June 1993.
%
UPDATE: Regarding the issue of model order selection (size of the autocorrelation matrix to use), this turns out to be more complex than my comments below suggest.
Marple suggests on page 374 that
The idea of separating eigenvectors into signal and noise subspaces based upon an examination of either the eigenvalues of the autocorrelation matrix or the singular values of the data matrix does not work well in practice, especially with short sample records. The AIC (Akaike Information Criterion) order-selection criterion first introduced in Chap. 8 has been extended by Wax and Kailath [1985] to handle the subspace separation problem.
