1 ## Copyright (C) 2006 Quentin Spencer <qspencer@ieee.org>
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 ## b = firls(N, F, A);
17 ## b = firls(N, F, A, W);
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.
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.
27 ## A specifies the amplitude of the desired response at each band edge.
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.
33 ## The least squares optimization algorithm for computing FIR filter
34 ## coefficients is derived in detail in:
36 ## I. Selesnick, "Linear-Phase FIR Filter Design by Least Squares,"
37 ## http://cnx.org/content/m10577
39 function coef = firls(N, frequencies, pass, weight, str);
41 if nargin<3 || nargin>6
44 weight = ones(1, length(pass)/2);
49 weight = ones (size (pass));
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.");
59 error("This feature is implemented yet");
63 w = kron(weight(:), [-1; 1]);
64 omega = frequencies * pi;
65 i1 = 1:2:length(omega);
66 i2 = 2:2:length(omega);
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));
77 ## The vector b is derived from solving the integral:
81 ## b = / W(w) D(w) cos(kw) dw
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
92 ## / (ax + b) cos(nx) dx = --- cos (nx) + ------ sin(nx)
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);
102 ## Having computed the components Q and b of the matrix equation,
103 ## solve for the filter coefficients.
105 coef = [ a(end:-1:2); 2 * a(1); a(2:end) ];