Possible fine-detail error in Matlab fft?

13 vues (au cours des 30 derniers jours)
Tran Trung Luu
Tran Trung Luu le 15 Mar 2021
Commenté : David Freiman le 19 Fév 2022
Dear all,
I often work with FFT and I like the fft function from Matlab a lot:
However, recently I encoutered a possible problem with the current implementation of Matlab:
In the standard calculation of the frequency axis that I extracted from the above link:
f = Fs*(0:(L/2))/L;
Fs is the sampling frequency and L is the length of the signal, it seems that the frequency step (the distance between two consecutive frequencies in the frequency axis) is Fs/L. This is usually not a problem if L is large enough, for example in this case. However, when I deal with small L, for example L = 5, then I realize that the frequency step is actually not the inverse of the total time window, which it should theoretically be. If I want to do this theoretically perfect, I had to use Fs/(L-1) but then this is not compatible with how the fft is done in Matlab, i.e. I got inconsistency of frequency peaks' position. It seems to me that the fft was not implemented with "perfect theoretical consideration".
Could you please advise/comment? Or you think I missed something?
Thanks and best regards,
Trung

Réponses (1)

Prabhan Purwar
Prabhan Purwar le 21 Mar 2021
Hi,
As the frequency resolution is defined as Fs/L, both the implementation and example seem rightly scaled and working as expected. If you can provide the calculations and steps you have tried, documents you have referred it will be helpful in resolving and understanding the issue more clearly.
Thanks
  2 commentaires
Tran Trung Luu
Tran Trung Luu le 25 Mar 2021
Hi Prabhan,
Thanks for your message.
I fully understand that "both the implementation and example seem rightly scaled and working as expected". I said this is a potential "fine-detail" problem. The problem is not big, but the implementation is probably not perfect, and that's why I want to discuss to find the perfect solution.
Here is a more detailed description:
Theoretically speaking, fft in Matlab is the fast DFT implementation of the Fourier Transform. From the intrinsic formalism for Fourier Transform, we know that the frequency step (df) is reciprocal of the time window (T) and vise versa, the frequency window (F) is the reciprocal of the time step (dt). If you consider now a discretized time axis or frequency axis, they are reprensented by N points, from t_1 to t_N, and the time step is dt = t_2 - t_1. Similarly we have f_1 to f_N and df = f_2 - f_1;
The standard description of frequency step/resolution in Matlab is Fs/L. This is not really wrong, but let's consider the following case: if we have only three points, N = 3. In this case, if you take L = N = 3, then the calculated fft result matches well the frequency grid calculated using Fs/L. However, it is physically wrong. The reason is that, in this case, the time window is (N-1)*dt which is only 2dt, not 3dt. So the actual frequency step, which it should theoretically be, should be 1/(2dt). But the real one calculated with the given formula is 1/(3dt).
In short, the discrepancy comes from the small number of the samples. If L is more than 1000, for example, then it does not matter much if we do Fs/1000 or Fs/999. For engineering purposes, it is ok in most cases. However, for scientists like me, it is important to make sure it is perfect.
Let me know how you think.
Thanks and best regards,
Trung
David Freiman
David Freiman le 19 Fév 2022
Hi Trung,
I recently had a similar question and came across this post. I believe the resolution to this issue lies in the fact that for a discrete time signal as you describe above, with N=3, the time window is actually 3 * dt, because each sample point represents a time of dt.
In the past I have made the mistake of calculating the time window as t(N) - t(1). As you pointed out above, this would barely be noticeable for a large N of, say, 1000. However, when recently checking some code with a very low number of samples, N, I realized that the time window should be calculated as:
T = N * dt = N / Fs
So delta F would be 1 / T = Fs / N
Please let me know if you have any questions or see any errors in my logic.
Thank you,
Dave

Connectez-vous pour commenter.

Tags

Produits


Version

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by