Talk:Kaiser window

DTFT of Kaiser window

The DTFT of Kaiser window is probably incorrect.

I compared in MatLab the plot of the formula of W_c to the plot of fftshift(fft(w_c)) and they look very different...

where do you find this formula??

fft computes the DFT, not the DTFT. If you want to approximate the DTFT you will need to zero-pad the input a lot. —Steven G. Johnson (talk) 15:21, 16 October 2008 (UTC)

yes, also making fftshift(fft(w_c),N),that is the N-point FFT, padded with zeros if w_c has less than N points, the plot is very very different.

try it... —Preceding unsigned comment added by 62.101.65.6 (talk) 17:18, 16 October 2008 (UTC)

if you want this is a link to a m-file that show this problem:

—Preceding unsigned comment added by 62.101.65.6 (talk) 17:55, 16 October 2008 (UTC)

Figure 2 and Figure 3 look a lot more alike if you remove the real() operator in the $g_{k}(l,:)\,$ formula and plot the abs() of $g_{k}(l,:)\,$ and $g_{f}(l,:).\,$ The horizontal scales appear to be the same, but the vertical scales appear to be slightly different. Also note that $W_{K}(\omega )\,$ is real-valued, which requires its inverse transform to be symmetrical about the Y-axis.

So apparently most, if not all, of the problem is that the $W_{K}(\omega )\,$ formula shown is actually the DTFT of $\omega _{0}(n),\,$ as defined at Window_function#Window_examples.

--Bob K (talk) 03:00, 17 October 2008 (UTC)

The (N+1) factor is also wrong (should be N). I fixed both errors, and now the formula gives the same real and imaginary results as the zero-padded DFT.

--Bob K (talk) 03:48, 17 October 2008 (UTC)

ok!! thanks a lot —Preceding unsigned comment added by 62.101.65.6 (talk) 09:15, 17 October 2008 (UTC)

Oops. After some sleep, I now notice that changing the N+1 factor to N appears to fix only the green and yellow curves (a=4.5787 & 6.8955). But now the red and blue curves (a= 0 & 1.8950) don't match. Close but no cigar, so I'm going to change N back to N+1 in the article until we know what's right.

--Bob K (talk) 11:52, 17 October 2008 (UTC)

yes there's still something a little wrong...

moreover i try to calculate the zero of the abs() of DTFT. the first should be at

$w={\frac {2\cdot \alpha }{N+1}}$,

when the argument of arcsin is 0 (also the denominator goes to 0, but it goes more "slowly"), but MatLab doesn't agree with me....

mystery of Kaiser —Preceding unsigned comment added by 62.101.65.6 (talk) 13:14, 17 October 2008 (UTC)

Got it:

${\frac {(N)\cdot \sinh \left({\sqrt {\alpha ^{2}-\left({\frac {(N)\cdot \omega }{2}}\right)^{2}}}\right)}{I_{0}(\alpha )\cdot {\sqrt {\alpha ^{2}-\left({\frac {(N)\cdot \omega }{2}}\right)^{2}}}}}\ \cdot \ e^{-j\omega {\frac {N}{2}}}$

is the DTFT of the sequence:

${\frac {I_{0}\left(\alpha {\sqrt {1-\left({\frac {2n}{N}}-1\right)^{2}}}\right)}{I_{0}(\alpha )}}\quad \quad \underbrace {0\leq n<N} _{{\mbox{not}}\quad 0\leq n\leq N}$

Therefore:

${\frac {(N+1)\cdot \sinh \left({\sqrt {\alpha ^{2}-\left({\frac {(N+1)\cdot \omega }{2}}\right)^{2}}}\right)}{I_{0}(\alpha )\cdot {\sqrt {\alpha ^{2}-\left({\frac {(N+1)\cdot \omega }{2}}\right)^{2}}}}}\ \cdot \ e^{-j\omega {\frac {N+1}{2}}}$

is the DTFT of the sequence:

${\frac {I_{0}\left(\alpha {\sqrt {1-\left({\frac {2n}{N+1}}-1\right)^{2}}}\right)}{I_{0}(\alpha )}}\quad \quad 0\leq n\leq N,$

NOT:

${\frac {I_{0}\left(\alpha {\sqrt {1-\left({\frac {2n}{N}}-1\right)^{2}}}\right)}{I_{0}(\alpha )}}\quad \quad 0\leq n\leq N.$

--Bob K (talk) 14:59, 17 October 2008 (UTC)

mmm, i come to a different conclusion (unfortunately....)

the right formula for the old window ($0\leq n\leq N$ is:

$\underbrace {\frac {(N)\cdot \sinh \left({\sqrt {\alpha ^{2}-\left({\frac {N\cdot \omega }{2}}\right)^{2}}}\right)}{I_{0}(\alpha )\cdot {\sqrt {\alpha ^{2}-\left({\frac {N\cdot \omega }{2}}\right)^{2}}}}} _{W_{K}(\omega )}\ \cdot \ e^{-j\omega {\frac {N}{2}}}\quad \quad {\mbox{for }}-\pi \leq \omega \leq \pi ,$

but it's only an approximation; try this formula with N=50, the two DTFT are the same. \\ probably the real DTFT formula is incalculable (sure i'm not able to do...)

your first window ($0\leq n<N$) isn't symmetric for N even... —Preceding unsigned comment added by Belox (talk • contribs) 15:24, 17 October 2008 (UTC)

We agree on everything except the inclusion of w(N) or not. Using your .m file (appropriately modified), and comparing the abs() values, I get the best agreement with w(N) excluded. Why did you suggest trying N=50? The inclusion or exclusion of one point is much less significant, and it just makes the differences much harder to see. But I reached the same conclusion as the N=12 and N=11 cases.

--Bob K (talk) 12:10, 18 October 2008 (UTC)

Based on my tests so far, $W_{K}(\omega )\,$ is not just an approximation. It's right on the money, provided you exclude w(N) from the window function. I'm talking about its magnitude, not the unimportant linear phase factor. The symmetry issue for N-even is also unimportant. It only matters for the academic exercise of trying to create a function centered at n=0 whose DTFT is real-valued. But most actual windowing is done like Window_function#Window_examples,with an even number of points and a window centered at (N-1)/2; i.e. the visible peak is flat-topped at both w(N/2-1) and w(N/2+1). So the DTFT generally has linear phase and therefore is not real-valued.

--Bob K (talk) 12:42, 18 October 2008 (UTC)

To create the same type of even-numbered, flat-topped window from the w(n) function in this article, we need to shift the peak from N/2 to (N-1)/2; i.e. we need the function:

$w\left(n+{\frac {1}{2}}\right)={\frac {I_{0}\left(\alpha {\sqrt {1-\left({\frac {2(n+.5)}{N}}-1\right)^{2}}}\right)}{I_{0}(\alpha )}}\quad \quad 0\leq n\leq N-1,$

where N is even. And this function has the 2-sided symmetry that you were worried about.

--Bob K (talk) 14:12, 18 October 2008 (UTC)

I don't think that a fractional index (n+0.5) can be accepted in a discrete function.

I found this:

http://demonstrations.wolfram.com/KaiserWindowTransform/

it seem that the equation of DTFT is in fact the fourier transform of the continuous function $w_{k}(t)$ and not the true DTFT of $w_{k}(n)$ —Preceding unsigned comment added by Belox (talk • contribs) 20:23, 19 October 2008 (UTC)

Good point!

OK, if that means:

$\underbrace {\frac {I_{0}\left(\alpha {\sqrt {1-\left({\frac {2t}{N}}-1\right)^{2}}}\right)}{I_{0}(\alpha )}} _{w(t)}\quad {\stackrel {\mathcal {F}}{\Longleftrightarrow }}\quad {\frac {N\cdot \sinh \left({\sqrt {\alpha ^{2}-\left({\frac {N\cdot \omega }{2}}\right)^{2}}}\right)}{I_{0}(\alpha )\cdot {\sqrt {\alpha ^{2}-\left({\frac {N\cdot \omega }{2}}\right)^{2}}}}}\cdot e^{-j\omega {\frac {N}{2}}}$

