octave_packages/signal-1.1.3/chebwin.m

   1 ## Copyright (C) 2002 André Carezia <acarezia@uol.com.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:  chebwin (L, at)
  17 ##
  18 ## Returns the filter coefficients of the L-point Dolph-Chebyshev window
  19 ## with at dB of attenuation in the stop-band of the corresponding
  20 ## Fourier transform.
  21 ##
  22 ## For the definition of the Chebyshev window, see
  23 ##
  24 ## * Peter Lynch, "The Dolph-Chebyshev Window: A Simple Optimal Filter",
  25 ##   Monthly Weather Review, Vol. 125, pp. 655-660, April 1997.
  26 ##   (http://www.maths.tcd.ie/~plynch/Publications/Dolph.pdf)
  27 ##
  28 ## * C. Dolph, "A current distribution for broadside arrays which
  29 ##   optimizes the relationship between beam width and side-lobe level",
  30 ##   Proc. IEEE, 34, pp. 335-348.
  31 ##
  32 ## The window is described in frequency domain by the expression:
  33 ##
  34 ##          Cheb(L-1, beta * cos(pi * k/L))
  35 ##   W(k) = -------------------------------
  36 ##                 Cheb(L-1, beta)
  37 ##
  38 ## with
  39 ##
  40 ##   beta = cosh(1/(L-1) * acosh(10^(at/20))
  41 ##
  42 ## and Cheb(m,x) denoting the m-th order Chebyshev polynomial calculated
  43 ## at the point x.
  44 ##
  45 ## Note that the denominator in W(k) above is not computed, and after
  46 ## the inverse Fourier transform the window is scaled by making its
  47 ## maximum value unitary.
  48 ##
  49 ## See also: kaiser
  50
  51 function w = chebwin (L, at)
  52
  53   if (nargin != 2)
  54     print_usage;
  55   elseif !(isscalar (L) && (L == round(L)) && (L > 0))
  56     error ("chebwin: L has to be a positive integer");
  57   elseif !(isscalar (at) && (at == real (at)))
  58     error ("chebwin: at has to be a real scalar");
  59   endif
  60
  61   if (L == 1)
  62     w = 1;
  63   else
  64                                 # beta calculation
  65     gamma = 10^(-at/20);
  66     beta = cosh(1/(L-1) * acosh(1/gamma));
  67                                 # freq. scale
  68     k = (0:L-1);
  69     x = beta*cos(pi*k/L);
  70                                 # Chebyshev window (freq. domain)
  71     p = cheb(L-1, x);
  72                                 # inverse Fourier transform
  73     if (rem(L,2))
  74       w = real(fft(p));
  75       M = (L+1)/2;
  76       w = w(1:M)/w(1);
  77       w = [w(M:-1:2) w]';
  78     else
  79                                 # half-sample delay (even order)
  80       p = p.*exp(j*pi/L * (0:L-1));
  81       w = real(fft(p));
  82       M = L/2+1;
  83       w = w/w(2);
  84       w = [w(M:-1:2) w(2:M)]';
  85     endif
  86   endif
  87
  88   w = w ./ max (w (:));
  89 endfunction