Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Kaiser windowing in FftSharp #35

Closed
ADD-eNavarro opened this issue Oct 6, 2021 · 10 comments
Closed

Kaiser windowing in FftSharp #35

ADD-eNavarro opened this issue Oct 6, 2021 · 10 comments

Comments

@ADD-eNavarro
Copy link

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?

@swharden
Copy link
Owner

swharden commented Oct 6, 2021

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 👍

@swharden
Copy link
Owner

swharden commented Oct 6, 2021

For example the Hanning window is just:

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;
}

@ADD-eNavarro
Copy link
Author

Absolutelly!

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

@swharden
Copy link
Owner

swharden commented Oct 6, 2021

@swharden
Copy link
Owner

swharden commented Oct 7, 2021

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

@swharden swharden closed this as completed Oct 7, 2021
@swharden
Copy link
Owner

swharden commented Oct 7, 2021

... 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:

@swharden swharden reopened this Oct 7, 2021
@swharden
Copy link
Owner

swharden commented Oct 7, 2021

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

@swharden
Copy link
Owner

swharden commented Oct 7, 2021

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);

windows

@ADD-eNavarro
Copy link
Author

ADD-eNavarro commented Oct 7, 2021

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.

@swharden
Copy link
Owner

swharden commented Oct 8, 2021

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");

image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants