Kaiser windowing in FftSharp

Hi!

Great work here! I am looking forward to use your FftSharp library, but I would need to have a Kaiser windowing implementation working. Do you plan on implementing that? Or else, do you know some other library I could use tith Kaiser windowing implemented?

Hi @ADD-eNavarro

Thanks for suggesting this! I can add this feature this evening.

The quickstart on the main readme shows how to use a window. A window is just a double[] so it should be easy to use a Kaiser or other custom window with the library how it is now 👍

For example the Hanning window is just:

FftSharp/src/FftSharp/Window.cs

Lines 67 to 75 in 6f51c47

    
           public static double[] Hamming(int pointCount) 
        
           { 
        
               double[] window = new double[pointCount]; 
        
               for (int i = 0; i < pointCount; i++) 
        
                   window[i] = 0.54 - 0.46 * Math.Cos(2 * Math.PI * i / pointCount); 
        
               return window; 
        
           }

Absolutelly!

Here's an implementation in Python, may be helpful if you're going to implement Kaiser window in your awesome tool.

Note to self: https://numpy.org/doc/stable/reference/generated/numpy.kaiser.html

So this turned out to be way harder than I expected...

Calculating the Kaiser window is easy. It's a modified zero-order Bessel function. In python it's just this (source):

w = np.i0(beta * np.sqrt(1 - (2 * n / (N - 1) - 1) ** 2)) / np.i0(beta)

The hard part is figuring out what that np.i0() is doing. It looks like it's generating a "Modified Bessel function of the first kind, order 0." https://numpy.org/doc/stable/reference/generated/numpy.i0.html

Snaking through the numpy source code I land at a C call and am struggling to drill down to the actual source code of that function.

There are some citations to math books published in the 60s but I don't have access to those.

One trail led to https://metacpan.org/dist/Math-Cephes/view/lib/Math/Cephes.pod#y0:-Bessel-function-of-the-second-kind,-order-zero which says this... basically it depends on j0(), another function we have to hunt down

DESCRIPTION:
Returns Bessel function of the second kind, of order
zero, of the argument.
 
The domain is divided into the intervals [0, 5] and
(5, infinity). In the first interval a rational approximation
R(x) is employed to compute
  y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
Thus a call to j0() is required.
 
In the second interval, the Hankel asymptotic expansion
is employed with two rational functions of degree 6/6
and 7/7.

I'm probably going to put this down and close this issue because this is taking a lot of my time and I want to devote my firepower elsewhere. I appreciate your suggestion @ADD-eNavarro, and would love to get a Bessel filter added to FftSharp (and then it will be easy to add a Kaiser window). If you can come-up with a function that returns Bessel values with C#, I'd be happy to pick this up again.

Based on the Numpy docs I'd expect it to look like this, so if you get it to work we know we did it right 👍

double[] BesselFilter(int pointCount)
{
    /* magic */
};

Console.WriteLine(string.join(", ", BesselFilter(4)));

1, 1.26606588, 2.2795853 , 4.88079259

... okay I kept looking and got a lead, so I'm opening this issue back up. Accord.NET is GPL-licensed so we have to be very careful not to let that code get into the FftSharp library or it would lose its ability to be released under the MIT license...

But here's a zero-order Bessel function of the first kind which according to How to Create a Configurable Filter Using a Kaiser Window is what we need:

Accord.Math/Functions/Bessel.cs#L376-L406

I found a repository that has a C++ Bessel function generator that is MIT licensed. Thanks @oygx210!

What a roller coaster... but it's in now and published on NuGet (version 1.0.12)

Let me know if it doesn't do the trick! Good luck with your project 👍 🚀

https://github.com/swharden/FftSharp#windowing

double[] window = FftSharp.Window.Kaiser(audio.Length, beta: 15);

Wow, such a commitment!
Thank you very much for your time, sorry I wasn't responsive, I guess we don't live in the same part of the planet, since you've worked through my night. I'll be testing this right away, be back in a little while with results :^)

Lately: Actully, I'm using Kaiser window to resample, but I'm having trouble adapting to C# the code from resampy, so still no real results.
I'll be posting when I have some.

Following-up, I found this diagram helpful

Plot plt = new Plot(500, 400);
plt.Palette = Palette.Tsitsulin;

double[] xs = ScottPlot.DataGen.Range(-1, 1, .01, true);
for (int i = 0; i < 15; i++)
{
    int beta = 3 + i * i;
    var sp = plt.AddScatter(xs, FftSharp.Window.Kaiser(xs.Length, beta), label: $"beta {beta}");
    sp.MarkerSize = 0;
    sp.LineWidth = 3;
    sp.Color = System.Drawing.Color.FromArgb(200, sp.Color);
}

plt.Title("Kaiser Window");
plt.Legend(enable: true, location: Alignment.UpperRight);
plt.SaveFig("kaiser.png");

	public static double[] Hamming(int pointCount)
	{
	double[] window = new double[pointCount];

	for (int i = 0; i < pointCount; i++)
	window[i] = 0.54 - 0.46 * Math.Cos(2 * Math.PI * i / pointCount);

	return window;
	}