Mistake in Python implementation of Pisarenko's harmonic decomposition

Question

I'm trying to implement Pisarenko's Harmonic decomposition in Python and can't get a proper pseudo spectrum revealing the three frequencies of the signal I'm analysing. I tried a few things (using np.cov instead of the manual autocorrelation computed with a matrices product, transpose without conjugate, checking every array dimension) without spotting the remaining mistake(s). My script, which computes the pseudo-spectrum for each individual eigenvector instead of considering only eigvecs[:, np.argmin(np.abs(eigvals))]:

import scipy
import numpy as np
import matplotlib.pyplot as plt

######################################################################################################
# CREATE TRIPLE FREQUENCY SIGNAL AND PLOT IT
freq_sampl = 1e6

freqs = np.array([freq_sampl/10, freq_sampl/5, freq_sampl/4])
hyp_nbr_freqs = 3

sampl_time = 0.001
t = np.arange(0, sampl_time, step=1/freq_sampl)

sigs = np.exp(1j*2*np.pi*freqs[0]*t)
sigs += np.exp(1j*2*np.pi*freqs[1]*t)
sigs += np.exp(1j*2*np.pi*freqs[2]*t)

noise_std = 0.1
noise = (np.random.normal(loc=0.0, scale=noise_std, size=sigs.shape) +
         1j*np.random.normal(loc=0.0, scale=noise_std, size=sigs.shape))
sigs += noise

ax1 = plt.subplot(411)
plt.plot(t,sigs)
plt.title("Noisy signal in time domain")

######################################################################################################
# PLOT FFT TO CHECK FOR EXPECTED FREQUENCY CONTENT
sigsfft = np.fft.fftshift(np.fft.fft(sigs))
sigsfftfreq = np.fft.fftshift(np.fft.fftfreq(sigs.shape[0], d=1/freq_sampl))

ax2 = plt.subplot(412)
plt.plot(sigsfftfreq/freq_sampl, np.abs(sigsfft))
plt.title("Noisy signal in Fourier domain")

######################################################################################################
# PLOT ALL EIGENVECTOR PSEUDO-SPECTRUMS
nbr_sampls = sigs.shape[0]
ar_order = hyp_nbr_freqs + 1 # determines dimensions of autocorr matrix used for eigendecomposition

# accumulate shifted ar_order-sized vectors of successive samples
sampl_matrix = sigs[0:ar_order][np.newaxis, :]
for start_idx in range(1,nbr_sampls-ar_order):
    sampl_matrix = np.concatenate((sampl_matrix, sigs[start_idx:start_idx+ar_order][np.newaxis, :]))

autocor = np.matmul(np.transpose(np.conjugate(sampl_matrix)), sampl_matrix) / start_idx
eigvals, eigvecs = scipy.linalg.eigh(autocor)

freqs = np.array([idx * freq_sampl/1000 for idx in range(500)]) # list of freq to evaluate pseudo spectrum at
def steervec_func(freqs, ar_order):
    ar_order = np.arange(ar_order)
    return np.exp(1j*2*np.pi*freqs[:,np.newaxis]*ar_order[np.newaxis,:])

steervec = steervec_func(freqs, ar_order)
for eig_idx in range(eigvals.shape[0]):

    eig_vec = eigvecs[:, eig_idx] # plotting all of them instead of just np.argmin(np.abs(eigvals)) eig vector
    pseudo_spect = np.squeeze(np.matmul(np.conjugate(steervec), eig_vec[:,np.newaxis]))

    ax3 = plt.subplot(413)
    plt.plot(freqs / freq_sampl, 1/np.abs(np.squeeze(pseudo_spect))**2)
    plt.title("Pseudo-spectrums (steervec * eigvec)")

    # use FFT of zero-padded eigvec for alternative proposed by Baddioes
    zpad_eigvec = np.pad(eig_vec, pad_width=(0, pseudo_spect.shape[0] - eig_vec.shape[0]))
    pseudo_spect_fft = np.fft.fft(zpad_eigvec)
    pseudo_spect_fft_freqs = np.fft.fftfreq(zpad_eigvec.shape[0], 1/freq_sampl)

    ax4 = plt.subplot(414)
    plt.plot(pseudo_spect_fft_freqs, 1/(np.abs(pseudo_spect_fft) ** 2))
    plt.title("Pseudo-spectrums (FFT padded eigvec)")
