X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fsignal-1.1.3%2Ffirls.m;fp=octave_packages%2Fsignal-1.1.3%2Ffirls.m;h=e79e6fc7679575a227d658473b79a7b01d4426cb;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05

diff --git a/octave_packages/signal-1.1.3/firls.m b/octave_packages/signal-1.1.3/firls.m
new file mode 100644
index 0000000..e79e6fc
--- /dev/null
+++ b/octave_packages/signal-1.1.3/firls.m
@@ -0,0 +1,109 @@
+## 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.
+##
+## 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 = [];
+    end
+  end
+  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
+
+    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) ];
+
+  end
+
+endfunction