FFTfreq: possible off-by-one error for one-sided FFTs
Closed this issue · 7 comments
Could it be possible that there is some offset when computing the sample frequency by calling FftSharp.Transform.FFTfreq()?
These are the steps I've followed and the results obtained:
- Create a 1-cycle 2-Hz-sinus wave at a sampling rate of 16 Hz (see 16 points 2 Hz sinus .txt) (The offset is less noticeable the higher the sampling frequency).
- Compute the FFT using:
signalFFT = FftSharp.Transform.FFTmagnitude(signalWindowed);
freq = FftSharp.Transform.FFTfreq(sampleFreq, signalFFT.Length);
plotFFT.Plot.AddScatterLines(freq, signalFFT);- The plot shows the frequency of the signal at 1.78 Hz instead of 2 Hz.
Possible explanation:
- The array returned by
FftSharp.Transform.FFTmagnitude()has one more point than half its size:new double[input.Length / 2 + 1](2n-1 +1). - Inside function
FftSharp.Transform.FFTfreq(), the frequency is computed as:
if (OneSided)
double fftPeriodHz = sampleRate / pointCount / 2;
else
double fftPeriodHz = sampleRate / pointCount;- The problem is that
pointCountis passed as the length of the array returned byFftSharp.Transform.FFTmagnitude()(2n-1 +1), so it has one more point than what is required to computefftPeriodHz.
So, i can think of two options:
- Subtract 1 from the array's length when calling
FftSharp.Transform.FFTfreq(). This might be counterintuitive and not obvious. - Perform the subtraction inside
FftSharp.Transform.FFTfreq():
if (OneSided)
double fftPeriodHz = sampleRate / (pointCount-1) / 2;
else
double fftPeriodHz = sampleRate / (pointCount-1);This last option produces the correct result for the sinus wave mentioned at the beginning. I hesitated before doing any PR since I'm not sure which one would you consider more appropriate.
Hi @arthurits!
Regarding the two possible options, here are my thoughts:
Option 1: User supplies Length - 1
Subtract 1 from the array's length when calling FftSharp.Transform.FFTfreq(). This might be counterintuitive and not obvious.
I agree with you this would be confusing.
However, an overload that accepts the array itself could achieve this without requiring the user to have special knowledge. Adding this overload may be useful either way.
Option 2:
Perform the subtraction inside FftSharp.Transform.FFTfreq(): This last option produces the correct result for the sinus wave mentioned at the beginning.
This sounds like the way to go.
What does Python do?
When in doubt I often try to mimic what commonly used Python packages do.
It would be ideal to analyze sample data with Python's function and write a test to confirm our values match theirs exactly. If you create a PR you don't have to go through the trouble of doing all that, but I'll probably add it in before merging.
Relevant files (note to self):
- https://github.com/swharden/FftSharp/tree/main/dev/data
- https://github.com/swharden/FftSharp/blob/main/dev/python/fft.py
- https://github.com/swharden/FftSharp/blob/main/src/FftSharp.Tests/VsNumpy.cs
Pull Request
I hesitated before doing any PR since I'm not sure which one would you consider more appropriate.
If you want to create a PR, I'd be happy to work on it and merge it in shortly! If not, that's okay too, I'll fix this on my own within the next few days 👍
Thanks for your feedback @swharden! 👍
PR has been created implementing your suggestions, except the Python tests.
I realized that options 1 and 2 can't both be implemented at the same time, so I modified the PR.
Option 1: User supplies Length - 1
However, an overload that accepts the array itself could achieve this without requiring the user to have special knowledge. Adding this overload may be useful either way.
I believe the overload option to be the most practical one. This overload should compute the correct pointCount and call the original FFTfreq.
Option 2:
This sounds like the way to go.
However, the subtraction can't be implemented here because the user could call the overload passing the array which computes the correct length and calls the orginal FFTfreq function where no subtraction should be done.
In this case, the XML function documentation should clearly state what value is expected for the pointCount parameter.
Good points! I'm glad you liked the overload idea, and PR #50 is looking great
BTW I found I already implemented "official" Python/numpy fftfreq and comparison, so I think the correct behavior (or at least the common behavior shared with a really big Python library) is already locked-in
FftSharp/src/FftSharp.Tests/VsNumpy.cs
Lines 79 to 89 in 454708e
Lines 1 to 512 in 454708e
correct behavior (or at least the common behavior shared with a really big Python library) is already locked-in
Actually, this discussion has been about one-sided FFT frequency calculation. That's not the default behavior of Python/numpy, and the one-sided behavior hasn't been sufficiently tested, so it's still worth this closer inspection
Wow, that got way more complicated than I expected!
This fix modifies the math of FftFreq() only for half-size FFTs so it doesn't affect the full FFT math.
This is the new behavior. What do you think? If you like it I'll merge and publish a new package 🚀
// generate signal with 2 Hz sine wave
int sampleCount = 16;
double sampleRateHz = 16;
double sinFrequencyHz = 2;
double[] samples = Enumerable.Range(0, sampleCount)
.Select(i => Math.Sin(2 * Math.PI * sinFrequencyHz * i / sampleRateHz))
.ToArray();
// get FFT
double[] fft = FftSharp.Transform.FFTmagnitude(samples);
double[] fftKnown = { 0, 0, 1, 0, 0, 0, 0, 0, 0 };
Assert.That(fft, Is.EqualTo(fftKnown).Within(1e-10));
// calculate FFT frequencies both ways
double[] fftFreqKnown = { 0, 1, 2, 3, 4, 5, 6, 7, 8 };
double[] fftFreq = FftSharp.Transform.FFTfreq(sampleRateHz, fft.Length); // <---- original overload
double[] fftFreq2 = FftSharp.Transform.FFTfreq(sampleRateHz, fft); // <---- new overload
Assert.That(fftFreq, Is.EqualTo(fftFreqKnown));
Assert.That(fftFreq2, Is.EqualTo(fftFreqKnown));
This looks really good! Documentation of the functions is clear enough and the point about the one-sided FFT is just on the spot.
I'd say it's ready to be merged.
Nice work! 🤩