/usr/share/octave/packages/signal-1.2.2/firls.m is in octave-signal 1.2.2-1build1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 | ## Copyright (C) 2006 Quentin Spencer <qspencer@ieee.org>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
## b = firls(N, F, A);
## b = firls(N, F, A, W);
##
## FIR filter design using least squares method. Returns a length N+1
## linear phase filter such that the integral of the weighted mean
## squared error in the specified bands is minimized.
##
## F specifies the frequencies of the band edges, normalized so that
## half the sample frequency is equal to 1. Each band is specified by
## two frequencies, to the vector must have an even length.
##
## A specifies the amplitude of the desired response at each band edge.
##
## W is an optional weighting function that contains one value for each
## band that weights the mean squared error in that band. A must be the
## same length as F, and W must be half the length of F. N must be
## even. If given an odd value, firls increments it by 1.
##
## The least squares optimization algorithm for computing FIR filter
## coefficients is derived in detail in:
##
## I. Selesnick, "Linear-Phase FIR Filter Design by Least Squares,"
## http://cnx.org/content/m10577
function coef = firls(N, frequencies, pass, weight, str);
if (nargin < 3 || nargin > 6)
print_usage;
elseif (nargin == 3)
weight = ones (1, length(pass)/2);
str = [];
elseif (nargin == 4)
if ischar (weight)
str = weight;
weight = ones (size (pass));
else
str = [];
endif
endif
if length (frequencies) ~= length (pass)
error("F and A must have equal lengths.");
elseif 2 * length (weight) ~= length (pass)
error("W must contain one weight per band.");
elseif ischar (str)
error("This feature is implemented yet");
else
N += mod (N, 2);
M = N/2;
w = kron (weight(:), [-1; 1]);
omega = frequencies * pi;
i1 = 1:2:length (omega);
i2 = 2:2:length (omega);
## Generate the matrix Q
## As illustrated in the above-cited reference, the matrix can be
## expressed as the sum of a Hankel and Toeplitz matrix. A factor of
## 1/2 has been dropped and the final filter coefficients multiplied
## by 2 to compensate.
cos_ints = [omega; sin((1:N)' * omega)];
q = [1, 1./(1:N)]' .* (cos_ints * w);
Q = toeplitz (q(1:M+1)) + hankel (q(1:M+1), q(M+1:end));
## The vector b is derived from solving the integral:
##
## _ w
## / 2
## b = / W(w) D(w) cos(kw) dw
## k / w
## - 1
##
## Since we assume that W(w) is constant over each band (if not, the
## computation of Q above would be considerably more complex), but
## D(w) is allowed to be a linear function, in general the function
## W(w) D(w) is linear. The computations below are derived from the
## fact that:
## _
## / a ax + b
## / (ax + b) cos(nx) dx = --- cos (nx) + ------ sin(nx)
## / 2 n
## - n
##
cos_ints2 = [omega(i1).^2 - omega(i2).^2; ...
cos((1:M)' * omega(i2)) - cos((1:M)' * omega(i1))] ./ ...
([2, 1:M]' * (omega(i2) - omega(i1)));
d = [-weight .* pass(i1); weight .* pass(i2)] (:);
b = [1, 1./(1:M)]' .* ((kron (cos_ints2, [1, 1]) + cos_ints(1:M+1,:)) * d);
## Having computed the components Q and b of the matrix equation,
## solve for the filter coefficients.
a = Q \ b;
coef = [ a(end:-1:2); 2 * a(1); a(2:end) ];
endif
endfunction
|