Merge pull request #50 from arthurits/main

FFTfreq: fix off-by-one error

Author: swharden

dev/quickstart/fft-windowed.png

@@ -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));
+        }
+    }
+}

dev/quickstart/fft.png

@@ -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<ArgumentException>(() => 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<ArgumentException>(() => 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<ArgumentException>(() => FftSharp.Transform.FFTmagnitude(signal));
+            }
+            else
+            {
+                Assert.DoesNotThrow(() => FftSharp.Transform.FFTmagnitude(signal));
+            }
+        }
     }
 }

dev/quickstart/periodogram.png

@@ -16,9 +16,9 @@
     <PackageReleaseNotes>https://github.com/swharden/FftSharp/releases</PackageReleaseNotes>
     <IncludeSymbols>true</IncludeSymbols>
     <DocumentationFile>FftSharp.xml</DocumentationFile>
-    <Version>1.1.5</Version>
-    <AssemblyVersion>1.1.5.0</AssemblyVersion>
-    <FileVersion>1.1.5.0</FileVersion>
+    <Version>1.1.6</Version>
+    <AssemblyVersion>1.1.6.0</AssemblyVersion>
+    <FileVersion>1.1.6.0</FileVersion>
     <NoWarn>1591</NoWarn>
     <DebugType>portable</DebugType>
     <IncludeSymbols>true</IncludeSymbols>
@@ -34,15 +34,15 @@
   </ItemGroup>
 
   <ItemGroup>
-    <PackageReference Include="Microsoft.SourceLink.GitHub" Version="1.0.0">
+    <PackageReference Include="Microsoft.SourceLink.GitHub" Version="1.1.1">
       <IncludeAssets>runtime; build; native; contentfiles; analyzers; buildtransitive</IncludeAssets>
       <PrivateAssets>all</PrivateAssets>
     </PackageReference>
   </ItemGroup>
 
   <ItemGroup Condition="'$(TargetFramework)' == 'netstandard2.0'">
     <PackageReference Include="System.ValueTuple" Version="4.5.0" />
-    <PackageReference Include="System.Memory" Version="4.5.4" />
+    <PackageReference Include="System.Memory" Version="4.5.5" />
   </ItemGroup>
 
 </Project>

src/FftSharp.Tests/FftFreqTests.cs

@@ -122,13 +122,16 @@ private static int BitReverse(int value, int maxValue)
         /// <summary>
         /// Calculate sample frequency for each point in a FFT
         /// </summary>
+        /// <param name="sampleRate">Sample rate (Hz) of the original signal</param>
+        /// <param name="pointCount">Number of points to generate (typically the length of the FFT)</param>
+        /// <param name="oneSided">Whether or not frequencies are for a one-sided FFT (containing only real numbers)</param>
         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
             }
         }
 
+        /// <summary>
+        /// Calculate sample frequency for each point in a FFT
+        /// </summary>
+        /// <param name="sampleRate">Sample rate (Hz) of the original signal</param>
+        /// <param name="fft">FFT array for which frequencies should be generated</param>
+        /// <param name="oneSided">Whether or not frequencies are for a one-sided FFT (containing only real numbers)</param>
+        public static double[] FFTfreq(double sampleRate, double[] fft, bool oneSided = true)
+        {
+            return FFTfreq(sampleRate, fft.Length, oneSided);
+        }
+
         /// <summary>
         /// Return the distance between each FFT point in frequency units (Hz)
         /// </summary>
@@ -294,6 +308,9 @@ public static double[] FFTmagnitude(double[] input)
         /// <param name="input">real input</param>
         public static void FFTmagnitude(Span<double> destination, Span<double> 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");