FFT and estimating frequencies not on discrete points

Assume there is a real function f and its spectrum F. If f is discrete
sampled, we have samples at some t = st * k, where st is sampling
period and k = 0, 1, 2, ..., N - 1. Also assume every frequency in f
is within [-fc..fc] range, fc = Nyqist frequency. If f is a simple
sine function whose minimums, zeros and maximums are at the sampling
points, F will be very "good", without leakage and just by looking at
some point of F we can see what is the frequency m of f =
sin(2*pi*m).

However, if this is not the case, i.e. extremal points of f do not lay
on the sampled points (k = 0, 1, 2, ..., N - 1), this will not be the
case - there will be leakage. Assume that I calculated F from the
known f and I have N/2 bins of F at my hand. Bin n - 1 is for
frequency m - j, bin n is for frequency m and bin n + 1 for frequency
m + j. How can I calculate the power for frequency e.g. m - j / 2 or m
+ j / 3 or such? That is, calculate the power at some frequency for
which I don't have F defined, which I presume would allow me to see
which frequency is the "main" one (let away special cases).

wimxa@yahoo.com wrote:

> Assume there is a real function f and its spectrum F. If f is discrete
> sampled, we have samples at some t = st * k, where st is sampling
> period and k = 0, 1, 2, ..., N - 1. Also assume every frequency in f
> is within [-fc..fc] range, fc = Nyqist frequency. If f is a simple
> sine function whose minimums, zeros and maximums are at the sampling
> points, F will be very "good", without leakage and just by looking at
> some point of F we can see what is the frequency m of f =
> sin(2*pi*m).

> However, if this is not the case, i.e. extremal points of f do not lay
> on the sampled points (k = 0, 1, 2, ..., N - 1), this will not be the
> case - there will be leakage.

You can do an M function least squares fit for M <= N.

> Assume that I calculated F from the
> known f and I have N/2 bins of F at my hand. Bin n - 1 is for
> frequency m - j, bin n is for frequency m and bin n + 1 for frequency
> m + j. How can I calculate the power for frequency e.g. m - j / 2 or m
> + j / 3 or such? That is, calculate the power at some frequency for
> which I don't have F defined, which I presume would allow me to see
> which frequency is the "main" one (let away special cases).

You lose any advantage of the FFT doing this.  For N functions that are
different from the basis functions of the FFT there is an N by N matrix
that you can multiply the result of the FFT by to get the result you want.

-- glen

>If f is a simple
>sine function whose minimums, zeros and maximums are at the sampling
>points, F will be very "good"

The "minimums, zeros, maximums" requirement is misleading (although not
exactly wrong). I can get a clean peak out of an FFT on a sinewave that
you wouldn't recognize by eye as samples of a sine wave.
For example, take "3" in the following example:
http://www.elisanet.fi/mnentwig/webroot/FFT_interpolation_how_does_it_work/index.html

sines/cosines with a full number of periods within the FFT lenght are
orthogonal to each other. Any frequency "in-between" is not.

For arbitrary frequencies, you may have some luck using the techniques
below, read the instructions (regarding windowing etc):
http://www.dspguru.com/howto/tech/peakfft.htm

On Sep 23, 9:26 am, wi...@yahoo.com wrote:
> Assume there is a real function f and its spectrum F. If f is discrete
> sampled, we have samples at some t = st * k, where st is sampling
> period and k = 0, 1, 2, ..., N - 1. Also assume every frequency in f
> is within [-fc..fc] range, fc = Nyqist frequency. If f is a simple
> sine function whose minimums, zeros and maximums are at the sampling
> points, F will be very "good", without leakage and just by looking at
> some point of F we can see what is the frequency m of f =
> sin(2*pi*m).
>
> However, if this is not the case, i.e. extremal points of f do not lay
> on the sampled points (k = 0, 1, 2, ..., N - 1), this will not be the
> case - there will be leakage.

There won't be any leakage for any sinusoid which is strictly
periodic in your N point DFT aperture, even if the extrema
points do not line up with the sample instants.  Any phase
sinusoid will work, as long as it is strictly periodic.  It
will show up as some ratio between the real and imaginary
components of one DFT bin.

> Assume that I calculated F from the
> known f and I have N/2 bins of F at my hand. Bin n - 1 is for
> frequency m - j, bin n is for frequency m and bin n + 1 for frequency
> m + j. How can I calculate the power for frequency e.g. m - j / 2 or m
> + j / 3 or such? That is, calculate the power at some frequency for
> which I don't have F defined, which I presume would allow me to see
> which frequency is the "main" one (let away special cases).

Bin frequencies in a DFT are orthogonal.  The representation
of bin frequencies only is unique.  The same is not true for
non-bin frequencies, so you may have to make some assumptions
about the lack of nearby interference or noise in your
estimation of any non-bin frequencies.

A DFT of length N can be looked at as a DFT of a rectangular
window of length N on your signal.  In the frequency domain,
this is equivalent to convolving the DFT of your frequency
with a Sinc function.  So, in the absence of interference,
you can try fitting a pair of mirrored Sinc functions to the
bins nearest your non-bin frequency of interest (for more
than 2 points, do the fit in the complex domain).  The peak
of this best-fit Sinc function will give you the amplitude.
You can also try Sinc interpolation, which may give you a
reasonable approximation.



IMHO. YMMV.
--
rhn AT nicholson d.0.t C-o-M
  http://www.nicholson.com/rhn/dsp.html

On Sun, 23 Sep 2007 16:26:40 -0000, wimxa@yahoo.com wrote:

