The Kaiser Window

Music 420 / EE 265
	Winter 1992
	Julius O. Smith III
CCRMA, Stanford University

Definitions

Constants

M = 27;		(* window length *)
NN = 512;	(* DFT length *)
beta = 5;	(* Kaiser window time-bandwidth/2 parameter *)

nl = -(M-1)/2;	(* Lowest window time index *)
nu = (M-1)/2;	(* Highest window time index *)

I0 [x] - The 0th-Order Modified Bessel Function of the 1st Kind

I0[x_] := BesselI[0,x];

Plot[I0[x],{x,0,beta}];

Comparison of I0[x] to Exp[x]

?Series

Series[Exp[x],{x,0,4}]

Series[I0[x],{x,0,4}]

Note that the series for I0[x] is term-wise square of the series for Exp[x/2]:

Series[Exp[x/2],{x,0,4}]

Plot[{I0[x], Exp[x/2]},{x,0,beta}];

The FFT

?Fourier

Experimentation will show that Mathematica's Fourier function is smart about powers of 2.  If you give it a list of numbers which is a power of 2 in length, execution speed is much faster.  Here are two examples

Fourier[Range[4096]];	(* Takes about two seconds *)

Fourier[Range[4097]];	(* Takes about ten seconds *)

One thing we have to watch out for is alternate definitions of the Fourier transform.  Let's check it out:

Fourier[{1,0,0,0}]  (* See what it does to an impulse *)

Fourier[{0,1,0,0}]  (* Add a unit-sample delay *)

Chop[%]

The above indicate that Fourier normalizes by (1/Sqrt[NN]) (which we do not want) and it uses Exp[j 2 Pi k n / NN] as the kernel instead of Exp[-j 2 Pi k n / NN].  If we agree only to look at real sequences in the time domain, we may compensate as follows:

fft[x_] := Chop[N[Sqrt[Length[x]] Conjugate[Fourier[x]]]]

fft[{1,0,0,0}]  (* See what it does to an impulse *)

fft[{0,1,0,0}]  (* Look at Exp[-j 2 Pi k n / NN] *)

That's better

The Kaiser Window

kaiserUN[n_] := If[ n < nl || n > nu, 0,
		N[I0[beta Sqrt[1 - (n/nu)^2]]/I0[beta]]];

kscale = 1/N[Sum[kaiserUN[n],{n,nl,nu}]]

kaiser[n_] := kscale kaiserUN[n];

Display the Kasier window itself:

Plot[kaiser[n],{n,nl,nu}];

Display the Kasier window and the zero padding to be used in the DFT:

Plot[kaiser[n],{n,-NN/2,NN/2-1}];

The Kaiser Window Transform

Theoretical Approximation

The following closed-form expression for the transform appears in Harris.  It can be thought of as a "hyberbolic sinc function" (which I suppose should be denoted "sinch") when beta is zero.  The beta parameter serves to flatten the function near the origin.

Clear[Wh,WhUN];
WhUN[o_] := (M/I0[beta]) Sinh[z]/z /. z->Sqrt[beta^2-(M o/2)^2];
Wh[o_] := WhUN[o]/WhUN[0];

Wht = Table[N[Wh[2 Pi i / NN]],{i,-NN/2,NN/2 - 1}];

ListPlot[Wht, PlotRange->All];

Whtdb = 20 Log[10,Abs[Wht]];

pWhtdb = ListPlot[Whtdb];

Computed Window Transform

Note carefully the way the real, symmetric window must be arranged in the FFT input buffer in order to obtain a real transform:

wu = Table[kaiser[n],{n,0,nu}];
wl = Table[kaiser[n],{n,nl,-1}];
zp = 0 * Range[NN-M];
w = Flatten[{wu,zp,wl}];
ListPlot[w, PlotLabel->"Input to FFT"];

W = fft[w];	(* Kaiser window transform *)

Wdb = 20 Log[10,Abs[W]];

ListPlot[Wdb];

tmp1 = Table[Wdb[[i]],{i,1,NN/2}];
tmp2 = Table[Wdb[[i]],{i,NN/2+1,NN}];
Wdbro = Flatten[{tmp2,tmp1}];

pWdbro = ListPlot[Wdbro];

Compare theoretical transform to measured transform:

