1 ## Copyright (C) 2002 André Carezia <andre@carezia.eng.br>
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
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
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/>.
16 ## Usage: qp_kaiser (nb, at, linear)
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.
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.
28 ## By using a third non-null argument, the minimum-phase calculation is
31 function h = qp_kaiser (nb, at, linear = 0)
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");
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);
50 # size of window (rounded to next odd
52 N = (at - 8) / (2.285*bandwidth);
58 beta = 0.1102 * (at - 8.7);
60 beta = 0.5842 * (at - 21)^0.4 + 0.07886 * (at - 21);
65 # squared in freq. domain
68 # multiplied by ideal lowpass filter
70 hideal = 1/nb * sinc(n/nb);
71 hcomp = wsquared .* hideal;
73 # extract square-root of response and
74 # compute minimum-phase version
76 Hsqr = sqrt(abs(fft(hcomp,Ndft)));
80 h = [fliplr(h) h(1) h];
82 Hmin = Hsqr .* exp(-j*imag(hilbert(log(Hsqr))));
86 # truncate and fix amplitude scale