Estimating period of low frequency oscillations: autocorrelation vs. Frequency approaches

Question

I have a noisy signal around $15\text{s}$ long sampled a $32\,\text{Hz}$.
I am trying to estimate the peak frequency or period for a low frequency component with expected frequency in between $0.08\,\text{Hz}$ and $0.2\,\text{Hz}$.

Empirically I find easier to find such a peak in the time domain by measuring the lag at peak autocorrelation, rather than by identifying peaks in DFT derived representations.
Is there any way I can justify this from a theoretical perspective?

Answer 1

If we can assume white noise and the signal itself is of constant frequency over the full 15 second duration, then the autocorrelation would be suboptimal in determining the frequency compared to the DFT. Given the relatively narrow frequency range, I would recommend neither DFT nor autocorrelation but to do a correlation with frequency tones directly as a maximum likelihood estimator, where we can then resolve to any resolution desired (as limited by SNR) using:

$$Corr = \sum_{n=0}^{N-1}x[n]e^{-j2 \pi n f/f_s}$$

Where $x[n]$ is the noisy waveform, $f$ is a test frequency in Hz for correlation, and $f_s$ is the sampling rate in Hz.

Note how this looks very similar to the DFT, with the difference being the use of any arbitrary frequency for $f$. So instead of correlating a noisy signal with a time delayed copy of itself (the auto-correlation), we correlate the noisy signal with a noise-free copy of what we are looking for (a pure tone). Repeat the correlation with tones stepped by the desired precision and the maximum likelihood estimate will be at the maximum correlation result. Ultimately the post-correlation SNR will limit the ability to discern a single peak when the precision exceeds what the SNR will allow, which results in a statistical uncertainty in the frequency estimate consistent with the post correlation SNR.

To see why the autocorrelation approach would be suboptimal in comparison, consider the OP's case of a 15 second duration signal with a 0.2 Hz repetition hidden in the noise. With 0.2 seconds there are three repetitions and we expect a correlation at every 5 second offset. The autocorrelation at offset $\tau = 0$ would be the correlation using all samples in the 15 second duration, providing the maximum correlation over any of the other peaks. However, the next correlation peak with the 5 second offset would only have a 10 second waveform overlap and therefore does not use all the samples to maximize the processing gain or SNR for that next peak.

In contrast, the correlation formula given above is using every sample to test for each $f$ and therefore (under the conditions stated) would provide the optimum SNR for estimation. The DFT would also provide an optimum SNR but would be limited by the granularity of the DFT bin size, or require further interpolation using the result for multiple bins within vicinity of the peak.

Further, note the significance that the noisy signal is correlated with a noise-free reference (exponential tone) when using the suggested correlation, in contrast to correlating the noisy signal with itself with the autocorrelation which increases noise or reduces SNR in the result. The noise increase is quantified (for each multiplication in the correlation) statistically at MattL's answer here.

Answer 2

Similar approach to Dan's but a different way to go about it. First lets define what exactly we mean by "peak" frequency. I suggest it is the frequency that minimizes the square error between the original signal and the sine wave of that frequency. The error is

$$E = \sum \left( A\cdot \cos(\omega t + \varphi) - x[n] \right)^2 $$

and it can be minimized with respect to amplitude $A$, frequency $\omega$ and phase $\varphi$

Unfortunately this is a non-linear minimization problem for $\omega$ and $\varphi$ . However you can easily calculate the partial derivatives and then solve it iteratively with steepest decent (or a similar method). So this would roughly look like this:

1. Get a first estimate $A_0$, $\omega_0$ and $\varphi_0$ from a DFT
2. Calculate current error $E[i]$
3. Calculate partial derivatives.
4. Determine $A_{i+1}$ directly (the partial derivative is linear in A)
5. Estimate $\omega_{i+1}$ and $\varphi_{i+1}$ from the derivatives and the current error
6. Go back to step 2: repeat until the error or the partial derivative become "small enough" or movement of the frequency is less than your desired precision.

Answer 3

I'd just add to these two fantastic answers that since the frequencies of interest are so low compared to your sampling frequency, I imagine you can also decimate your signal quite a bit before trying either methods.

---

EDIT: Implementing @DanBoschen's answer
Here is the result of a simulation with the OP's specifications:

It's very robust to SNR (try increasing the noise power!)

And code to reproduce:

clear all; close all;

% signal parameters
fs = 32e3;
T = 15;
N = fs*T;

% time vector
t = (0:1/fs:T-1/fs);

% target and probe frequencies
target_freqs = [0.08,0.12,0.13,0.18,0.2];
probe_freqs = (0.06:0.005:0.3);

% Will add noise to our signal
noise = randn(1,N);

figure(1)
hold on;
cnt = 1;
for k = target_freqs

    % signal at low frequency + noise
    x = sin(2*pi*k*t);
    sig = x+noise;
    
    % option: decimate 
    R = 40000; fsnew = fs/R; Nnew = N/R;
    sig = decimate(sig,R);
    
    % window
    sig = hamming(Nnew)'.*sig;

    %Dan Boschen
    tnew = (0:1/fsnew:T-1/fsnew);
    Wmat= exp(-1i*2*pi*probe_freqs'*tnew);
    Corr = abs(sig*Wmat');
    [~,indmax] = max(Corr);
    temp = stem(probe_freqs,Corr,'.');
    hold on
    pltLegend(cnt) = stem(probe_freqs(indmax), Corr(indmax), 'linewidth', 2, 'Color', get(temp, 'Color'));
    cnt = cnt+1;

end
grid on
xlabel('probe frequencies (Hz)')
ylabel('Corr magnitude')
legendStrings = "f_{target} = " + string(target_freqs);
legend(pltLegend,legendStrings);

Answer 4

Empirically I find easier to find such a peak in the time domain by measuring the lag at peak autocorrelation, rather than by identifying peaks in DFT derived representations. Is there any way I can justify this from a theoretical perspective?

The biggest difference between Fourier transform and autocorrelation is that autocorrelation will include any harmonic frequencies, while Fourier transform only includes the single frequency sine wave.

As such, Fourier will yield better results if you have a pure sine wave disturbed by noise, as it discards any noise at different frequency. Autocorrelation works well if your signal of interest is periodic, but not purely sinusoidal - the higher frequency harmonics can provide more accuracy to the calculation.

As a comparison, here are the Fourier transform (using FFT) results for a sine wave and a square wave signal:

In the case of a pure sine wave, the autocorrelation peak is quite smooth and easily affected by noise. For square wave, the autocorrelation peak is sharp and probably more precise than the frequency obtained from FFT results.

The FFT still works well, because square wave has a strong fundamental frequency component. If the signal was more impulse-like, such as an ECG heart beat signal, results are different:

In autocorrelation results, the peak at 628 sample lag is still readily distinguishable. But in FFT results, there is no peak at f = 26. This is because the impulse signal in this case has very little energy at the fundamental frequency, but has strong harmonic components.