Show[pWhtdb, pWdbro];

Main Lobe Width based on Theoretical Approximation

numeratorW[o_] = Sinh[z] /. z->Sqrt[anyBeta^2-(anyM o/2)^2];
roots = Solve[numeratorW[omega]==0,omega];

roots[[1]] (* Only positive, real, root *)

Practical Design Formulas

Window Time-Bandwidth Parameter Beta

The input is the desired sidelobe attenuation in dB, e.g. 80

goodBeta[a_] := If[a<21,0,
		  If[a>50, 0.1102 (a-8.7),
		  	0.5842 (a-21)^0.4 + 0.07886 (a-21)]];

Window Transform Transition Width

deltaFN[a_] - The input is the desired sidelobe attenuation in dB, and the output is normalized as transition bandwidth in Hz times window length in seconds.  For the rectangular window (sinc window transform), this product is 1 if the transition bandwidth is taken as the distance from the peak of the main lobe to the first zero crossing of the sinc function.  Thus, the normalized transition width parameter will always be greater than 1, and the larger it is, the more sidelobe rejection we expect.

deltaFN[a_] := If[a<21,0.9222, (a-7.95)/14.36]

deltaBins[] - The input is the desired sidelobe attenuation in dB, and the output is in units of FFT bins at the current FFT size NN and window length M.

deltaBins[a_] := deltaFN[a] NN / M

High-Quality Audio Example

att2 = 100;	(* desired stop-band attenuation in dB *)

beta = goodBeta[att2]

The above value checks with Kaiser's paper

deltaFN[att2]

The above value also checks with Kaiser's paper

wu2 = Table[kaiser[n],{n,0,nu}];
wl2 = Table[kaiser[n],{n,nl,-1}];
zp = 0 * Range[NN-M];
w2 = Flatten[{wu2,zp,wl2}];
W2 = fft[w2];
Wdb2 = 20 Log[10,Abs[W2]];

pw2 = ListPlot[Flatten[{wl2,wu2}]];

pWdb2 = ListPlot[Wdb2];

Beta apears to be underestimated somewhat.  The observed sidelobe attenuation is less than 80 dB.

deltaBins[att2]

The transition width evidently accounts for the entire main lobe.

Evaluation and the Window Method for FIR Digital Filter Design

These results makes sense when you consider that Kaiser's application was the window method for FIR digital filter design (a topic we will take up if there is time).  In that application, the filter is obtained by multiplying the window function times the ideal impulse response (which is a sinc function in the case of an ideal lowpass filter).  That means the window transform gets convolved with the ideal frequency response.  As a result, the transition width of the lowpass filter is the whole main lobe width of the window.  (Think about the convolution of the window transform with a rectangular pulse in the frequency domain representing the ideal frequency response of a lowpass filter.)  Moreover, the stopband of the lowpass filter designed in this way will be determined by the tail integral of the window transform, not the window transform itself, provided that the main lobe width is narrow compared with the desired lowpass filter passband.  (Again, think about the convolution operation pictorially to see this directly.)  The tail integral is smaller because the area under the oscillating sidelobes partially cancels out under the integration.  When the desired filter passband is narrow compared with the window transform main lobe, then the lowpass filter frequency response starts to look like a delta function compared with the window transform.  That means the convolution of the two will tend to reproduce the window transform.  In this situation, the transition width is half a main lobe width, the passband with is zero, and the stop-band attenuation is approximately given by the window sidelobe attenuation rather than the tail integral of the sidelobes.  A good reference on the window method for digital filter design is Rabiner and Gold '75: Theory and Application of Digital Signal Processing, Prentice-Hall.

Voice-Quality Audio Example

att3 = 60;	(* desired stop-band attenuation in dB *)

beta = goodBeta[att3]

wu3 = Table[kaiser[n],{n,0,nu}];
wl3 = Table[kaiser[n],{n,nl,-1}];
zp = 0 * Range[NN-M];
w3 = Flatten[{wu3,zp,wl3}];
W3 = fft[w3];
Wdb3 = 20 Log[10,Abs[W3]];

pw3 = ListPlot[Flatten[{wl3,wu3}]];

pWdb3 = ListPlot[Wdb3];

Comparison of Examples

Show[pw2,pw3];

Show[pWdb2,pWdb3];