octave_packages/signal-1.1.3/kaiserord.m

   1 ## Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
   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: [n, Wn, beta, ftype] = kaiserord(f, m, dev [, fs])
  17 ##
  18 ## Returns the parameters needed for fir1 to produce a filter of the
  19 ## desired specification from a kaiser window:
  20 ##       n: order of the filter (length of filter minus 1)
  21 ##       Wn: band edges for use in fir1
  22 ##       beta: parameter for kaiser window of length n+1
  23 ##       ftype: choose between pass and stop bands
  24 ##       b = fir1(n,Wn,kaiser(n+1,beta),ftype,'noscale');
  25 ##
  26 ## f: frequency bands, given as pairs, with the first half of the
  27 ##    first pair assumed to start at 0 and the last half of the last
  28 ##    pair assumed to end at 1.  It is important to separate the
  29 ##    band edges, since narrow transition regions require large order
  30 ##    filters.
  31 ## m: magnitude within each band.  Should be non-zero for pass band
  32 ##    and zero for stop band.  All passbands must have the same
  33 ##    magnitude, or you will get the error that pass and stop bands
  34 ##    must be strictly alternating.
  35 ## dev: deviation within each band.  Since all bands in the resulting
  36 ##    filter have the same deviation, only the minimum deviation is
  37 ##    used.  In this version, a single scalar will work just as well.
  38 ## fs: sampling rate.  Used to convert the frequency specification into
  39 ##    the [0, 1], where 1 corresponds to the Nyquist frequency, fs/2.
  40 ##
  41 ## The Kaiser window parameters n and beta are computed from the
  42 ## relation between ripple (A=-20*log10(dev)) and transition width
  43 ## (dw in radians) discovered empirically by Kaiser:
  44 ##
  45 ##           / 0.1102(A-8.7)                        A > 50
  46 ##    beta = | 0.5842(A-21)^0.4 + 0.07886(A-21)     21 <= A <= 50
  47 ##           \ 0.0                                  A < 21
  48 ##
  49 ##    n = (A-8)/(2.285 dw)
  50 ##
  51 ## Example
  52 ##    [n, w, beta, ftype] = kaiserord([1000,1200], [1,0], [0.05,0.05], 11025);
  53 ##    freqz(fir1(n,w,kaiser(n+1,beta),ftype,'noscale'),1,[],11025);
  54
  55 ## TODO: order is underestimated for the final test case: 2 stop bands.
  56 ## TODO:     octave> ftest("kaiserord") # shows test cases
  57
  58 function [n, w, beta, ftype] = kaiserord(f, m, dev, fs)
  59
  60   if (nargin<2 || nargin>4)
  61     print_usage;
  62   endif
  63
  64   ## default sampling rate parameter
  65   if nargin<4, fs=2; endif
  66
  67   ## parameter checking
  68   if length(f)!=2*length(m)-2
  69     error("kaiserord must have one magnitude for each frequency band");
  70   endif
  71   if any(m(1:length(m)-2)!=m(3:length(m)))
  72     error("kaiserord pass and stop bands must be strictly alternating");
  73   endif
  74   if length(dev)!=length(m) && length(dev)!=1
  75     error("kaiserord must have one deviation for each frequency band");
  76   endif
  77   dev = min(dev);
  78   if dev <= 0, error("kaiserord must have dev>0"); endif
  79
  80   ## use midpoints of the transition region for band edges
  81   w = (f(1:2:length(f))+f(2:2:length(f)))/fs;
  82
  83   ## determine ftype
  84   if length(w) == 1
  85     if m(1)>m(2), ftype='low'; else ftype='high'; endif
  86   elseif length(w) == 2
  87     if m(1)>m(2), ftype='stop'; else ftype='pass'; endif
  88   else
  89     if m(1)>m(2), ftype='DC-1'; else ftype='DC-0'; endif
  90   endif
  91
  92   ## compute beta from dev
  93   A = -20*log10(dev);
  94   if (A > 50)
  95     beta = 0.1102*(A-8.7);
  96   elseif (A >= 21)
  97     beta = 0.5842*(A-21)^0.4 + 0.07886*(A-21);
  98   else
  99     beta = 0.0;
 100   endif
 101
 102   ## compute n from beta and dev
 103   dw = 2*pi*min(f(2:2:length(f))-f(1:2:length(f)))/fs;
 104   n = max(1,ceil((A-8)/(2.285*dw)));
 105
 106   ## if last band is high, make sure the order of the filter is even.
 107   if ((m(1)>m(2)) == (rem(length(w),2)==0)) && rem(n,2)==1, n = n+1; endif
 108 endfunction
 109
 110 %!demo
 111 %! Fs = 11025;
 112 %! for i=1:4
 113 %!   if i==1,
 114 %!     subplot(221); bands=[1200, 1500]; mag=[1, 0]; dev=[0.1, 0.1];
 115 %!   elseif i==2
 116 %!     subplot(222); bands=[1000, 1500]; mag=[0, 1]; dev=[0.1, 0.1];
 117 %!   elseif i==3
 118 %!     subplot(223); bands=[1000, 1200, 3000, 3500]; mag=[0, 1, 0]; dev=0.1;
 119 %!   elseif i==4
 120 %!     subplot(224); bands=100*[10, 13, 15, 20, 30, 33, 35, 40];
 121 %!     mag=[1, 0, 1, 0, 1]; dev=0.05;
 122 %!   endif
 123 %!   [n, w, beta, ftype] = kaiserord(bands, mag, dev, Fs);
 124 %!   d=max(1,fix(n/10));
 125 %!   if mag(length(mag))==1 && rem(d,2)==1, d=d+1; endif
 126 %!   [h, f] = freqz(fir1(n,w,ftype,kaiser(n+1,beta),'noscale'),1,[],Fs);
 127 %!   hm = freqz(fir1(n-d,w,ftype,kaiser(n-d+1,beta),'noscale'),1,[],Fs);
 128 %!   plot(f,abs(hm),sprintf("r;order %d;",n-d), ...
 129 %!        f,abs(h), sprintf("b;order %d;",n));
 130 %!   b = [0, bands, Fs/2]; hold on;
 131 %!   for i=2:2:length(b),
 132 %!     hi=mag(i/2)+dev(1); lo=max(mag(i/2)-dev(1),0);
 133 %!     plot([b(i-1), b(i), b(i), b(i-1), b(i-1)],[hi, hi, lo, lo, hi],"c;;");
 134 %!   endfor; hold off;
 135 %! endfor
 136 %!
 137 %! %--------------------------------------------------------------
 138 %! % A filter meets the specifications if its frequency response
 139 %! % passes through the ends of the criteria boxes, and fails if
 140 %! % it passes through the top or the bottom.  The criteria are
 141 %! % met precisely if the frequency response only passes through
 142 %! % the corners of the boxes.  The blue line is the filter order
 143 %! % returned by kaiserord, and the red line is some lower filter
 144 %! % order.  Confirm that the blue filter meets the criteria and
 145 %! % the red line fails.
 146
 147 %!# XXX FIXME XXX extend demo to show detail at criteria box corners