#### Estimating Frequency Over A Short Capture

by

, 2nd November 2014 at 21:26 (12691 Views)

I was recently asked how to estimate the frequency of a noisy sine wave over a short capture. In this context, 'short' means the total capture durationrelative to the period of the sine wave. So we may have, at most, a few periods of the sine wave (or, at worst, less than one period). Here is an example (download raw data: sine_data_40k.zip):

Clealy, there is less than one period available. To illustrate the difficulty with this problem, I will first describe a popular approach for estimating frequency over 'long' captures.

FFT-based Frequency Estimation ('Long' Captures)

Whenever we are interested in the frequency content of a signal, the Fast Fourier Transform (FFT) is often an excellent tool to use. In practice, when working with real-world (finite and imperfect) data, rather than using the raw output of the FFT, it is often preferable to enhance the FFT's estimate of Power Spectral Density (PSD) using Welch's method, since this can significantly reduce noise (at the expense of reduced resolution).

Here is some Matlab code showing how to use Welch's method to estimate frequency over a 'long' capture (download: pwelch_test.zip):

Code Matlab M - [expand] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 close all; clear all; clc; % Assume we capture Nsamps samples at 1kHz sample rate Nsamps = 32768; fsamp = 1000; Tsamp = 1/fsamp; t = (0:Nsamps-1)*Tsamp; % Assume the noisy signal is exactly 123Hz fsig = 123; signal = sin(2*pi*fsig*t); noise = 1*randn(1,Nsamps); x = signal + noise; % Plot time-domain signal subplot(2,1,1); plot(t, x); ylabel('Amplitude'); xlabel('Time (secs)'); axis tight; title('Noisy Input Signal'); % Choose FFT size and calculate spectrum Nfft = 8192; [Pxx,f] = pwelch(x,gausswin(Nfft),Nfft/2,Nfft,fsamp); % Plot frequency spectrum subplot(2,1,2); plot(f,Pxx); ylabel('PSD'); xlabel('Frequency (Hz)'); grid on; % Get frequency estimate (spectral peak) [~,loc] = max(Pxx); FREQ_ESTIMATE = f(loc); title(['Frequency estimate = ',num2str(FREQ_ESTIMATE),' Hz']);

The code outputs this figure:

Although the input looks like a noisy mess, we can see an extremely narrow spike in the spectrum at 123.0469 Hz (i.e. a tiny error of 0.0469 Hz).

However, if we try the same approach with the 'short' data capture (using this code: short_pwelch.zip), we get:

Now the frequency spike is invisible because it is squeezed right down in the first bin after DC. The reason is that, even though I have eliminated any averaging in the Welch method, the frequency bin spacing (i.e. resolution) is still only:

Therefore, even though the input signal is somewhere in the region of 35kHz, it is impossible for the FFT to provide an answer anywhere between 0 and 50kHz.

It is worth reiterating that the 'short' capture is short in terms of its total duration compared to the period of the sine wave. In fact, there are significantlymoresamples in the 'short' capture than the 'long' capture, but this doesn't ultimately help.

Frequency Estimation Via Curve Fitting

A far more effective method with 'short' captures is curve fitting. In a previous blog, I wrote about polynomial curve fitting. Well, in this case, we want to do sinusoidal curve fitting. This is easy to do in Matlab, using the 'sin1' option with the fit() function (code: short_sine_fit.zip):

Code Matlab M - [expand] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 close all; clear all; clc; % Load data load sine_data_40k.dat x = sine_data_40k; % Data info Tsamp = 5e-10; fsamp = 1/Tsamp; Nsamps = length(x); t = (0:Nsamps-1).'*Tsamp; % Get sine fit F = fit(t, x, 'sin1'); % Plot data and sine fit figure(); plot(F,t,x); grid on; xlabel('Time (secs)'); ylabel('Amplitude'); title(['Estimated frequency = ', num2str(F.b1/(2*pi)), ' Hz']);

The code outputs this graph:

The frequency estimate of 35526.7 Hz is very much more accurate than the FFT-based method.