Then:

${\frac {I_{0}\left(\alpha {\sqrt {1-\left({\frac {2n}{N}}-1\right)^{2}}}\right)}{I_{0}(\alpha )}}\quad {\stackrel {DTFT}{\Longleftrightarrow }}\quad \sum _{k=-\infty }^{\infty }{\frac {N\cdot \sinh \left({\sqrt {\alpha ^{2}-\left({\frac {N\cdot (\omega -k*2\pi )}{2}}\right)^{2}}}\right)}{I_{0}(\alpha )\cdot {\sqrt {\alpha ^{2}-\left({\frac {N\cdot (\omega -k*2\pi )}{2}}\right)^{2}}}}}\cdot e^{-j\omega {\frac {N}{2}}}$

And you are right... the k=0 term (which is what we have been plotting) is just an approximation.

But also note that we are not restricted to sampling w(t) at t = 0, 1, 2, ...

We can also sample it at t = .5, 1.5, 2.5, ...

That is why I can write w(n+.5). No problem.

--Bob K (talk) 22:36, 19 October 2008 (UTC)

And still, we're talking about an approximation. To force the $w_{n}$ sequence to be of length N, we would have to define w(t) as:

${\frac {I_{0}\left(\alpha {\sqrt {1-\left({\frac {2t}{N}}-1\right)^{2}}}\right)}{I_{0}(\alpha )}}\cdot \operatorname {rect} \left({\frac {t-{\frac {N}{2}}}{N}}\right),$

And its transform would be:

$\left({\frac {N\cdot \sinh \left({\sqrt {\alpha ^{2}-\left({\frac {N\cdot \omega }{2}}\right)^{2}}}\right)}{I_{0}(\alpha )\cdot {\sqrt {\alpha ^{2}-\left({\frac {N\cdot \omega }{2}}\right)^{2}}}}}*N\cdot \operatorname {sinc} \left({\frac {N\omega }{2\pi }}\right)\right)\cdot e^{-j\omega {\frac {N}{2}}}$

--Bob K (talk) 23:22, 19 October 2008 (UTC)

ok, but i think that also the phase should be periodic. —Preceding unsigned comment added by Belox (talk • contribs) 20:10, 20 October 2008 (UTC)

When N/2 is an integer, $e^{-j\omega {\frac {N}{2}}}=e^{-j(\omega +k\cdot 2\pi ){\frac {N}{2}}}$. Also, see the "Shift in time" row of the table at DTFT#Properties.

--Bob K (talk) 20:59, 20 October 2008 (UTC)

ok, you done a great work! —Preceding unsigned comment added by Belox (talk • contribs) 23:13, 27 October 2008 (UTC)

I am added word "approximate" to Fourier transform there, because it mentioned as approximate in original paper cited there.

Interestingly, exact Fourier transform of that approximation is by itself excellent window function, much more attractive than original KB window, because it have exact closed form in frequency domain. It is because that new window have intrinsically bound support (without need of artificial cutting of it by rect() in time domain), hence inverse transform reconstructs it exactly into frequency domain. — Preceding unsigned comment added by Lexey73 (talk • contribs) 09:43, 23 December 2021 (UTC)

two different alphas

