Take the 2-minute tour ×

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?

share|improve this question

2 Answers 2

up vote 4 down vote accepted

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.

share|improve this answer
I would add for the OP, the size of the autocorrelation matrix $R_{xx}$, as well as what exactly a 'slice' is. I do not think it is currently clear from the answer. Otherwise, its a good answer. –  Mohammad Jan 28 '13 at 15:16

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

if T==1
   signal = signal(:);
   [T,N] = size(signal);

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 (kootsoop@ieee.org)
% 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.

share|improve this answer
Peter, now that this subject is brought up, how does one choose the size of the autocorrelation matrix $R_{xx}$? For a signal of length $N$, it seems to me that the full autocorrelation matrix should always be of size $2N-1$ x $2N-1$ correct? –  Mohammad Jan 28 '13 at 15:46
As the matlab code above shows, for a single sinusoid you only need two coefficients. See Dave C's answer. –  Peter K. Jan 28 '13 at 16:37
Peter, I am asking a more general question. I know that for one real signal we need 2 co-efficients, and for one complex signal we need one co-efficient. I am asking, in general, how big one chooses to make the autocorrelation matrix when running Pisarenko or MUSIC? –  Mohammad Jan 28 '13 at 16:52
@Mohammad: My memory is that for M real-valued sinusoids, you need an autocorrelation matrix of size 2M+1. That is, 2M for the sinusoids and 1 for the "rest" (the noise). I need to go back to Kay / Marple to check for sure. Will try to review them this evening. –  Peter K. Jan 28 '13 at 17:00
Peter, thank you, I would appreciate that. –  Mohammad Jan 28 '13 at 17:13

Your Answer


