octave_packages/signal-1.1.3/firls.m

   1 ## Copyright (C) 2006 Quentin Spencer <qspencer@ieee.org>
   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 ## b = firls(N, F, A);
  17 ## b = firls(N, F, A, W);
  18 ##
  19 ##  FIR filter design using least squares method. Returns a length N+1
  20 ##  linear phase filter such that the integral of the weighted mean
  21 ##  squared error in the specified bands is minimized.
  22 ##
  23 ##  F specifies the frequencies of the band edges, normalized so that
  24 ##  half the sample frequency is equal to 1.  Each band is specified by
  25 ##  two frequencies, to the vector must have an even length.
  26 ##
  27 ##  A specifies the amplitude of the desired response at each band edge.
  28 ##
  29 ##  W is an optional weighting function that contains one value for each
  30 ##  band that weights the mean squared error in that band. A must be the
  31 ##  same length as F, and W must be half the length of F.
  32 ##
  33 ## The least squares optimization algorithm for computing FIR filter
  34 ## coefficients is derived in detail in:
  35 ##
  36 ## I. Selesnick, "Linear-Phase FIR Filter Design by Least Squares,"
  37 ## http://cnx.org/content/m10577
  38
  39 function coef = firls(N, frequencies, pass, weight, str);
  40
  41   if nargin<3 || nargin>6
  42     print_usage;
  43   elseif nargin==3
  44     weight = ones(1, length(pass)/2);
  45     str = [];
  46   elseif nargin==4
  47     if ischar(weight)
  48       str = weight;
  49       weight = ones (size (pass));
  50     else
  51       str = [];
  52     end
  53   end
  54   if length (frequencies) ~= length (pass)
  55     error("F and A must have equal lengths.");
  56   elseif 2 * length (weight) ~= length (pass)
  57     error("W must contain one weight per band.");
  58   elseif ischar(str)
  59     error("This feature is implemented yet");
  60   else
  61
  62     M = N/2;
  63     w = kron(weight(:), [-1; 1]);
  64     omega = frequencies * pi;
  65     i1 = 1:2:length(omega);
  66     i2 = 2:2:length(omega);
  67
  68     ## Generate the matrix Q
  69     ## As illustrated in the above-cited reference, the matrix can be
  70     ## expressed as the sum of a Hankel and Toeplitz matrix. A factor of
  71     ## 1/2 has been dropped and the final filter coefficients multiplied
  72     ## by 2 to compensate.
  73     cos_ints = [omega; sin((1:N)' * omega)];
  74     q = [1, 1./(1:N)]' .* (cos_ints * w);
  75     Q = toeplitz (q(1:M+1)) + hankel (q(1:M+1), q(M+1:end));
  76
  77     ## The vector b is derived from solving the integral:
  78     ##
  79     ##           _ w
  80     ##          /   2
  81     ##  b  =   /       W(w) D(w) cos(kw) dw
  82     ##   k    /    w
  83     ##       -      1
  84     ##
  85     ## Since we assume that W(w) is constant over each band (if not, the
  86     ## computation of Q above would be considerably more complex), but
  87     ## D(w) is allowed to be a linear function, in general the function
  88     ## W(w) D(w) is linear. The computations below are derived from the
  89     ## fact that:
  90     ##     _
  91     ##    /                          a              ax + b
  92     ##   /   (ax + b) cos(nx) dx =  --- cos (nx) +  ------ sin(nx)
  93     ##  /                             2                n
  94     ## -                             n
  95     ##
  96     cos_ints2 = [omega(i1).^2 - omega(i2).^2; ...
  97                  cos((1:M)' * omega(i2)) - cos((1:M)' * omega(i1))] ./ ...
  98                  ([2, 1:M]' * (omega(i2) - omega(i1)));
  99     d = [-weight .* pass(i1); weight .* pass(i2)] (:);
 100     b = [1, 1./(1:M)]' .* ((kron (cos_ints2, [1, 1]) + cos_ints(1:M+1,:)) * d);
 101
 102     ## Having computed the components Q and b of the  matrix equation,
 103     ## solve for the filter coefficients.
 104     a = Q \ b;
 105     coef = [ a(end:-1:2); 2 * a(1); a(2:end) ];
 106
 107   end
 108
 109 endfunction