0
$\begingroup$

I generated a time history of 200 points using the equation $$x(t) = A\cos(2\pi f t)+ B\sin(2\pi f t) $$

with $A = 1$, $B = 0.1$ and $f = 17.2 \mathrm{Hz}$ and a sample rate of 400 samples per second. The numerical Fourier transform of the time history has a modulus plot of:

Mathematica graphics

The question is how to get the frequency, amplitude and phase from this plot? This is the typical case where the frequency falls between bins. As there are three unknowns I would imagine you could get the precise values from, say, the three largest points which I call premax, max and postmax. However, complex values are available for all the points so many methods are possible.

This question arose as an output from a previous question where it was suggested I ask afresh and that several answers were available. I have also asked this question on Mathematica Stackexchange here where it was suggested I ask on Signal Processing.

Many answers are available. Exact rather than approximate answers are preferable. Some, like using a Prony series or curve fitting a sine wave, do not use the Fourier transform. These are also interesting but I would like to use the Fourier data.

If there is noise on the data then some sort of fitting is probably needed depending on the extent of the noise. All suggestions welcome. Thanks

Edit

A bit more in response to the answers which I still have to fully investigate

There is an exact analytic formula for the Fourier transform of a sampled sine wave. For my data it plots as

Mathematica graphics

This plot could also be obtained by calculating in-between values for the original data. I also calculated the frequency at the maximum of the analytic curve. To 25 significant figures it is

17.22726210702718945079653

which is significantly different to value of 17.2 used in the calculation of the time history.

Some of the answers suggest using parabolic interpolation. I show the results of this here

Mathematica graphics

I have tried parabolic interpolation using the complex values of the Fourier transform and the absolute values. The absolute values show a better maximum than the complex values but they are still well off. I still have other suggestions to try.

Methods based on time domain analysis look more promising than using Fourier which is disappointing.

$\endgroup$
  • $\begingroup$ Hugh, in general it is frowned upon to cross-post the same question across multiple stack exchange sites... though that question you asked on Mathematica.SE is quite old, so I don't think it's strictly cross-posting! :-) $\endgroup$ – Peter K. Sep 23 '15 at 17:09
  • $\begingroup$ Peter K. I thought this was reasonable since the issue was raised again yesterday and it was also suggested on Mathematica. Hope this is not an issue but if you insist I will remove it. $\endgroup$ – Hugh Sep 23 '15 at 17:13
  • $\begingroup$ I wrote the note so that other's won't complain about cross-posting. I think your question here is fair game. If I were that worried I just would have closed it. $\endgroup$ – Peter K. Sep 23 '15 at 17:14
  • $\begingroup$ That's very odd. The parabolic interpolations are among the fastest and most accurate. Have you tried using some of the code on Eric's website? I suspect there is something wrong with your implementation or something amiss with your data. $\endgroup$ – Peter K. Sep 24 '15 at 11:31
  • $\begingroup$ @PeterK. I have just put a quadratic through three points. Is it more subtle? If you look you can see it goes through the points. $\endgroup$ – Hugh Sep 24 '15 at 13:10
2
$\begingroup$

For a lower quality frequency estimate, but still more precise than just using the (center) frequency of the peak magnitude DFT result bin, try using parabolic interpolation (using 3 neighboring points) of the spectral peak.

You can't get the phase from a DFT magnitude-only plot. For phase you will need an estimator that uses the complex DFT results.

For a higher quality estimate, including a phase estimate, try using windowed-Sinc interpolation of the complex DFT result components, perhaps followed by successive approximation to find a better estimate of the magnitude peak near the peak DFT bin. The quality of the estimate will vary with the width of the window.

Equivalent to Sinc interpolation would be to use a much longer zero-padded DFT or FFT to get more closely spaced frequency points.

There are a few added comments regarding the above frequency estimation methods on my blog: http://www.nicholson.com/rhn/dsp.html#1

