On using eigenvectors to estimate a signals' fundamental frequency, via MUSIC

Question

Context:

(Disclaimer: This is NOT a comm problem).

I am trying to estimate the fundamental frequency of a real, periodic signal. This signal, was constructed by match filtering a raw signal, to that of a pulse. (the matched filter). The resultant signal has the following characteristics:

• It is periodic. (Fundamental is 1/period), and this is what I am trying to estimate.

• It is non-stationary in time. Specifically, the amplitudes of the periodic pulses can vary in amplitude. (e.g, one pulse can be low, while another is high, and the next low again, and one after that medium, etc).

• I believe it is stationary in frequency, (in so much as you accept changing amplitudes, but not changing bands).

• It has harmonic distortion. What I mean here is that, (and correct me if I am wrong), but that the individual pulses within the signal are not sinusoids, but are 'funky' shapes like a gaussian, triangle-ish, half-a-parabola, etc.

I am trying to estimate the fundamental frequency of this signal.

Of course, sometimes the raw signal is nothing but noise, but it still goes through the path and gets matched filtered anyway. (More on that later).

What I have tried:

Now, I am aware of a multitude of fundamental frequency estimators such as

1. The auto-correlation method
2. YIN, and all its dependencies
3. FFT method.

etc,

• YIN: I have not tried YIN yet.

• FFT Method: The FFT method will give you all the harmonics and the fundamental, but I have noticed that it can be finicky especially with this non-stationary business, as the fundamental is not always the highest peak. Very quickly, you find yourself trying to ascertain which of the many peaks is the fundamental, and it becomes a hard problem.

• Autocorrelation: The autocorrelation method seems to do better than the FFT method, but it is still sensitive to the amplitude irregularities of the time-domain signal. The auto-correlation method measures the distance between the center lobe, to the next highest lobe. That distance corresponds to the fundamental. However in non-stationary cases, this secondary lobe can be way too low, and you might miss it in some thresholding scheme.

It then occurred to me that maybe I can use a subspace method like MUSIC to estimate the fundamental. Upon testing this, I found that it really does give some very nice results - it peaks - robustly - and even in non-stationary cases - at frequencies corresponding to the fundamental of your signal. (Set the number of signals you are looking for to 2, and it will retrieve the fundamental - i.e, pick the 2 highest eigenvectors (corresponding to the highest values of the eigenvalues) of the signals' covariance matrix, discard them, and construct the noise subspace from the remaining, project your hypothesis complex sinusoids against them, take the reciprocal, and voila, a nice pseudo-spectrum).

Questions & Problems:

1. That being said, I would still like to understand why this works better.
2. In MUSIC we discard the signal subspace and use the noise subspace. It seems to me that the eigenvectors of the signal subspace actually are some sort of 'best fit' - they are in fact optimal matched filters. So: Why not simply use the signal subspace eigenvectors directly? (I know its not MUSIC anymore but why is using noise subspace better then?)
3. Lastly, the final problem is that although this method seems to work much more robustly for non-stationary signals (as defined above), the problem is that now I am ALWAYS getting an answer - even when there is nothing but noise in the system! (I mentioned above that the raw pre-matched filtered signal can be just white noise sometimes, when you do not have a periodic signal present).

What ways might exist to counteract this? I have tried looking at the eigenvalues and there is some more 'curvature' in their decay in cases where there is just noise VS cases where there is a signal, but I fear it might not be robust enough.

Bonus:

4. When are the eigenvectors of a covariance matrix sinusouds VS something else? What determines whether or not they are sinusoids or not? Why arent they square waves? Or insert-other-shape-here signals?

Answer 1

The autocorrelation matrix is diagonalized by sinusoids when the process is stationary, this follows from the fact that the covariance operator is a convolution for a stationary process. A more rigorous proof is that $$f(t,s)=Cov(X(t),X(s))=Cov(X(t-u),X(s-u))=f(t-u,s-u)$$ which in particular means that $f(t,s)=f(t-s,0)$ which is also a positive seimidefinite function of $t-s$, hence by Bochner's Theorem we have $$Cov(X(s),X(t))=\int_{-\infty}^{\infty}e^{i(s-t)x}d\mu(x)$$ which proves the claim.

The intuition is that an autocorrelation matrix estimated for some finite set of observations in a signal asymptotically behaves likes a circulant matrix because the correlation depends only on the time differences not the absolute positions and circulant matrices have discrete sinusoids as their eigenvectors (since they are convolution operators). There are plenty of proofs of this and this is a sketchy intuition.

The set of autocorrelation functions that are diagonalized by sinusoids are exactly those that correspond to stationary processes, but many other processes's autocorrelation functions will be approximately diagonalized by sinusoids over some interval. These processes correspond to those that can approximated by stationary processes over an interval. More details are here.

General non-stationary processes can have autocorrelation functions that need not be diagonalized by sinusoids.

Locally stationary processes will either have a slowly changing spectrum and/or a small number of well-spaced abrupt changes in the spectrum. Speech, animal noises, music, and many other natural sounds fit this description. The reason why subspace identification algorithms work, as I understand it, is that some form local stationarity (not rigorous) generally holds for the types of signals we analyze.