You could use the Gaussian ratio technique for frequency estimation:
http://fftguru.com/fftguru.com.tutorial4.pdf
and associated C++ code
But your low frequency creates a problem. As noted in the classic 1978 windows paper by Harris (http://web.mit.edu/xiphmont/Public/windows.pdf), nearby negative frequencies can add constructively/destructively with positive frequencies. Thus, it's important not to operate too closely to either the 0 or N/2 frequency bins - and with inputs from .45 to 4 Hz, you are definitely doing that.
So the first thing to do is to modulate your input up to say, 20 Hz or so. Since I don't usually work on comms stuff on a day to day basis, I dusted off one of my old comms books and read up about ssb. As you only have a single input sinusoid (with an additional 90 degree shifted version), this made things a bit easier.
The following C++ code demonstrates the process: 1) generate cos/sin sinusoids of a specific frequency (note that I include pi in the time parameter below - when generating signals containing multiple sinusoids, it's easier to type pi in only once) - phi is an arbitrary phase variable that can be changed; 2) ssb modulate the signals up to 20 Hz using the phasing method - the addition operation chooses the upper sideband, and you get a real signal as a result - ssb gets things away from the 0 frequency bin; 3) compute the Gaussian window; 4) multiply the data record by the window; 5) compute FFT and print out results; 6) compute frequency and amplitude estimates.
// *****************************************************************************************
// ssb_freq_est - 1) generate a sinusoidal signal between .45 Hz and 4 Hz, 2) bandshift up
// by 20 Hz and obtain upper sideband by ssb, 3) multiply data record by a Gaussian window,
// 3) take fft, 4) locate the highest amplitude squared output(s), 5) compute the frequency
// of the tone, written by Kevin J. McGee at fftguru.com
#include <cstdio> //
#include <cstdlib> //
#include <iostream> //
#include <cmath> //
using namespace std; //
void fft_recur(long, double* r, double* i);
int main (int nNumberofArgs, char* pszArgs[ ] ) {
const long N = 64;
double sample_rate = 64., r[N] = {0.}, i[N] = {0.}, s[N] = {0.}, w[N] = {0.} ;
long n ; double phi = 0., alpha = 6., t, twopi = 6.2831853071795865;
double fbin = 0., amp = 0., c;
for (n = 0; n < N; n++) { // generate input data
t = twopi*n/sample_rate;
r[n] = 10.*cos(1.5*t + phi) ;
i[n] = 10.*sin(1.5*t + phi) ;
} // end for
// create ssb - upper sideband
for (n = 0; n < N; n++) { // modulate up to 20 Hz
t = twopi*n/sample_rate;
r[n] = r[n]*sin(20.*t) ;
i[n] = i[n]*cos(20.*t) ;
r[n] = r[n] + i[n]; // + chooses upper sideband
} // end for
//Gaussian window computed using HARRIS 78 windows paper
for (n = -N/2; n < N/2; n++) {
w[n + N/2] = exp(-2.*(n*alpha/N)*(n*alpha/N)) ;
} // end for
// multiply data record by Gaussian window
for (n = 0; n < N; n++) {
r[n] = w[n]*r[n] ;
i[n] = 0. ;
}
fft_recur(N, r, i);
for (n = 0; n < N/2 + 1; n++) {
r[n] = r[n]*r[n] + i[n]*i[n] ; // compute magnitude squared
printf("%2d\t%9.7f\n",n,r[n]);
}
cout << "\n\n\n";
// compute frequency and amplitude
// in code below, c is a parameter related to the Gaussian window beamwidth
// - fbin is a fractional bin number (ie: if input sinusoid is at 23.5 Hz, then fbin should
// compute to be 23.5 - actual frequency is then fbin*sample_rate/N )
for (n = 16; n < 26; n++) { // compute for bins 16 to 25
c = 8.*(alpha/twopi)*(alpha/twopi);
fbin = n +.5 +.5*c*(.5*log(r[n + 1]) - .5*log(r[n]));
amp = sqrt(c)*exp(.5*log(r[n]) + ((fbin-n)*(fbin-n))/c)*(2.0*sqrt(twopi/2.)/N) ;
printf("%2d\t%9.7f\t%9.7f\n",n,fbin*sample_rate/N,amp);
}
system ("PAUSE");
return 0;
} // end main
//******************** fft_recur ***********************
void fft_recur(long N, double *r, double *i) {
long h, i1, j = 0, k, i2 = N/2, l1, l2 = 1;
double c = -1.0, s = 0.0, t1, t2, u1, u2;
for (h = 0; h < N-2; h++) { // ***** bit reverse starts here ***
if (h < j) {
t1 = r[h]; r[h] = r[j]; r[j] = t1;
t2 = i[h]; i[h] = i[j]; i[j] = t2;
}
k = i2;
while (k <= j) {
j = j - k; k = k/2;
}
j = j + k;
} //****** bit reverse done ******
for (k = 1; k < N; k = k*2) {
l1 = l2; l2 = l2*2;
u1 = 1.0; u2 = 0.0;
for (j = 0; j < l1; j++) {
for (h = j; h < N; h = h + l2) {
i1 = h + l1;
t2 = (r[i1] - i[i1])*u2 ;
t1 = t2 + r[i1]*(u1 - u2) ;
t2 = t2 + i[i1]*(u1 + u2) ;
r[i1] = r[h] - t1;
i[i1] = i[h] - t2;
r[h] = r[h] + t1;
i[h] = i[h] + t2;
} // end for over h
t1 = u1 * c - u2 * s;
u2 = u1 * s + u2 * c;
u1 = t1; //x = u1 - u2; y = u1 + u2;
} // end for over j
s = - sqrt((1.0 - c) / 2.0);
c = sqrt((1.0 + c) / 2.0);
} // end for over k
} // end fft *******************************************************
With an input of 1.5 Hz, and ssb'd by 20 Hz, the program computes the following fft output (bins 10 to 30):
10 0.0000000
11 0.0000000
12 0.0000001
13 0.0000112
14 0.0008969
15 0.0416534
16 1.1179498
17 17.3406472
18 155.4460465
19 805.3144096
20 2411.1413024
21 4172.0668447
22 4172.0668527
23 2411.1412963
24 805.3144131
25 155.4460450
26 17.3406477
27 1.1179496
28 0.0416535
29 0.0008969
30 0.0000112
You can see that the output shows something equidistant between bins 21 and 22, which is where it should be for a 1.5 Hz/20 Hz ssb signal. The frequency/amplitude calculation results in (bins 16 to 25):
16 21.5000001 10.0000015
17 21.5000000 9.9999998
18 21.5000000 10.0000000
19 21.5000000 10.0000000
20 21.5000000 10.0000000
21 21.5000000 10.0000000
22 21.5000000 10.0000000
23 21.5000000 9.9999999
24 21.5000000 10.0000002
25 21.5000001 9.9999991
Since the peak(s) in the previously shown FFT output were at bins 21 and 22, the frequency/amplitude calculation of bin 21 above is the one to use (the bin 21 output used bins 21 and 22 in the calculation). As you can see, it's quite accurate, as you would expect from a good algorithm and zero input noise.
The above code was not at all optimized - there are a lot of efficiencies that could be done (eg: use LUTs instead of sin/cos functions, pre-compute some variables, etc). And depending on the amount of noise you have, things may not be optimal for your problem (ie: you may need to use more points in the data record to meet your accuracy requirements, or ssb bandshift differently, or shift to a different part of the spectrum, or choose a different Gaussian beamwidth, or replace the FFT with a faster, more accurate one, etc.
But I've already spent too much time on the above code, and I don't know what kind of problems you may encounter.
EDIT: Corrected the stated bin numbers (22 changed to 21). Also, run the code with very low frequencies (eg: .0000001), and you'll see just how good this technique can be.