I just noticed the 27-Feb-2018 change described as "Replaced function definitions with actual definitions from a paper by Kaiser. Reverted images to ones that use alpha as defined by Kaiser." It replaces πα with just "α". Because of that change, we currently have two different alphas (in Window_function#Kaiser_window and Kaiser window), which is unnecessarily confusing. The α used in Window_function#Kaiser_window makes it very convenient to compare the DPSS windows (which are parameterized with α) to their Kaiser approximations, for equal values of α. That Kaiser formula comes from a legitimate source: Harris, Fredric J. (Jan 1978), p 73, which attributes it to two 1961 Bell System Journals by Slepian, Pollak, and Landau. I also followed this link: https://ccrma.stanford.edu/~jos/sasp/Kaiser_DPSS_Windows_Compared.html and found that the product πα is referred to as β (not just "α"). Both α and β are displayed on the figures (a totally avoidable inconvenience), but the text refers only to α (the same one as Window_function#Kaiser_window). For example: "The following Matlab comparison of the DPSS and Kaiser windows illustrates the interpretation of α as the bin number of the edge of the critically sampled window main lobe". Is there an equivalently useful interpretation of β (which is the "α" in Kaiser window)? And this link: https://ccrma.stanford.edu/~jos/sasp/Kaiser_Window.html uses β (not just "α") in the formula. It also concedes "Sometimes the Kaiser window is parametrized by α, where β=πα." And finally, https://www.mathworks.com/help/signal/ug/kaiser-window.html makes no mention of any α: "The statement kaiser(n,beta) computes a length n Kaiser window with parameter beta." Therefore, the preponderance of evidence (so far) is for either showing πα or β (not just "α"). My preference is πα, because that α is comparable to the DPSS α.
--Bob K (talk) 14:06, 22 March 2019 (UTC)

Kaiser's original paper, as published in Selected Papers in Digital Signal Processing II, IEEE Press, New York, 1975, uses just "α". --MindlessXD (talk) 22:58, 25 October 2019 (UTC)

Yes, that was explained with the 27-Feb-2018 edit. And Harris (1978) introduces the "πα" convention, where α has a physical significance (bin number of the edge of the window main lobe). Following the lead of Mathworks and ccrma.stanford.edu, we have chosen to avoid confusion by defining Harris' πα (and Kaiser's α) as β.

--Bob K (talk) 19:14, 29 October 2019 (UTC)

Window Formula

I noticed that windows made with the previously in this article shown formula

${\frac {I_{0}\left(\pi \alpha {\sqrt {1-\left({\frac {2k}{N}}-1\right)^{2}}}\right)}{I_{0}(\pi \alpha )}}$

were different from the ones created in Matlab using the function kaiser. Thus I checked the references of J.F. Kaiser and found that he defined the window formula without the $\pi$.

${\frac {I_{0}\left(\alpha {\sqrt {1-\left({\frac {2k}{N}}-1\right)^{2}}}\right)}{I_{0}(\alpha )}}$

I corrected the formula and the two figures in the article so that it fits the definitions from Kaiser. In fact it does not change the shape of the window, but it requires now by 3.14 larger $\alpha$ constants to get the same window as before. --RetoGalli 20:09, 19 September 2007 (UTC)

Hi Reto, I was a bit confused about the ${\alpha }$, because I could exactly reproduce the upper plot with matlab, but could not reproduce the lower plot until I read your note, and multiplied the ${\alpha }$ by ${\pi }$. Maybe it is worth putting a note in the text. (Heiko, 27.Feb 2008, 16:15). —Preceding unsigned comment added by 141.83.61.65 (talk) 15:17, 27 February 2008 (UTC)

Heiko, I do not exactly understand your comment. If I calculate a window with the lower function as it is in the article right now this looks like this in matlab:

N = 64; alpha = 2; n = 0:N; wk1 = besseli(0,alpha*sqrt(1-(2*n./N-1).^2))./besseli(0,alpha);

The result of this calculation is identical with the result of the two matlab functions window or kaiser:

wk2 = kaiser(N+1,alpha); wk3 = window(@kaiser,N+1,alpha);

The equation as it is presented in the article gives the same result as the matlab functions and is consitent with the paper from J.F.Kaiser. I do not see what you would like to have in an additional note in the text. But if you see a way to make the whole thing more clear, feel free to add a note. --RetoGalli (talk) 20:32, 28 February 2008 (UTC)

