Pisarenko Harmonic Decomposition
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.