plt.show()

This returns "flat" pseudo spectrums for each eigenvector when I rely on np.matmul(np.conjugate(steervec), eig_vec[:,np.newaxis]); one of the latter (the one corresponding to the min. eigenvalue, i.e. eigvecs[:, np.argmin(np.abs(eigvals))]) should however reveal the three individual sine waves within the signal:

This implementation should be equivalent to what was proposed in this dsp.SE answer and corresponds to this other dsp.SE answer. I was also a bit surprised to see a "weird" formula ( Rss2 + sqrt(Rss2.^2 + 8*Rss1.^2) ) associated with the same method in Peter K.'s answer but I didn't investigate much in that direction, assuming the wikipedia pseudo-spectrum equation should work:

$$ Pseudo-spectrum(\omega) = f(2*\pi*f) = \frac{1}{|e^H(\omega) v_{min}|^2} $$

EDIT: -------

As Baddioes said in his answer:

The reason that Pisarenko's method works is that the assumption is the covariance matrix has dimensionality one greater than the signal subspace, so that the entire noise subspace is contained in one eigenvalue.

so I want to emphasize we have here 3 sine waves in sigs and a 4x4 autocorrelation matrix autocor due to ar_order = hyp_nbr_freqs + 1, hence the 4 pseudo-spectrums on the last plot, one of which being expected to reveal the 3 frequencies since one eigenvector is expected to represent the noise subspace.

EDIT 2: -------

Added a second set of pseudo spectrums (last graph on updated image), one per eigvec, this time computed with Baddioes' suggestion i.e. the FFT of the zero-padded (right-side) eigvec:

# use FFT of zero-padded eigvec for alternative proposed by Baddioes
zpad_eigvec = np.pad(eig_vec, pad_width=(0, pseudo_spect.shape[0] - eig_vec.shape[0]))
pseudo_spect_fft = np.fft.fft(zpad_eigvec)
pseudo_spect_fft_freqs = np.fft.fftfreq(zpad_eigvec.shape[0], 1/freq_sampl)

This produces three correctly spaced peaks though I must say I'm not sure I fully understand why this kinda works (for right-side zero-padding only ?) and why my peaks aren't located at the correct frequencies (missed a conjugate / fftshift somewhere perhaps ?).

EDIT 3: -------

Modified the last graph on the updated image with the $VV^H \rightarrow V^* V^T$ change to get the proper peaks location, which is what Baddioes suggested. This translates to the following change in my script:

# OLD
autocor = np.matmul(np.transpose(np.conjugate(sampl_matrix)), sampl_matrix) / start_idx
# NEW
autocor = np.matmul(np.transpose(sampl_matrix), np.conjugate(sampl_matrix)) / start_idx

In light of another post (cf. "non-overlapping"), I also removed the overlap in the signal slices used to compute the autocorrelation matrix:

sampl_matrix_no_overlap = sigs[0:ar_order][np.newaxis, :]
for start_idx in range(ar_order, nbr_sampls-ar_order, ar_order):
    sampl_matrix_no_overlap = np.concatenate((sampl_matrix, sigs[start_idx:start_idx+ar_order][np.newaxis, :]))

but this did not change anything in the pseudo-spectrums apparently (? is this actually necessary ? added a comment on the other post but no luck yet).

I also realized the "weird" formula ( Rss2 + sqrt(Rss2.^2 + 8*Rss1.^2) ) associated with the same method in Peter K.'s answer corresponds to the specific "Pisarenko for single tone signal" case, as described in the An exact analysis of Pisarenko's single-tone frequency estimation algorithm paper. So this point was actually off-topic here.

