Estimation de fréquence
Principe
The goal of frequency estimation is, given N successive observations of a pure sinusoidal or exponential tone, contaminated by AWGN noise, to estimate as accurately as possible the frequency of the signal. Note that the two problems (real and complex cases) are not equivalent.
These scripts (see links at the right of this page) contain the SCILAB implementation of various estimators, based both on time series analysis and on DFT interpolation. You can also find in the attached PDF note (link at the right) a short description of the available estimators and the results I obtained by simulating them.
Installation
First, if not already done, install a recent version of Scilab on your computer (official SCILAB download page). Then, download the two scripts (links at the right of this page) and execute them inside Scilab (to do this, just open them with the Scilab editor Scinotes, and press F5). Now you can use the estimation functions as described below.
Script API
estimator = get_freqestim(name)Build an estimator object. Supported estimators are the following (see attached PDF file for a short description of each one):
'kc', 'lrp', 'pscfd', 'argmax(dft)', 'quadratic',
'Gaussian', 'McEachern', 'Jacobsen', 'Candan',
'Candan 2', 'Cedron', 'Cedron cplx'.
Note that the name specification is not case sensitive.
nf = freqestim(estimator, x)Estimate the normalized frequency nf (nf = f/fs), using the specified estimator and the input signal x.
test_freqestim(lst,cols,form)Test a list of frequency estimators: plots the mean square error versus SNR for each estimator. lst must be a list of estimators objects, cols is a vector of color specification (same as plot third argument), and form specify real ('r') or complex ('c') test.
Example 1: complex tone frequency estimation
fs = 1e3; // 1 KHz sampling frequency
f = 0.21 * fs; // F = 0.21 Fs = 210 Hz
N = 128; // 128 samples
// Get the Jacobsen estimator, designed for complex tones
estimator = get_freqestim('Jacobsen');
// Generate the test signal
t = linspace(0,(N-1)/fs,N);
x = exp(2*%pi*%i*f*t);
// Proceed to estimation
f = freqestim(estimator, x)
disp(f)
0.2100002
Example 2: real tone frequency estimation
fs = 1e3; // 1 KHz sampling frequency
f = 0.21 * fs; // F = 0.21 Fs = 210 Hz
N = 128; // 128 samples
// Get the Cedron estimator, designed for real tones
estimator = get_freqestim('Cedron');
// Generate the test signal (real signal this time)
t = linspace(0,(N-1)/fs,N);
x = sin(2*%pi*f*t);
// Proceed to estimation
f = freqestim(estimator, x)
disp(f)
0.21
Example 3: evaluation of different estimations algorithm
This example shows how to use the function test_freqestim so as to compare different algorithms (here Gaussian, Candan 2, Cedron and PSCFD estimators, for complex tones).
// Build the list of frequency estimators
lst= list(get_freqestim('Gaussian'),get_freqestim('Candan 2'),...
get_freqestim('Cedron'), get_freqestim('PSCFD'));
// Curve colors: red+stars, cyan+diamonds, magenta, blue+circles
cols = ['r-*','c-d','m-','b-o'];
// Test and plot MSE = f(SNR), complex case
test_freqestim(lst, cols, 'c')
Here is the figure obtained:
Below are more complete graphs, both for complex and real tones detection:
