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

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:

1. 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?
• Mohammad- Can you please make a few edits/clarifications? I can be a stickler for terminology, but it's important for future visitors. In addition to 'nice & clean', can say harmonic distortion. Instead of repetitive, can you say periodic. Stationary can refer to time-varying statistics or time-varying spectrum. Can you clarify? The autocorrelation method is an alias for the Yule-Walker method. When you say 'number of signals' is this real sinusoids or complex exponentials? Can you use largest value eigenvalue? Rank has other meanings in linear algebra. Same with 'highest variance'... Commented May 27, 2012 at 16:39
• ...(cont) One important thing (and I'll note this in my answer when you clarify), is that the MUSIC method is a noise subspace method. So, ideally, the signal subspace eigenvectors, the ones with the largest value eigenvalues, are not used. Also, your signal is a sum of sinusoids if it is periodic. If it's periodic, it can be defined by a Fourier series, which is a sum of discrete sinusoids. Commented May 27, 2012 at 16:41
• @Bryan Sorry for the delay in getting back (long weekend), I will actually revamp the entire question soon and let you know - thanks! Commented May 29, 2012 at 18:23
• @Bryan I have finally revamped the entire post, added your suggestions, and also clarified a lot of the context/problem. Please see. By all means let me know if I can clarify anything else. Commented Jun 1, 2012 at 15:44
• @Mohammad Can you discern whether a signal is present or not by the "strength" of the eigenvectors- i.e. the eigenvalues? Commented Jun 1, 2012 at 17:23

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.

• The answers are for everyone's benefit. It's worth nothing that Bochner's theorem is a generalization of the more familiar Wiener-Khinchin. $\mu$ is the spectral density.
– Emre
Commented May 31, 2012 at 20:09
• @MarkS Thank you a lot. I have some follow ups: 1) Based on this, can we say that a process is stationary in as much as the eigenvectors of its covariance matrix are sinusoidal? Can this be a sort of measure of stationarity? 2) You mention "...and circulant matrices have discrete sinusoids as their eigenvectors (since they are convolution operators)..." I am unclear as to what this means - what operators? Can you please clarify. 3) When you say "The set of autocorrelation functions" you are talking about the rows of the covariance matrix? Thanks again. Commented Jun 1, 2012 at 16:14
• @Mohammad Cheers: 1) Yes, this can be loosely thought of as a measure of stationarity. 2) A circulant matrix is formed from all cyclic permutations of a vector, so multiplying a circulant matrix by another vector is a convolution between those two vectors. 3) An autocorrelation function Corr(s,t) is the autocorrelation between X(s) and X(t) for some random process X. I call it a function because I am wanting to simultaneously handle the continuous and discrete case. The sample autocorrelation matrix can be seen as a discrete approximation to this function. Commented Jun 1, 2012 at 19:42
• @Emre thanks for pointing out the Wiener–Khinchin_theorem, I learned my Fourier analysis first on groups and was never formally introduced to it in a signal processing class. Commented Jun 1, 2012 at 19:48