Note that near DC or Fs/2 (closer than the width of your interpolation kernel), you may need to use a mirrored interpolation kernel. Note that hybrid methods are possible: e.g. using Sinc interpolation to up the resolution by 4X or 10X, then using parabolic or higher-degree polynomial interpolation to estimate the peak location, then using Sinc interpolation to estimate the magnitude and phase at that estimated peak frequency. etc. Note that sinusoidal phase estimation in noise will often be more accurate when the phase is referenced to the center of the DFT window (N/2), not to sample point 0. Referencing phase to the center of a data vector can be done using some form of fftShift.

The question does not ask about frequency estimators that do not involve a DFT, but note that in zero noise, only 3 or 4 unaliased points are required to solve for 3 unknowns (frequency, magnitude and phase of a pure unmodulated sinusoid). See: http://www.claysturner.com/dsp/3pointfrequency.pdf and http://www.claysturner.com/dsp/4pointfrequency.pdf

$\endgroup$
  • 2
    $\begingroup$ Hey! Nice references to Clay's stuff. :-) $\endgroup$ – Peter K. Sep 23 '15 at 20:56
  • $\begingroup$ @hotpaw2 Thanks. Your blog is a very extensive list of methods. Most helpful. Again I note that time domain methods and not DFT based methods seem to be attractive. The Clay results are new to me but very simple. $\endgroup$ – Hugh Sep 24 '15 at 8:23
1
$\begingroup$

Eric Jacobsen has done quite a bit of work on this topic.

Check out his first and second frequency estimation pages.

There are also some algorithms here that may be of interest.

Let me know if any of these are of any interest, and I can elaborate.

$\endgroup$
-2
$\begingroup$

Finding the peak frequency

Your frequency resolution $\Delta f$ is directly related to the amount of samples and the aquisition frequency like this: $\Delta f = \frac{f_A}{N}$ with $N$ being the amount of samples of your signal, and $f_A$ being the sample aquisition frequency.
Thus, to increase your frequency resolution you can either decrease the aquisition frequency (Make sure to not drop below the nyquist frequency of your signal though), or increase the amount of samples.
There are 2 ways to accomplish the latter:

  1. Increase Aquisition Time
    This is pretty obvious. If you observe the signal for a longer period, you get more information about the signal, but aquiring more samples isn't always possible.
  2. Zero-Padding
    By adding zeros to your signal you can artifically elongate your aquisition time, thus generating more samples and lowering $\Delta f$. Zero-padding your signal doesn't change the information content, and is thus the most precise way to find the peak frequency.

Phase and Amplitude

Your image shows a magnitude plot of your DFT. This means that for each complex value $Y_k$ of the DFT only $\left|Y_k\right|$ is plotted. Assuming that you actually do get complex values in mathematica, and this is just the way the plotting method operates, here's what you do:
You can obtain the amplitude of a given frequency-value $Y_k$ by taking the absolute $\left|Y_k\right|$. To get the phase, you get the complex argument $\arg{Y_k}$.

$\endgroup$
  • $\begingroup$ Thank you for your practical comments with which I agree. However, if you use either of your methods you will still end up with a cluster of points around the frequency of interest. Certainly, the resolution will be better but you still have to do something with the premax, max and postmax to get the frequency and amplitude. This final postprocessing is my concern which perhaps you did not spot. $\endgroup$ – Hugh Sep 24 '15 at 8:16
  • $\begingroup$ -1: Your "Little Tip on the side" is incorrect. FFT algorithms can be done for any prime radix signal length, not just $2^N$. And, in fact, most common libraries that perform the FFT are not limited to $2^N$ signal lengths. $\endgroup$ – Peter K. Sep 24 '15 at 11:29
  • $\begingroup$ I agree with you Peter K. $\endgroup$ – Lukas Sep 24 '15 at 18:17
  • $\begingroup$ To support Peter K.'s assertion, have a look at <techonline.com/electrical-engineers/education-training/…>. $\endgroup$ – Richard Lyons Sep 25 '15 at 12:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.