octave_packages/signal-1.1.3/qp_kaiser.m

   1 ## Copyright (C) 2002 André Carezia <andre@carezia.eng.br>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## Usage:  qp_kaiser (nb, at, linear)
  17 ##
  18 ## Computes a finite impulse response (FIR) filter for use with a
  19 ## quasi-perfect reconstruction polyphase-network filter bank. This
  20 ## version utilizes a Kaiser window to shape the frequency response of
  21 ## the designed filter. Tha number nb of bands and the desired
  22 ## attenuation at in the stop-band are given as parameters.
  23 ##
  24 ## The Kaiser window is multiplied by the ideal impulse response
  25 ## h(n)=a.sinc(a.n) and converted to its minimum-phase version by means
  26 ## of a Hilbert transform.
  27 ##
  28 ## By using a third non-null argument, the minimum-phase calculation is
  29 ## ommited at all.
  30
  31 function h = qp_kaiser (nb, at, linear = 0)
  32
  33   if (nargin < 2)
  34     print_usage;
  35   elseif !(isscalar (nb) && (nb == round(nb)) && (nb >= 0))
  36     error ("qp_kaiser: nb has to be a positive integer");
  37   elseif !(isscalar (at) && (at == real (at)))
  38     error ("qp_kaiser: at has to be a real constant");
  39   endif
  40
  41                                 # Bandwidth
  42   bandwidth = pi/nb;
  43
  44                                 # Attenuation correction (empirically
  45                                 # determined by M. Gerken
  46                                 # <mgk@lcs.poli.usp.br>)
  47   corr = (1.4+0.6*(at-20)/80)^(20/at);
  48   at = corr * at;
  49
  50                                 # size of window (rounded to next odd
  51                                 # integer)
  52   N = (at - 8) / (2.285*bandwidth);
  53   M = fix(N/2);
  54   N = 2*M + 1;
  55
  56                                 # Kaiser window
  57   if (at>50)
  58     beta = 0.1102 * (at - 8.7);
  59   elseif (at>21)
  60     beta = 0.5842 * (at - 21)^0.4 + 0.07886 * (at - 21);
  61   else
  62     beta = 0;
  63   endif
  64   w = kaiser(N,beta);
  65                                 # squared in freq. domain
  66   wsquared = conv(w,w);
  67
  68                                 # multiplied by ideal lowpass filter
  69   n = -(N-1):(N-1);
  70   hideal = 1/nb * sinc(n/nb);
  71   hcomp = wsquared .* hideal;
  72
  73                                 # extract square-root of response and
  74                                 # compute minimum-phase version
  75   Ndft = 2^15;
  76   Hsqr = sqrt(abs(fft(hcomp,Ndft)));
  77   if (linear)
  78     h = real(ifft(Hsqr));
  79     h = h(2:N);
  80     h = [fliplr(h) h(1) h];
  81   else
  82     Hmin = Hsqr .* exp(-j*imag(hilbert(log(Hsqr))));
  83     h = real(ifft(Hmin));
  84     h = h(1:N);
  85   endif
  86                                 # truncate and fix amplitude scale
  87                                 # (H(0)=1)
  88   h = h / sum(h);
  89
  90 endfunction