Sorry, my comment was unclear. What I meant is this: As you said, the Kaiser windows in the upper plot can be reproduced with matlab wk=kaiser(128,alpha); with alpha=1,2,4,8 and 16. But I needed Kaiser-Bessel derived windows (KBD), and for this I encoded the KBD formula given in the lower wiki section:

 N=100;kw=kaiser(N,alpha);fac=sum(kw(1:N));kw_cum=cumsum(kw(1:N));

and

 sqr_sum=sqrt( kw_cum./fac );kbd_win=[sqr_sum ; sqr_sum((N:-1:1))];

where kbd_win should reproduce the windows in the KBD plot, if alpha=0.5,2.0,8.0, and 32. However, to reproduce them exactly, I needed to enter alpha_2=pi*alpha. So, if I'm not completly confused (which is absolutely possible), the ${\alpha }$'s used in the two plots differ by ${\pi }$. I guess, the lower plot has been done with the C++ Code from the Bosi-Site.(Heiko, 29.Feb 2008, 13:34). —Preceding unsigned comment added by 141.83.61.65 (talk) 13:00, 29 February 2008 (UTC)

Ok, I got it. When I corrected the window formula last year, I did not check the section about KBD. I corrected now also the KBD plot accordingly. I think it should be consistent now. --RetoGalli (talk) 21:29, 29 February 2008 (UTC)

I just got burned by this factor of π, since I designed my filter with MATLAB, but implemented the filter according to the definition here. So, feel free to edit or move what I wrote, but I really think that the factor of π difference should be prominently noted in the article since different texts use different definitions of the shape factor. Derekt75 (talk) 04:32, 12 January 2010 (UTC)

Window Center

There is a discrepancy between here and http://wiki.riteme.site/wiki/Window_function#Kaiser_windows on the expression of Kaiser window. I found this while implementing some C++ code. This page's expression centers the peak at 0 while the other page is at N/2. As a newbie on DSP I hesitate to make the change but I did find it confusing, especially since this page mentions that the window is centered at N/2. Rud 18:17, 28 September 2006 (UTC)

My own comparison of the two expressions concludes that the other one spans [0,N-1], whose center is (N-1)/2, which is a little different than your observation. And this one spans [0,N], whose center is N/2, not 0. The two expressions are not different enough to behave as you say, so I advise you to check your code again. --Bob K 19:14, 28 September 2006 (UTC)

Specifically, check for this mistake (note the missing "1" term):

$w_{k}={\frac {I_{0}\left(\pi \alpha {\sqrt {1-(2k/N)^{2}}}\right)}{I_{0}(\pi \alpha )}}$

--Bob K 19:23, 28 September 2006 (UTC)

Or more likely you put it in the wrong place, like this:

$w_{k}={\frac {I_{0}\left(\pi \alpha {\sqrt {1-({\frac {2k}{N-1}})^{2}}}\right)}{I_{0}(\pi \alpha )}}$

What you want is this:

$w_{k}={\frac {I_{0}\left(\pi \alpha {\sqrt {1-({\frac {2k}{N}}-1)^{2}}}\right)}{I_{0}(\pi \alpha )}}$

--Bob K 19:29, 28 September 2006 (UTC)

Nearly Optimal

In what sense is the Kaiser window a "nearly optimal window function"? What does that mean? --Bob K 16:18, 12 December 2005 (UTC)

From Oppenheim, Schafer, and Buck:

The trade-off between the main-lobe width and side-lobe area can be quantified by seeking the window function that is maximally concentrated around ω=0 in the frequency domain. The issue was considered in depth in a series of classic papers by Slepian et al. (1961). The solution found in this work involves prolate spheroidal wave functions, which are difficult to compute and therefore unattractive for filter design. However, Kaiser (1966, 1974) found that a near-optimal window could be formed using the zeroth-order modified Bessel function of the first kind, a function that is much easier to compute.

