Spectrum estimation with noisy data, unbiased estimation of amplitudes needed

[omg!]

hi everyone,

i would be grateful for suggestions for the following problem:
i have a set of measurements, discrete samples from a continuous function of time. with each measurement comes a estimate of its uncertainty (standard error). by doing a discrete Fourier transform, i have seen some frequency content in the amplitude spectrum that is distinctly visible above the noise. however, due to the significant noise, the amplitude of the peaks cannot be taken as the real amplitude of those oscillations. do you happen to be familiar with a method that can be used to estimate the true amplitude of the oscillations without bias? the gained precision in amplitude estimation may come with the sacrifice of frequency resolution, which is not that important to me.

i have already found some papers concerning spectrum estimation, but most of them use fairly complicated algorithms that i don't feel like implementing. I'm wondering if there is a basic solution that i have overlooked. thanks for any input! cheers

 

[berkeman]

hi everyone,

i would be grateful for suggestions for the following problem:
i have a set of measurements, discrete samples from a continuous function of time. with each measurement comes a estimate of its uncertainty (standard error). by doing a discrete Fourier transform, i have seen some frequency content in the amplitude spectrum that is distinctly visible above the noise. however, due to the significant noise, the amplitude of the peaks cannot be taken as the real amplitude of those oscillations. do you happen to be familiar with a method that can be used to estimate the true amplitude of the oscillations without bias? the gained precision in amplitude estimation may come with the sacrifice of frequency resolution, which is not that important to me.

i have already found some papers concerning spectrum estimation, but most of them use fairly complicated algorithms that i don't feel like implementing. I'm wondering if there is a basic solution that i have overlooked. thanks for any input! cheers

I am not an expert in this, but it seems like one approach would be to use the noisy DFT to identify discrete frequencies of interest, and then go back and filter the time domain data for each particular frequency of interest (digital bandpass filters), and then for each of those filtered datasets, do a DFT to get a better idea of the actual amplitude of each frequency component of interest...

 

[marcusl]

I am not an expert in this, but it seems like one approach would be to use the noisy DFT to identify discrete frequencies of interest, and then go back and filter the time domain data for each particular frequency of interest (digital bandpass filters), and then for each of those filtered datasets, do a DFT to get a better idea of the actual amplitude of each frequency component of interest...

Actually you can show formally that the DFT is equivalent to a bank of linear phase digital filters with adjacent passbands, so there's usually little to be gained by iterating in the way you outline.

do you happen to be familiar with a method that can be used to estimate the true amplitude of the oscillations without bias? the gained precision in amplitude estimation may come with the sacrifice of frequency resolution, which is not that important to me.

i have already found some papers concerning spectrum estimation, but most of them use fairly complicated algorithms that i don't feel like implementing. I'm wondering if there is a basic solution that i have overlooked. thanks for any input! cheers

Look into Welch's method of averaged and windowed periodograms (see Oppenheim and Schafer for a detailed treatment). It is the standard "non-parametric" technique for computing power spectral density using the FFT. Little is assumed about the nature of the signals or noise, so it is straightforward and pretty forgiving. Dealing with small signals close to the noise level, however, is a tough problem and you may not do very well with this approach (or any other for that matter). Unfortunately things get complicated from here.

If you know something about the signal, then you can do much better with parametric spectral estimators (fancy words for modeling). For example, if you know that three continuous sinusoids are present in the noise, then you can use one of several auto-regressive (AR) methods of order three. One of the best performing, and most complex, is Burg's maximum entropy spectral estimator. If you know the characteristics of the noise, then that can be useful in a model also. Estimators are available for other situations, too, like if the system has both poles and zeros. The drawbacks of parametric approaches are a) they are complicated and b) they work well when your assumptions about the signals are correct and poorly otherwise (AR performs badly if, for instance, you have broad spectral features instead of sinusoids, or non-stationary signals where the frequency or some other parameter changes with time).

Finally if you know a lot about your signal of interest, then you can implement a matched filter. This is used in situations like radar where the shape of the returned echo is known (since it's a copy of the waveform that you transmitted a short time earlier). A matched filter will dig out a signal that is buried in noise, if you know precisely what you are looking for.

 

[berkeman]

Very interesting. Thanks marcus.

 

[omg!]

thank you very much for the suggestion to use a band pass filter. i have tried it, but it merely filters out the band passed part of the spectrum without changing the amplitudes (see attachment). the data was generated with sines with amplitudes 0.37 and 0.24, added with normally distributed random numbers of std dev 1.

DFT seems to produce amplitudes that are slightly biased. the mean amplitude of the 0.37*sin term of 10000 repetitions is 0.3756 with std deviation 0.0439. i can live with the slight bias, so i guess what I'm really looking for is a method that gives me unbiased or slightly biased amplitude estimates with lower standard error than DFT.

i will definitely look into welch's method, thanks for the suggestion. it's also very convenient since it is already implemented and ready to use in matlab.

i do know a lot about the noise. the noise in every data point is normally distributed (with known variance), but the noise terms happen to be slightly correlated with each other. maybe this could be a problem too?
unfortunately, i don't know anything about the signal. in fact, i wish there was no signal at all! you see, my problem is to check whether signals are present, and in the unfortunate case where there are, i need to estimate their amplitudes. i will also look into the burg method, but I'm afraid my limited knowledge of statistics will not get me very far. thank you very much for your detailed explanation, i certainly learned a lot.

 

[marcusl]

It's helpful to see your results. How long are your data sequences? Your signal components are nicely visible, so ignore some of the concerns in my earlier post.

EDIT: The FFT-based methods ("periodogram") have bias built in that you can't do much about. (I think they are asymptotically unbiased, however, meaning that the bias approaches zero for infinitely long data sequences. You'll have to check on that.) You can reduce the variance by dividing the data record into shorter records, computing the periodogram of each, and averaging. Also applying a window to each short record, and (optionally) overlapping the records further improves performance. This is Welch's procedure. There is a tradeoff since short records have poor frequency resolution, which can hurt your estimate in other ways. At the end of the day, you have to play around a little to get good results.

 

[omg!]

the graph that i showed previously used randomly generated data. data from an actual measurement is attached with this post. note that the data is not stationary, it doesn't have zero mean. i need to know whether there are vibrations "on top" of the curve, so i first subtract the moving average (window size 11), and then do the dft.

thank you very much for your help!