1

Been experimenting with FFT on a generated sinusoid and found something strange that doesn't seem to be described anywhere (though I may be missing something of course).

A sinusoid that exactly corresponds to a bin gives what we all know: an amplitude and a phase shift of the sinusoid in that bin.

However, for a frequency that falls between two bins the phase shift of the strongest bin expressed in polar angle seems to be shifted by almost exactly the amount by which the peak frequency is shifted from the strongest bin number. The phase shifts of the neighbouring bins also seem to show some correlation but it is not as strong.

For example, if the distance between bins is $\Delta = 44Hz$, I generate $f = 440Hz$, I get a clear sinusoid at bin $i = 10$ that yields phase angle $\phi = 0$ (the sinusoid itself is shifted by a quarter period for simplicity, hence $0$).

If I generate say

$f = 440Hz + k·\Delta$, where $k \in [-0.5, 0.5]$

then $\phi$ for the strongest bin 10 looks like $k·\pi$.

Some sample measurements:

$k=-0.5 \Rightarrow \phi=-0.499·\pi$

$k=-0.25 \Rightarrow \phi=-0.252·\pi$

$k=-0.15 \Rightarrow \phi=-0.151·\pi$

$k=-0.05 \Rightarrow \phi=-0.054·\pi$

The deviation gets bigger as you move towards the bin, though it's pretty good near the middle. Is this some kind of a parabolic correlation?

So what's the math behind this? Why do these values seem correlated?

And if they are correlated, could this be used for finding the peak frequency with some precision?

mojuba
  • 131
  • 5

1 Answers1

2

For a set of exact equations that describe how a pure real tone works in a DFT check out equations (23)-(25) in my blog article:

The phenomenon you have found can be seen clearer if you consider the complex definition of a real valued sinusiodal.

$$ \cos( \theta ) = \frac{e^{i\theta}+e^{-i\theta}}{2} $$

So, when you are near a peak bin, one of the numerator values is dominant, and you have an approximation of a complex tone. Using equation (24) in this article shows how a the distance to the bin $\delta$ rotates the bin value.

Hope this helps.

To answer your last question: Yes, I've done it different ways with more coming.

Check out:


Your $k$ is the fraction of a bin distance from a bin on the frequency scale of cycles per frame.

Thus, from (20):

$$ k = -\frac{N}{2\pi} \delta $$

The negative sign is unfortunate, I should have done it differently in the article.

$$ \delta = -\frac{2\pi}{N} k $$

Your $\phi$ is the angle of the nearest complex bin value, so I will use $\theta$ in the formula for the phase angle from the signal definition in (24):

$$ \phi = \arg(X[]) = -\delta (N-1)/2 + \theta $$

$$ \phi = \frac{2\pi}{N} k (N-1)/2 + \theta $$

$$ \phi = \frac{N-1}{N} k \pi + \theta $$

So, for N being large, and $\theta$ being zero, $\phi$ looks like $k \pi$.

Remember, this is based on using a complex tone as an approximation for the real value tone, so it won't be exact. I am not sure why your ratios are slightly larger than one when they should be slightly smaller.

Cedron Dawg
  • 7,560
  • 2
  • 9
  • 24
  • Thanks! That's a lot of information to process but I can't see an immediate answer to my question. What do you think would be the formula for peak frequency based on the phase shifts alone? There should be one apparently. – mojuba Nov 14 '19 at 20:34
  • 1
    @mojuba From (20)

    $$ f = k - \frac{N}{2\pi} \delta $$

    From (24)

    $$ \arg(X[k]) = -\delta (N-1)/2 + \phi $$

    Arg is the angle of the bin value. If you haven't been able to put these together by tomorrow, I will post the details. Remember, this is treating your signal as an approximation, which is why it isn't exact. Notice that a priori knowledge of $\phi$ is needed for this approach. This is why all the bin frequency formulas have more than one bin value used. The phase is then implicitly used.

    – Cedron Dawg Nov 14 '19 at 20:55
  • I still need some more explanation I'm afraid. So I have FFT, the frequency of interest happens to be between the bins and is known, now I want to know its phase. – mojuba Nov 15 '19 at 21:16
  • 1
    @mojuba Check out the update. To find the phase, look at my answer in https://dsp.stackexchange.com/questions/61947/amplitude-estimation-of-sinusoid-in-known-spiky-spectral-noise/61961#61961. I also strongly recommend that you implement the algorithm in https://www.dsprelated.com/showarticle/1284.php on the platform of your choice. – Cedron Dawg Nov 16 '19 at 16:12
  • Unfortunately this doesn't work. The $\frac{N-1}{N}$ correction isn't helping much, the error is too big for any practical application, even though seems close. I wonder if there's any good phase interpolation technique similar to that of the frequency. – mojuba Nov 20 '19 at 00:50
  • I suggested this method: https://dsp.stackexchange.com/a/61948/16836 which seems overkill, it uses single-frequency DFT with a window adjusted so that the frequency falls very close to a whole bin. I wonder if there's a better method... – mojuba Nov 20 '19 at 00:53
  • @mojuba The only way I know how to do it is to find the frequency first, then you find phase and amplitude from that. Goertzel is a way to calculate a single DFT bin, but you still need to know the frequency first. My answer here (https://dsp.stackexchange.com/questions/61947/amplitude-estimation-of-sinusoid-in-known-spiky-spectral-noise/61961#61961) should be instructive to you. The purpose of my answer in this question was to show you the foundation for your $k\pi$ observation. My "Two Bin Solution" article is the most efficient and accurate method I know. It is exact for a pure tone. – Cedron Dawg Nov 20 '19 at 04:24