As the peaks have been found, although not directly with the product I initally wanted to use np.matmul(np.conjugate(steervec), eig_vec[:,np.newaxis]), I'll accept @Baddioes' already very detailed answer and open a new question as suggested to understand how the FFT of the noise zero-padded eigenvector leads to the expected peaks and is equivalent to the initial approach. If someone spots the exact mistake in my product I'm still interested in a comment though !

Answer 1

By plotting all of the eigenvectors and not including their eigenvalues, you are trying to do MUSIC, but it won't work as I'll explain.

MUSIC (and Pisarenko's by extension) solve the equation \begin{equation} a^{H}(\omega)\mathbf{G}\mathbf{G}^{H}a(\omega) = 0 \end{equation} where $\mathbf{G}$ is a matrix where each of the columns are the eigenvectors corresponding to the noise subspace, and $a(\omega)$ (equivalent to $e(\omega)$) is the steering vector. The only solutions of this equation are the true frequency values. However, if you include not only the noise subspace, but also the signal subspace, you will end up with a flat spectrum for MUSIC because MUSIC whitens the eigenvectors contained in $\mathbf{G}$, and there are no solutions that satisfy the above equation.

There are two methods for overcoming this. The first way would be to not violate the signal/noise subspace assumption of MUSIC. Your MUSIC algorithm would then be \begin{align} S_{MUSIC}(\omega) &= \frac{1}{a^{H}(\omega)\mathbf{G}\mathbf{G}^{H}a(\omega)} \\ &= \frac{1}{\sigma_{0}^{-2}\sum_{k}a^{H}(\omega)g(k)g^{H}(k)a(\omega)} \\ &= \frac{1}{\sigma_{0}^{-2}\sum_{k}\left\lvert a^{H}(\omega)g(k)\right\rvert^{2}} \end{align}

The other option would be to compute an eigenvector spectral estimate, which is similar to MUSIC but doesn't whiten the clutter subspace. Letting $\mathbf{R} = \mathbf{V}\mathbf{\Sigma}\mathbf{V}^{H}$, and $\lambda_{k}$ being the $(k,k)$ entry of $\mathbf{\Sigma}$, the (this spectral estimate is \begin{align} S_{EV}(\omega) &= \frac{1}{a^{H}(\omega)\mathbf{G}\mathbf{\Sigma}_{G}^{-1}\mathbf{G}^{H}a(\omega)} \\ &= \frac{1}{\sum_{k}\lambda_{k}^{-1}a^{H}(\omega)g(k)g^{H}(k)a(\omega)} \\ &= \frac{1}{\sum_{k}\lambda_{k}^{-1}\left\lvert a^{H}(\omega)g(k)\right\rvert^{2}} \end{align} Since $\mathbf{R} = \mathbf{V}\mathbf{\Sigma}\mathbf{V}^{H}$, Eigenvector is a hybrid between Capon and MUSIC. Including all eigenvalues results in an approximate Capon method, but practically may have slightly better performance than Capon. Eigenvector will have slightly worse performance than MUSIC in scenarios where the model is matched, but will be much less sensitive to spurious estimates when there is a model mismatch (due to the more semi-parametric relationship with Capon).

The reason that Pisarenko's method works, ie \begin{equation} S_{Pisarenko}(\omega) = \frac{1}{\left\lvert a^{H}(\omega)v_{min}\right\rvert^{2}} \end{equation} is that the assumption is the covariance matrix has dimensionality one greater than the signal subspace, so that the entire noise subspace is contained in one eigenvalue. This is practically pretty hard to guarantee. Additionally, the statistical accuracy of MUSIC increases with increasing covariance matrix size, so we can assume Pisarenko's method will have far worse statistical accuracy than MUSIC in most practical applications.

EDIT Response to Eigenvector method questions in the comments

I realized I skipped some steps in my answer, so I have edited above as well as providing this edit.

It is correct that the Eigenvector (EV) method uses only the noise subspace, like MUSIC. This provides MUSIC like performance for sinusoid detection, but also gives better noise spectrum shaping since there is no whitening of the clutter. However, the methodology in getting to EV differs from MUSIC. I'll try to explain.

To understand EV, we have to understand Capon. Capon wishes to solve \begin{align} &\min_{h}E\left[\left\lvert h^{H}x\right\rvert^{2}\right] \; \text{s.t.} \; h^{H}z = 1 \\ & \min_{h}h^{H}\mathbf{R}h \; \text{s.t.}\;h^{H}z = 1 \end{align} The solution to this form is \begin{equation} h = \frac{\mathbf{R}^{-1}z}{z^{H}\mathbf{R}^{-1}z} \end{equation} producing the spectrum (after plugging in for h) \begin{equation} S = \frac{1}{z^{H}\mathbf{R}^{-1}z} \end{equation}

Many spectral analysis textbooks contain this derivation, so I won't provide it here. You can also see Capon's original IEEE paper (I believed published in 1969). The Stoica and Moses book also gives several different derivations. The point is that, in capon, we say $z=a(\omega)=\begin{bmatrix}1 & e^{j\omega} & \cdots & e^{j\omega(m-1)}\end{bmatrix}$. In EV, we wish to find some optimal matrix $\mathbf{C}$ such that $z=\mathbf{C}a(\omega)$. As the paper linked in the comments explains, this optimal matrix is the sum of the outer product of the eigenvectors in the noise subspace.

So, let's say we have an $N$-by-$N$ correlation matrix, and $n$ sources. We would then say \begin{align} \mathbf{C}_{opt} = \sum_{k=1}^{N-n}v(k)v^{H}(k) \end{align} This produces the spectrum (as shown in the paper) \begin{align} S_{EV}(\omega) &= \frac{1}{a^{H}(\omega)\mathbf{C}_{opt}\mathbf{R}^{-1}\mathbf{C}_{opt}^{H}a(\omega)}\\ &= \frac{1}{\sum_{k=1}^{N-n}\lambda_{k}^{-1}\left\lvert a^{H}(\omega)v(k)\right\rvert^{2}} \end{align}

So, while EV looks very similar in form to MUSIC, it is a completely different derivation. MUSIC again attempts to find solutions to the aforementioned quadratic form. EV attempts to find the optimal weighting of $z$ to maximize the spatial resolution of a source bearing within the Capon framework.

This has huge ramifications for the resulting spectra produced by EV and MUSIC, and is discussed at the end of the paper. Practically, we often don't know $n$, and approximate it with $n_{est}$. Briefly, if $n_{est} < n$, MUSIC will produce a spectrum with only $n_{est}$ peaks, whereas EV will produce a spectrum with $n$ peaks but with non-optimal resolution. Particularly, if $n_{est}=0$, EV produces the Capon method, i.e., the worst you will do with EV is Capon. If $n_{est}>n$, MUSIC will produce erratic peaks, due to its tendency to want to produce a spectrum with $n_{est}$ peaks. EV spectra with too large $n_{est}$ tend not to vary too much from spectra produced from $n_{est}=n$. The downside of EV is that it does not whiten the clutter spectrum, and so produces a typical noise spectrum shape. While this may often be desirable, in extremely low SNR cases, it may be more difficult to detect the peaks out of the noise if the noise spectrum isn't flat.

The reason for the relationship between Capon and EV is that \begin{align} \mathbf{R} &= \mathbf{V}\mathbf{\Sigma}\mathbf{V}^{H} \\ \mathbf{R}^{-1} &= \mathbf{V}\mathbf{\Sigma}^{-1}\mathbf{V}^{H} \end{align} The matrix $\mathbf{\Sigma}$ contains the eigenvalues of $\mathbf{R}$ on the main diagonal, and zeros everywhere else. If all eigenvalues are included ($n_{est}=0$), you will get Capon because of the definition of the eigendecomposition.