FIR Filter order estimation for Kaiser window design method

[radiolistener]

I'm trying to understand the proper way to estimate FIR filter order for Kaiser window design method.

Here is equation 7.30 from the book "Digital Signal Processing Using MATLAB":


So, we have empirical equation: M = (As - 7.95) / (2.285 * dw) + 1
where
As - stop band attenuation in dB (positive value)
dw - transition band width in angular frequency (radians)

dw = ws - wp

where
wp - passband edge angular frequency
ws - stopband edge angular frequency



Here is also similar equation from the book "Handbook for Digital Signal Processing" Sanjit K.Mitra and James F.Kaiser:


which is: M = (As - 7.95) / (14.36*dw / pi)

converting it for dw in radians, we get the same equation M = (As - 7.95) / (2.285*dw), but there is missing +1.

So, I'm a little bit confused, which one is correct?

Attempt to compare it with MATLAB function kaiserwin shows the significant difference in results:

Code: [Select]

Fs = 2;            % Sample rate [Hz]
Fstop = 0.5;       % Stopband [Hz]
Fpass = 0.4;       % Passband [Hz]
Ap = 1;            % Passband ripple [dB]
dw = 2 * pi * (Fstop-Fpass) / Fs;   % transtion band width in [rad]

for As = 10:10:160
   fir = kaiserwin(fdesign.lowpass('Fp,Fst,Ap,Ast', Fpass, Fstop, Ap, As, Fs)).Numerator;
   M = (As - 7.95) / (2.285 * dw);
   fprintf('As=%.0f dB: %d / %f\n', As, length(fir), M);
end



It looks close, but obviously not the same.

Unfortunately there is no source code for matlab kaiserwin, so I cannot check how it estimates filter order.
Any idea?


The second question is how to take into account passband ripple.
I found the following kaiserWin code implemented in R language:
https://rdrr.io/cran/FIACH/man/kaiserWin.html
https://rdrr.io/cran/FIACH/src/R/kaiserWin.R

Code: [Select]

#
# fh High pass frequency in Hz.
# fl Low pass frequency in Hz
# tw Transition width in Hz.
# sf Sampling frequency in Hz.
# d.sa Desired Stopband Attenuation specified in dB.
# d.pbr Desired Passband ripple specified in dB.
# type string of either "low","high","band" depending on what you wish to do.
#
kaiserWin<-function(fh=NULL,fl=NULL,tw,sf,d.sa,d.pbr,type){
 
  pbr<-10^(d.pbr/20)-1
  sbr<-10^(d.sa/-20)
  ripple<-min(c(pbr,sbr))
  e.sa<- -20*log10(ripple)
 
  df<-tw/sf
  N<-(e.sa-7.95)/(14.36*df)
  if(round(N)%%2==0){N<-round(N+1)}else{N<-round(N)}
 
 
  if(d.sa>=50)          {B<-.1102*(e.sa-8.7)}
  if(d.sa<50 & d.sa>21) {B<-.5842*(e.sa-21)^.4 + .07886*(e.sa-21)}
  if(d.sa<=21)          {B<-0}
 
  n<- -floor(N/2):floor(N/2)
  k<-B*sqrt(1-(2*n/(N-1))^2)
 
  kaiser.win<-besselI(k,0)/besselI(B,0)
  pulse<-sinc(fh=fh,fl=fl,tw=tw,sf=sf,type=type,n=n)
 
  fir<-pulse*kaiser.win
  return(fir)
}
But it confuse me even more.

First, where come from this equation?

Code: [Select]

pbr<-10^(d.pbr/20)-1
And why beta parameter for Kaiser window is calculated using if condition for Ast (d.sa), but calculate it based on combination of Ast and Ap (e.sa)? Is it mistake or there is some sense of that?

Any help or more article links on the subject are welcome.
Thanks

[aliarifat794]

For more detailed understanding, you can see "Digital Signal Processing: Principles, Algorithms, and Applications" by John G. Proakis and Dimitris G. Manolakis.

[gf]

I'm using kaiserord(). The referenced doc also contains examples. And kaiserord() is also available in Octave.

[radiolistener]

I'm using kaiserord(). The referenced doc also contains examples. And kaiserord() is also available in Octave.

kaiserord uses the same equation, but with less precision: n = max(1,ceil((A-8)/(2.285*dw)))

here it's source code:

Code: [Select]

function [n, w, beta, ftype] = kaiserord(f, m, dev, fs)

  if (nargin<2 || nargin>4)
    print_usage;
  endif

  ## default sampling rate parameter
  if nargin<4, fs=2; endif

  ## parameter checking
  if length(f)!=2*length(m)-2
    error("kaiserord must have one magnitude for each frequency band");
  endif
  if any(m(1:length(m)-2)!=m(3:length(m)))
    error("kaiserord pass and stop bands must be strictly alternating");
  endif
  if length(dev)!=length(m) && length(dev)!=1
    error("kaiserord must have one deviation for each frequency band");
  endif
  dev = min(dev);
  if dev <= 0, error("kaiserord must have dev>0"); endif

  ## use midpoints of the transition region for band edges
  w = (f(1:2:length(f))+f(2:2:length(f)))/fs;

  ## determine ftype
  if length(w) == 1
    if m(1)>m(2), ftype='low'; else ftype='high'; endif
  elseif length(w) == 2
    if m(1)>m(2), ftype='stop'; else ftype='pass'; endif
  else
    if m(1)>m(2), ftype='DC-1'; else ftype='DC-0'; endif
  endif

  ## compute beta from dev
  A = -20*log10(dev);
  if (A > 50)
    beta = 0.1102*(A-8.7);
  elseif (A >= 21)
    beta = 0.5842*(A-21)^0.4 + 0.07886*(A-21);
  else
    beta = 0.0;
  endif

  ## compute n from beta and dev
  dw = 2*pi*min(f(2:2:length(f))-f(1:2:length(f)))/fs;
  n = max(1,ceil((A-8)/(2.285*dw)));

  ## if last band is high, make sure the order of the filter is even.
  if ((m(1)>m(2)) == (rem(length(w),2)==0)) && rem(n,2)==1, n = n+1; endif

endfunction
It gives almost the same result as original M = (As - 7.95) / (2.285 * dw), the difference is < 1.
Interesting that matlab kaiserord implementation also gives the same results.

But this estimation is different from one used in matlab kaiserwin(). Why?

[Benta]

Is windowing even used today?
I though Parks-McLellan and Remez exchange had taken over, given the unlimited computing power available now.
Or am I wrong?

[radiolistener]

For more detailed understanding, you can see "Digital Signal Processing: Principles, Algorithms, and Applications" by John G. Proakis and Dimitris G. Manolakis.

Could you please clarify - what this book talks about the subject (FIR order esitmation for Kaiser window design method)?

[gf]

...
converting it for dw in radians, we get the same equation M = (As - 7.95) / (2.285*dw), but there is missing +1.
...

The question is whether M refers to filter order or number of taps. Number of taps = order+1. Kaiserord returns the order, but other functions or equations may return the number of taps.