Window: Add Kaiser window

resolves #35

dev/python/bessel.py

@@ -0,0 +1,10 @@
+import numpy as np
+xs = np.arange(20)
+print(xs)
+print(np.i0(xs))
+
+import matplotlib.pyplot as plt
+window = np.kaiser(50, 14)
+print(window)
+plt.plot(window)
+plt.show()

src/FftSharp.Tests/Window.cs

@@ -22,6 +22,7 @@ public void Test_Window_Functions()
             //plt.PlotScatter(xs, FftSharp.Window.BlackmanHarris(xs.Length), label: "BlackmanHarris");
             plt.PlotScatter(xs, FftSharp.Window.FlatTop(xs.Length), label: "FlatTop");
             plt.PlotScatter(xs, FftSharp.Window.Cosine(xs.Length), label: "Cosine");
+            plt.PlotScatter(xs, FftSharp.Window.Kaiser(xs.Length, 15), label: "Kaiser");
 
             // customize line styles post-hoc
             foreach (var p in plt.GetPlottables())
@@ -51,5 +52,69 @@ public void Test_Window_Reflection()
                 Assert.AreEqual(5, windowed.Length);
             }
         }
+
+        [TestCase(0, 1)]
+        [TestCase(1, 1)]
+        [TestCase(2, 2)]
+        [TestCase(3, 6)]
+        [TestCase(9, 362880)]
+        public void Test_Factorial_MatchesKnown(int k, int factorial)
+        {
+            Assert.AreEqual(factorial, FftSharp.Window.Factorial(k));
+        }
+
+        [Test]
+        public void Test_Bessel_MatchesPython()
+        {
+            /* expected values calculated with python:
+                >>> import numpy as np
+                >>> np.i0(np.arange(20))
+            */
+
+            double[] expected = {
+                1.00000000e+00, 1.26606588e+00, 2.27958530e+00, 4.88079259e+00,
+                1.13019220e+01, 2.72398718e+01, 6.72344070e+01, 1.68593909e+02,
+                4.27564116e+02, 1.09358835e+03, 2.81571663e+03, 7.28848934e+03,
+                1.89489253e+04, 4.94444896e+04, 1.29418563e+05, 3.39649373e+05,
+                8.93446228e+05, 2.35497022e+06, 6.21841242e+06, 1.64461904e+07,
+             };
+
+            double[] actual = FftSharp.Window.BesselZero(expected.Length);
+
+            for (int i = 0; i < expected.Length; i++)
+            {
+                double allowableError = .00001 * expected[i];
+                Assert.AreEqual(expected[i], actual[i], allowableError);
+            }
+        }
+
+        [Test]
+        public void Test_Kaiser_MatchesPython()
+        {
+            /* expected values calculated with python:
+                >>> import numpy as np
+                >>> np.kaiser(50, 14)
+            */
+            double[] expected = {
+                7.72686684e-06, 8.15094846e-05, 3.26000767e-04, 9.42588751e-04, 2.26624847e-03, 
+                4.80567914e-03, 9.27621459e-03, 1.66164301e-02, 2.79789657e-02, 4.46873500e-02, 
+                6.81537432e-02, 9.97574012e-02, 1.40689805e-01, 1.91778970e-01, 2.53311387e-01, 
+                3.24874218e-01, 4.05241756e-01, 4.92328143e-01, 5.83222700e-01, 6.74315433e-01, 
+                7.61509399e-01, 8.40504954e-01, 9.07130390e-01, 9.57685605e-01, 9.89261639e-01, 
+                1.00000000e+00, 9.89261639e-01, 9.57685605e-01, 9.07130390e-01, 8.40504954e-01, 
+                7.61509399e-01, 6.74315433e-01, 5.83222700e-01, 4.92328143e-01, 4.05241756e-01, 
+                3.24874218e-01, 2.53311387e-01, 1.91778970e-01, 1.40689805e-01, 9.97574012e-02, 
+                6.81537432e-02, 4.46873500e-02, 2.79789657e-02, 1.66164301e-02, 9.27621459e-03, 
+                4.80567914e-03, 2.26624847e-03, 9.42588751e-04, 3.26000767e-04, 8.15094846e-05, 
+            };
+
+            double[] actual = FftSharp.Window.Kaiser(51, 14);
+
+            for (int i = 0; i < expected.Length; i++)
+            {
+                double allowableError = .00001 * expected[i];
+                Assert.AreEqual(expected[i], actual[i], allowableError);
+            }
+        }
     }
 }

src/FftSharp/Window.cs

@@ -138,6 +138,64 @@ public static double[] Cosine(int pointCount)
             return window;
         }
 
+        public static int Factorial(int k)
+        {
+            int result = 1;
+
+            for (int i = 2; i <= k; i++)
+                result *= i;
+
+            return result;
+        }
+
+        public static double[] BesselZero(int pointCount)
+        {
+            double[] window = new double[pointCount];
+
+            for (int i = 0; i < pointCount; i++)
+                window[i] = I0(i);
+
+            return window;
+        }
+
+        public static double I0(double x)
+        {
+            // Derived from code workby oygx210/navguide:
+            // https://github.com/oygx210/navguide/blob/master/src/common/bessel.c
+
+            double ax = Math.Abs(x);
+            if (ax < 3.75)
+            {
+                double y = Math.Pow(x / 3.75, 2);
+                double[] m = { 3.5156229, 3.0899424, 1.2067492, 0.2659732, 0.360768e-1, 0.45813e-2 };
+                return 1.0 + y * (m[0] + y * (m[1] + y * (m[2] + y * (m[3] + y * (m[4] + y * m[5])))));
+            }
+            else
+            {
+                double y = 3.75 / ax;
+                double[] m = { 0.39894228, 0.1328592e-1, 0.225319e-2, -0.157565e-2, 0.916281e-2, -0.2057706e-1, 0.2635537e-1, -0.1647633e-1, 0.392377e-2 };
+                return (Math.Exp(ax) / Math.Sqrt(ax)) * (m[0] + y * (m[1] + y * (m[2] + y * (m[3] + y * (m[4] + y * (m[5] + y * (m[6] + y * (m[7] + y * m[8]))))))));
+            }
+        }
+
+        public static double[] Kaiser(int pointCount, double beta)
+        {
+            // derived from python/numpy:
+            // https://github.com/numpy/numpy/blob/v1.21.0/numpy/lib/function_base.py#L3267-L3392
+
+            int M = pointCount;
+            double alpha = (M - 1) / 2.0;
+
+            double[] window = new double[pointCount];
+
+            for (int n = 0; n < pointCount; n++)
+            {
+                window[n] = I0(beta * Math.Sqrt(1 - (Math.Pow((n - alpha) / alpha, 2)))) / I0(beta);
+            }
+
+            return window;
+        }
+
         public static string[] GetWindowNames()
         {
             return typeof(Window)