>Assume there is a real function f and its spectrum F. If f is discrete
>sampled, we have samples at some t = st * k, where st is sampling
>period and k = 0, 1, 2, ..., N - 1. Also assume every frequency in f
>is within [-fc..fc] range, fc = Nyqist frequency. If f is a simple
>sine function whose minimums, zeros and maximums are at the sampling
>points, F will be very "good", without leakage and just by looking at
>some point of F we can see what is the frequency m of f =
>sin(2*pi*m).
>
>However, if this is not the case, i.e. extremal points of f do not lay
>on the sampled points (k = 0, 1, 2, ..., N - 1), this will not be the
>case - there will be leakage. Assume that I calculated F from the
>known f and I have N/2 bins of F at my hand. Bin n - 1 is for
>frequency m - j, bin n is for frequency m and bin n + 1 for frequency
>m + j. How can I calculate the power for frequency e.g. m - j / 2 or m
>+ j / 3 or such? That is, calculate the power at some frequency for
>which I don't have F defined, which I presume would allow me to see
>which frequency is the "main" one (let away special cases).

As has been mentioned, regardless of the phase of the sine wave, if
the there are an integer number of cycles in the DFT aperture you'll
get the nice peak behavior.

You seem to imply that you may be doing only real-valued processing,
though, which changes a few things but not the generalities.

If you're really trying to better estimate the frequency of a tone
between bins, there are some good resources from folks here on
comp.dsp.   Matt Donadio's treatment on dspguru has already been
linked.   Here's Peter Kootsookos' page:

http://home.comcast.net/~kootsoop/freqalgs.htm

and mine:

http://www.ericjacobsen.org/fe.htm

Eric Jacobsen
Minister of Algorithms
Abineau Communications
http://www.ericjacobsen.org

On Sep 25, 11:29 am, Eric Jacobsen <eric.jacob...@ieee.org> wrote:
> On Sun, 23 Sep 2007 16:26:40 -0000, wi...@yahoo.com wrote:
> >Assume there is a real function f and its spectrum F. If f is discrete
> >sampled, we have samples at some t = st * k, where st is sampling
> >period and k = 0, 1, 2, ..., N - 1. Also assume every frequency in f
> >is within [-fc..fc] range, fc = Nyqist frequency. If f is a simple
> >sine function whose minimums, zeros and maximums are at the sampling
> >points, F will be very "good", without leakage and just by looking at
> >some point of F we can see what is the frequency m of f =
> >sin(2*pi*m).
>
> >However, if this is not the case, i.e. extremal points of f do not lay
> >on the sampled points (k = 0, 1, 2, ..., N - 1), this will not be the
> >case - there will be leakage. Assume that I calculated F from the
> >known f and I have N/2 bins of F at my hand. Bin n - 1 is for
> >frequency m - j, bin n is for frequency m and bin n + 1 for frequency
> >m + j. How can I calculate the power for frequency e.g. m - j / 2 or m
> >+ j / 3 or such? That is, calculate the power at some frequency for
> >which I don't have F defined, which I presume would allow me to see
> >which frequency is the "main" one (let away special cases).
>
> As has been mentioned, regardless of the phase of the sine wave, if
> the there are an integer number of cycles in the DFT aperture you'll
> get the nice peak behavior.
>
> You seem to imply that you may be doing only real-valued processing,
> though, which changes a few things but not the generalities.
>
> If you're really trying to better estimate the frequency of a tone
> between bins, there are some good resources from folks here on
> comp.dsp.   Matt Donadio's treatment on dspguru has already been
> linked.   Here's Peter Kootsookos' page:
>
> http://home.comcast.net/~kootsoop/freqalgs.htm
>
> and mine:
>
> http://www.ericjacobsen.org/fe.htm
>
> Eric Jacobsen
> Minister of Algorithms
> Abineau Communicationshttp://www.ericjacobsen.org

This one in particular is good, although he doesn't handle windows.

[1] Macleod, M.D., "Fast Nearly ML Estimation of the Parameters of
Real or Complex Single Tones or Resolved Multiple Tones," IEEE Trans.
Sig. Proc. Vol 46 No 1, January 1998, pp141-148

Just wondering...
If I go back to square one, it may turn out that FFT itself is a bad
approach. 
Simple carrier recovery / PLL / FM demodulator / ... might be what is
really needed.

-mn

On Sep 25, 12:10 pm, "mnentwig" <mnent...@elisanet.fi> wrote:
> Just wondering...
> If I go back to square one, it may turn out that FFT itself is a bad
> approach.
> Simple carrier recovery / PLL / FM demodulator / ... might be what is
> really needed.
>
> -mn

Bad in what way?  Maximum likelihood is hard to beat.

mnentwig wrote:
> Just wondering...
> If I go back to square one, it may turn out that FFT itself is a bad
> approach. 
> Simple carrier recovery / PLL / FM demodulator / ... might be what is
> really needed.

Why FM demod? Once the PLL is locked up, read it's frequency. If the 
frequency is variable, the PLL *is* the FM demod.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

>>Bad in what way? Maximum likelihood is hard to beat.
Huh? What does FFT have to do with maximum likelyhood?

From a practical point of view, I doubt think that a majority of carrier
recovery schemes relies on FFT, and that may be what the OP really needs.

>Why FM demod? Once the PLL is locked up, read it's frequency. If the 
>frequency is variable, the PLL *is* the FM demod.

That's right and that's my point. But search for PLL alone, and I doubt
something relevant will come up at the top. Search for FM demodulation
-and- PLL... well, you said it already.

-mn