Unfortunately, Oppenheim et al. is not very precise in defining what is meant by "maximally concentrated"; presumably it is in the references. Anyway, this is the sense in which that book calls it "near-optimal". If you have a chance to look up a more detailed explanation, that would be a nice addition to the article. —Steven G. Johnson 20:09, 12 December 2005 (UTC)

Frequency spectra

I replaced the frequency spectra plots by Bob K (old ones were Image:Kaiser_window_function_alpha_2.png and Image:Kaiser_window_function_alpha_3.png). I found the old ones hard to read because of the plotting format that was used — all sorts of Moire effects, etc. Also, it is nice to keep the same format as the time-domain plots we already have, in particular putting two α values on the same plot so that you can directly compare their spectra, and using a similar color and labelling scheme. (On the other hand, Bob's side-by-side plots look very nice in the context of the Window function article.) —Steven G. Johnson 22:19, 15 December 2005 (UTC)

Yours is definitely a better treatment of the "grass". All I can say in my defense is that in all my years of using window functions, no-one ever cared what the grass looks like below the peaks. And your overlay is essential to compare the lobe widths. I might have picked bolder colors, but good job! --Bob K 23:31, 15 December 2005 (UTC)

"Off-by-one error"

Not surprisingly, the strength of Wikipedia is also its glaring weakness... it takes constant vigilance to prevent erosion. Hopefully someday there will be a better way, but I don't know what it is. Thoughts? --Bob K 11:21, 7 May 2006 (UTC)

Wikipedia:Stable versions is one proposal to address this. Time will tell whether this works in practice. —Steven G. Johnson 17:49, 7 May 2006 (UTC)

Thanks. It's a start, at least. One thing I do not understand from that proposal is how one will know that a stable version exists? --Bob K 13:57, 8 May 2006 (UTC)

Kaiser Burnout

In the first line of the article, there is a reference to James Kaiser at Bell Laboratories. When I follow the link to James_Kaiser, I am redirected to the Kaiser_Burnout, an event in the american civil war. I don't think this is the correct James Kaiser reference. The link should be deleted. Simplimus (talk) 07:28, 27 November 2013 (UTC)

External links modified

Hello fellow Wikipedians,

I have just modified one external link on Kaiser window. Please take a moment to review my edit. If you have any questions, or need the bot to ignore the links, or the page altogether, please visit this simple FaQ for additional information. I made the following changes:

• Added archive https://web.archive.org/web/20090327144432/http://ccrma.stanford.edu/courses/422/projects/kbd/ to http://ccrma.stanford.edu/courses/422/projects/kbd/

When you have finished reviewing my changes, you may follow the instructions on the template below to fix any issues with the URLs.

This message was posted before February 2018. After February 2018, "External links modified" talk page sections are no longer generated or monitored by InternetArchiveBot. No special action is required regarding these talk page notices, other than regular verification using the archive tool instructions below. Editors have permission to delete these "External links modified" talk page sections if they want to de-clutter talk pages, but see the RfC before doing mass systematic removals. This message is updated dynamically through the template {{source check}} (last update: 5 June 2024).

• If you have discovered URLs which were erroneously considered dead by the bot, you can report them with this tool.
• If you found an error with any archives or the URLs themselves, you can fix them with this tool.

Cheers.—InternetArchiveBot (Report bug) 05:10, 5 December 2017 (UTC)

Apparently erroneous sign of Fourier transform

In doing some work based on the definition given for the FT of the window, I unified the two domains (|omega| <= omega_alpha and |omega| >= omega_alpha) into one expression: (2T/I_0(alpha))*sinc(alpha*i*sqrt(1-(omega/omega_alpha)^2)) .

Then I realised that I'd accidentally discarded the negative sign in the denominator of the expression for the |omega| <= omega_alpha domain. Although this is never explained in the article, omega_alpha appears to be the joining point where the main lobe meets the "ripples", and if you discard the 2T/I_0(alpha) coefficient, both expressions evaluate to unity at that point.

I therefore believe that this negative sign is wrong for these reasons:

• since the window is always nonnegative, its DC gain should be nonnegative
• the negative sign creates a discontinuity at omega_alpha rather than a smooth join; frequency-domain discontinuities imply infinite time-domain support, which is clearly inappropriate for a window function

Plotting the FT here: http://www.mathsisfun.com/data/function-grapher.php?func1=sinh(4*a*sqrt(1-x*x))%2F(4*a*sqrt(1-x*x))&func2=sin(4*a*sqrt(x*x-1))%2F(4*a*sqrt(x*x-1))&xmin=-6.950&xmax=6.950&ymin=-4.296&ymax=4.296 ...seems to confirm my thinking.

81.187.94.66 (talk) 10:44, 1 May 2018 (UTC)

Addendum: on looking at the reference for this formula, it seems that it was taken verbatim from James Kaiser's original paper. Although it seems unlikely that the paper is wrong, I can't spot the error in my reasoning. 81.187.94.66 (talk) 11:48, 1 May 2018 (UTC)

Perhaps I am misreading the equation in Kaiser's original paper? The horizontal line does seem a bit long to be a negative sign, but I'm not sure what else to make of it. I admit that other sources (example) for the equation online do not have the negative sign. --MindlessXD (talk) 04:19, 28 May 2018 (UTC)

Calculates KBD Windows of AAC

#include <stdio.h>
#include <math.h>

constexpr int fact(int N)
{
    return N ? N*fact(N-1) : 1;
}

constexpr double power(double x, int N)
{
    return N ? x*power(x, N-1) : 1.0;
}

constexpr double Pi = 3.1415926535897932384626433832795;

template <int n = 0>
constexpr double I(double x)
{
    double x2  = 0.5 * x;
    double x22 = x2 * x2;
    double S   = power (x2, n) / fact(n);
    double Sum = S;

    for (int i = 1; S; i++)
        Sum += (S *= x22 / i / (i + n));

    return Sum;
}

constexpr double KB (double alpha, double x)
{
    return I<0>(Pi*alpha*sqrt(1-power(x,2))) / I<0>(Pi*alpha);
}

constexpr double KB (double alpha, int n, double N)
{
    return KB(alpha, 2*n/N - 1);
}

double Sum(double x[], int N)
{
    double S = 0;
    for (int i = 0; i < N; i++)
        S += x[i];
    return S;
}

double KBD (double alpha, int n, int N)
{
    static double _Q [4096+1];
    static int    _N     = -1;
    static double _alpha = -1.0;

    if (alpha != _alpha || N != _N)
    {
        for (int i = 0; i <= N; i++)
            _Q[i] = KB (alpha, i, N);
        _alpha = alpha;
        _N     = N;
    }
    return sqrt (Sum(_Q,n+1) / Sum(_Q,N+1));
}

int main()
{
    for (int i = 0; i < 1024; i++)
        printf ("%4d %18.12f\n", i, KBD(4.0, i, 1024));
    printf ("\n");
    for (int i = 0; i < 128; i++)
        printf ("%4d %18.12f\n", i, KBD(6.0, i, 128));
    printf ("\n");
    return 0;
}

--2003:c3:671c:ac00:9951:83f4:72e:d01a

Excerpt from https://commons.wikimedia.org/wiki/File:Kbd-window.svg:

for k= 1:length(betas)
    w = besseli(0,betas(k)*sqrt(1-(2*(0:N)/N -1).^2))/besseli(0,betas(k));

    d = zeros(1,2*N);	% length 2N
    for n=0:N-1
    d(1+n) = sum(w(1+(0:n)));
    end

    for n=N:2*N-1
    d(1+n) = sum(w(1+(0:2*N-1-n)));
    end

    d = sqrt(d/sum(w(1+(0:N))));
    plot(0:2*N-1, d, "color", colors(k,:), "linewidth", 2)
end

--Bob K (talk) 14:20, 22 August 2021 (UTC)