diff --git a/dev/quickstart/fft-windowed.png b/dev/quickstart/fft-windowed.png index e4200cc..4704526 100644 Binary files a/dev/quickstart/fft-windowed.png and b/dev/quickstart/fft-windowed.png differ diff --git a/dev/quickstart/fft.png b/dev/quickstart/fft.png index 050e3ca..04ac119 100644 Binary files a/dev/quickstart/fft.png and b/dev/quickstart/fft.png differ diff --git a/dev/quickstart/periodogram.png b/dev/quickstart/periodogram.png index a3d37f0..60ac273 100644 Binary files a/dev/quickstart/periodogram.png and b/dev/quickstart/periodogram.png differ diff --git a/src/FftSharp.Tests/FftFreqTests.cs b/src/FftSharp.Tests/FftFreqTests.cs new file mode 100644 index 0000000..03e5dbf --- /dev/null +++ b/src/FftSharp.Tests/FftFreqTests.cs @@ -0,0 +1,78 @@ +using NUnit.Framework; +using System; +using System.Linq; + +namespace FftSharp.Tests +{ + internal class FftFreqTests + { + [Test] + public void Test_FftFreq_KnownValues() + { + // https://github.com/swharden/FftSharp/issues/49 + + // 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); + double[] fftFreq2 = FftSharp.Transform.FFTfreq(sampleRateHz, fft); + Assert.That(fftFreq, Is.EqualTo(fftFreqKnown)); + Assert.That(fftFreq2, Is.EqualTo(fftFreqKnown)); + + ScottPlot.Plot plt1 = new(400, 200); + plt1.AddSignal(samples, sampleRateHz); + TestTools.SaveFig(plt1, "signal"); + + ScottPlot.Plot plt2 = new(400, 200); + plt2.AddScatter(fftFreq, fft); + plt2.AddVerticalLine(2, System.Drawing.Color.Red, style: ScottPlot.LineStyle.Dash); + TestTools.SaveFig(plt2, "fft"); + } + + [Test] + public void Test_Freq_Lookup() + { + /* Compare against values generated using Python: + * + * >>> numpy.fft.fftfreq(16, 1/16000) + * array([ 0., 1000., 2000., 3000., 4000., 5000., 6000., 7000., + * -8000., -7000., -6000., -5000., -4000., -3000., -2000., -1000.]) + * + */ + + int sampleRate = 16000; + int pointCount = 16; + double[] signal = new double[pointCount]; + double[] fft = FftSharp.Transform.FFTmagnitude(signal); + + double[] freqsFullKnown = { + 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, + -8000, -7000, -6000, -5000, -4000, -3000, -2000, -1000 + }; + double[] freqsFull = FftSharp.Transform.FFTfreq(sampleRate, signal, oneSided: false); + double[] freqsFull2 = FftSharp.Transform.FFTfreq(sampleRate, signal.Length, oneSided: false); + Assert.That(freqsFull, Is.EqualTo(freqsFullKnown)); + Assert.That(freqsFull2, Is.EqualTo(freqsFullKnown)); + + double[] freqsHalfKnown = { + 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000 + }; + double[] freqsHalf = FftSharp.Transform.FFTfreq(sampleRate, fft, oneSided: true); + double[] freqsHalf2 = FftSharp.Transform.FFTfreq(sampleRate, fft.Length, oneSided: true); + Assert.That(freqsHalf, Is.EqualTo(freqsHalfKnown)); + Assert.That(freqsHalf2, Is.EqualTo(freqsHalfKnown)); + } + } +} diff --git a/src/FftSharp.Tests/FreqLookup.cs b/src/FftSharp.Tests/FreqLookup.cs deleted file mode 100644 index d47a510..0000000 --- a/src/FftSharp.Tests/FreqLookup.cs +++ /dev/null @@ -1,28 +0,0 @@ -using NUnit.Framework; -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; - -namespace FftSharp.Tests -{ - class FreqLookup - { - [Test] - public void Test_Freq_Lookup() - { - int sampleRate = 8000; - int pointCount = 10; - - // look up frequencies corresponding to FFT points - double[] freqs = FftSharp.Transform.FFTfreq(sampleRate, pointCount); - - // frequencies should span from 0 to nyquist frequency (half the sample rate) - Assert.AreEqual(new double[] { 0, 400, 800, 1200, 1600, 2000, 2400, 2800, 3200, 3600 }, freqs); - Assert.AreEqual(sampleRate / 2, freqs.Last() + freqs[1]); - - // frequencies should be spaced by the frequency period - Assert.AreEqual(400, FftSharp.Transform.FFTfreqPeriod(sampleRate, pointCount)); - } - } -} diff --git a/src/FftSharp.Tests/Transform.cs b/src/FftSharp.Tests/Transform.cs index bf2d137..7507c5f 100644 --- a/src/FftSharp.Tests/Transform.cs +++ b/src/FftSharp.Tests/Transform.cs @@ -155,5 +155,68 @@ public void Test_FftInput_ThrowsIfLengthIsNotPowerOfTwo(int length, bool shouldT Assert.DoesNotThrow(realSpanRFFT); } } + + [TestCase(0, true)] + [TestCase(1, false)] + [TestCase(2, false)] + [TestCase(4, false)] + [TestCase(8, false)] + [TestCase(16, false)] + [TestCase(32, false)] + [TestCase(64, false)] + public void Test_Fft_Complex_DifferentLengths(int pointCount, bool shouldFail) + { + Complex[] signal = new Complex[pointCount]; + if (shouldFail) + { + Assert.Throws(() => FftSharp.Transform.FFT(signal)); + } + else + { + Assert.DoesNotThrow(() => FftSharp.Transform.FFT(signal)); + } + } + + [TestCase(0, true)] + [TestCase(1, false)] + [TestCase(2, false)] + [TestCase(4, false)] + [TestCase(8, false)] + [TestCase(16, false)] + [TestCase(32, false)] + [TestCase(64, false)] + public void Test_Fft_Double_DifferentLengths(int pointCount, bool shouldFail) + { + double[] signal = new double[pointCount]; + if (shouldFail) + { + Assert.Throws(() => FftSharp.Transform.FFT(signal)); + } + else + { + Assert.DoesNotThrow(() => FftSharp.Transform.FFT(signal)); + } + } + + [TestCase(0, true)] + [TestCase(1, true)] + [TestCase(2, true)] + [TestCase(4, true)] + [TestCase(8, true)] + [TestCase(16, false)] + [TestCase(32, false)] + [TestCase(64, false)] + public void Test_Fft_Magnitude_DifferentLengths(int pointCount, bool shouldFail) + { + double[] signal = new double[pointCount]; + if (shouldFail) + { + Assert.Throws(() => FftSharp.Transform.FFTmagnitude(signal)); + } + else + { + Assert.DoesNotThrow(() => FftSharp.Transform.FFTmagnitude(signal)); + } + } } } \ No newline at end of file diff --git a/src/FftSharp/FftSharp.csproj b/src/FftSharp/FftSharp.csproj index 90eda11..c0f92d3 100644 --- a/src/FftSharp/FftSharp.csproj +++ b/src/FftSharp/FftSharp.csproj @@ -16,9 +16,9 @@ https://github.com/swharden/FftSharp/releases true FftSharp.xml - 1.1.5 - 1.1.5.0 - 1.1.5.0 + 1.1.6 + 1.1.6.0 + 1.1.6.0 1591 portable true @@ -34,7 +34,7 @@ - + runtime; build; native; contentfiles; analyzers; buildtransitive all @@ -42,7 +42,7 @@ - + diff --git a/src/FftSharp/Transform.cs b/src/FftSharp/Transform.cs index 725179a..5f07481 100644 --- a/src/FftSharp/Transform.cs +++ b/src/FftSharp/Transform.cs @@ -122,13 +122,16 @@ private static int BitReverse(int value, int maxValue) /// /// Calculate sample frequency for each point in a FFT /// + /// Sample rate (Hz) of the original signal + /// Number of points to generate (typically the length of the FFT) + /// Whether or not frequencies are for a one-sided FFT (containing only real numbers) public static double[] FFTfreq(double sampleRate, int pointCount, bool oneSided = true) { double[] freqs = new double[pointCount]; if (oneSided) { - double fftPeriodHz = sampleRate / pointCount / 2; + double fftPeriodHz = sampleRate / (pointCount - 1) / 2; // freqs start at 0 and approach maxFreq for (int i = 0; i < pointCount; i++) @@ -151,6 +154,17 @@ public static double[] FFTfreq(double sampleRate, int pointCount, bool oneSided } } + /// + /// Calculate sample frequency for each point in a FFT + /// + /// Sample rate (Hz) of the original signal + /// FFT array for which frequencies should be generated + /// Whether or not frequencies are for a one-sided FFT (containing only real numbers) + public static double[] FFTfreq(double sampleRate, double[] fft, bool oneSided = true) + { + return FFTfreq(sampleRate, fft.Length, oneSided); + } + /// /// Return the distance between each FFT point in frequency units (Hz) /// @@ -294,6 +308,9 @@ public static double[] FFTmagnitude(double[] input) /// real input public static void FFTmagnitude(Span destination, Span input) { + if (input.Length < 16) + throw new ArgumentException("This overload requires an input with at least 16 points"); + if (!IsPowerOfTwo(input.Length)) throw new ArgumentException("Input length must be an even power of 2");