FFTfreq: possible off-by-one error for one-sided FFTs #49
Comments
|
Hi @arthurits! Regarding the two possible options, here are my thoughts: Option 1: User supplies Length - 1
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:
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):
Pull Request
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 👍 |
swharden
changed the title
|
Thanks for your feedback @swharden! 👍 |
|
I realized that options 1 and 2 can't both be implemented at the same time, so I modified the PR.
I believe the overload option to be the most practical one. This overload should compute the correct
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 |
|
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 |
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 |
swharden
changed the title
|
Wow, that got way more complicated than I expected! This fix modifies the math of 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. |
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:
Possible explanation:
FftSharp.Transform.FFTmagnitude()has one more point than half its size:new double[input.Length / 2 + 1](2n-1 +1).FftSharp.Transform.FFTfreq(), the frequency is computed as: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:
FftSharp.Transform.FFTfreq(). This might be counterintuitive and not obvious.FftSharp.Transform.FFTfreq():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.
The text was updated successfully, but these errors were encountered: