Estimating number of sinusoid from spectral estimation

Suppose I am given a data X0 that have 50 realizations of a sum of an unknown number of sinusoids in additive white Gaussian noise.

I plotted the Spectral estimate as attached. How do I estimate the number of sinusoids in X0, their frequencies, and their relative amplitudes? is it just the peak? Edit:

n1 = 1; n2 = length(x)
Px = abs(fft(x(n1:n2),1024)).^2/(n2-n1+1);
Px(1)=Px(2); • can you maybe plot a single signal in time domain and it's PSD estimate, so that one gets a chance to assess it? Regarding your plot, all I dare to say with certainty is that it is, indeed, very colorful and nice to look at. Mar 9 '18 at 17:27
• also, the exact method of your spectral estimate would be interesting; what we're looking at is distinctively not line spectra, which one would ideally see from observing sums of sinusoids. Mar 9 '18 at 17:31
• @MarcusMüller I used a for loop for each realizations and used periodogram to plot out the resulting spectrum. I can average the total to make a single PSD estimate as attached Mar 9 '18 at 18:07
• how do you do that "single PSD estimate"? Mar 9 '18 at 20:55
• @marcus by summing each realization and averaging it Mar 9 '18 at 21:04

The best approach and the details of how to do it depend greatly on the particulars of your situation. If all your vectors are of the same length, and all your sinusoids have an integer number of cycles per vector (and therefore their frequencies correspond to a DFT bin) you can take the DFT of each vector. Add the vectors coherently if the phases are the same in each random realization, add the DFT magnitudes if the phases are not the same. Then set a threshold and find the points above the threshold.

If the frequencies are not on DFT bins but are spaced apart from each other, you need to modify the above approach to take account of the fact that the energy of each sinusoid is spread over multiple DFT bins. You can still use threshold-based detection but will need to interpolate to find the frequencies.

If the frequencies could be spaced very closely together - which, from your plot, appears to be the case - the DFT-based approach is difficult. The best approach may be a subspace-based technique like ESPRIT (Estimation of Signal Parameters Via Rotational Invariance Technique). This involves estimating the covariance matrix and computing its eigendecomposition, then using the eigenvalues to estimate the power of each sinusoid and the eigenvectors to estimate the frequencies. See, for examples: Richard Roy, Thomas Kailath: "Estimation of Signal Parameters via Rotational Invariance Techniques." IEEE Transactions on Acoustics, Speech, and Signal Processing. Vol 37, No 7, July 1989.

• The 50 vectors are of same length of 128. I have read the piece and I'm not entirely sure on how to go about the ESPRIT method. Is it possible to have a rough estimate of the frequency just by taking the frequency at the peaks?(0.18pi and 0.21pi) Also, is it safe to assume that there are 2 unknown sinusoids from the 2 peaks? Mar 9 '18 at 18:31
• Yes, you can certainly estimate the frequencies from the two peaks in the DFT. You can also estimate the amplitudes this way. The one complication is that it looks like the true frequencies may lie between DFT bins, so with a short window of 128 samples the two tones blur together in the frequency domain. The closer they get, the harder your task will be. If you need higher accuracy that one DFT bin, then interpolation may be useful. As for whether there are 2 unknown sinusoids, that seems likely unless there is another tone whose frequency is very close to the other two. Mar 9 '18 at 18:44
• Thank you, that helped a lot. Is the general idea for estimating the number of unknown sinusoids the number of peaks there are? Mar 9 '18 at 18:54
• Yes, exactly. It's simple in principle but sometimes complicated in practice, especially if you are trying to automate the process. The complications arise when (1) the frequency does not lie exactly on a DFT bin, as is generally the case, in which case you'll have to deal with the sin(x)/x spectrum of the windowed sinusoid; and (2) the frequencies of the two sinusoids are close enough that the spectra blur together. Mar 9 '18 at 20:29
• Would you consider the peak at 10dB as another peak? Mar 9 '18 at 21:05

Smoothing your data before the DFT will eliminate a lot of the noise and make your results better. This is particularly true for sinusoids with lower frequencies.

Also, for frequencies that fall between bins, the DFT values on the two bins will be nearly opposite each other on the complex plane. This indicator is lost when you are looking at the magnitudes only. This effect may be masked by the leakage from nearby tones (sinusoids) in the DFT.

Hope this helps a bit.

Ced