Re: Possible contribution for the "signal" package

[Søren Hauberg]

Hi

man, 28 09 2009 kl. 08:23 -0400, skrev Pete Gonzalez:
> We are offering to contribute an original implementation of the filter
> design algorithm described in this paper:
>
> I. W. Selesnick, M. Lang, and C. S. Burrus. A modified algorithm for
> constrained least square design of multiband FIR filters without
> specified transition bands. IEEE Trans. on Signal Processing,
> 46(2):497-501, February 1998.

Cool! I'm not really a signal processing guy (I do some image
processing, though), so I can't really comment on the actual algorithm.
That being said, I do think it sounds like a valuable contribution to
the 'signal' package.

Since the function should be part of a package, I think it should be
discussed at the Octave-Forge mailing list. So, I've moved the
discussion there.

As to the actual code, then I only have few comments/questions:

*) Why do you need the 'MallocArray' class? Can't you just use
   the 'Matrix' or 'ColumnVector' class?

*) You need to do some error checking in the input handling.
   Basically, you need to some

    if (error_state)
      {
        error ("something");
        return octave_value ();
      }
   
   after each call to '*_value ()' with '*' being e.g. 'double'.

Otherwise I guess things look good :-)

Søren


------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Pete Gonzalez]

Thanks Søren.  I have subscribed to the octave-dev list.  I will make
the changes you suggested and post an updated version.

Cheers,
-Pete


Søren Hauberg wrote:

> Hi
>
> man, 28 09 2009 kl. 08:23 -0400, skrev Pete Gonzalez:
>> We are offering to contribute an original implementation of the filter
>> design algorithm described in this paper:
>>
>> I. W. Selesnick, M. Lang, and C. S. Burrus. A modified algorithm for
>> constrained least square design of multiband FIR filters without
>> specified transition bands. IEEE Trans. on Signal Processing,
>> 46(2):497-501, February 1998.
>
> Cool! I'm not really a signal processing guy (I do some image
> processing, though), so I can't really comment on the actual algorithm.
> That being said, I do think it sounds like a valuable contribution to
> the 'signal' package.
>
> Since the function should be part of a package, I think it should be
> discussed at the Octave-Forge mailing list. So, I've moved the
> discussion there.
>
> As to the actual code, then I only have few comments/questions:
>
> *) Why do you need the 'MallocArray' class? Can't you just use
>    the 'Matrix' or 'ColumnVector' class?
>
> *) You need to do some error checking in the input handling.
>    Basically, you need to some
>
>     if (error_state)
>       {
>         error ("something");
>         return octave_value ();
>       }
>    
>    after each call to '*_value ()' with '*' being e.g. 'double'.
>
> Otherwise I guess things look good :-)
>
> Søren
>
>

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Pete Gonzalez]

Hi Søren,

Attached is an updated version of cl2bp.cc with the following changes:

- As requested, I added error_state checks for the input parameters.
BTW I noticed that a lot of shipping Octave code is missing these
checks; it should probably be audited at some point.

- The MallocArray class you asked about was intended to support LGPL
scenarios, i.e. where people can't use the GPL'ed oct.h.  That code is
now isolated with a CL2BP_STANDALONE switch.  In a normal build, the
native ColumnVector/Array classes are now used, but with C-style
indexing operators.  (I didn't convert x[i] to x(i) everywhere because I
think it would hinder readability.)

- We added a CL2BP_LOGGING switch for enabling/disabling diagnostics,
which is useful when testing for regressions

- We measured the speed of long-running calculations and inserted
OCTAVE_QUIT in the appropriate inner loops

- The input parameter checks have been adjusted

Let me know if you need anything else.  Thanks!

Cheers,
-Pete


Søren Hauberg wrote:

> Hi
>
> man, 28 09 2009 kl. 08:23 -0400, skrev Pete Gonzalez:
>> We are offering to contribute an original implementation of the filter
>> design algorithm described in this paper:
>>
>> I. W. Selesnick, M. Lang, and C. S. Burrus. A modified algorithm for
>> constrained least square design of multiband FIR filters without
>> specified transition bands. IEEE Trans. on Signal Processing,
>> 46(2):497-501, February 1998.
>
> Cool! I'm not really a signal processing guy (I do some image
> processing, though), so I can't really comment on the actual algorithm.
> That being said, I do think it sounds like a valuable contribution to
> the 'signal' package.
>
> Since the function should be part of a package, I think it should be
> discussed at the Octave-Forge mailing list. So, I've moved the
> discussion there.
>
> As to the actual code, then I only have few comments/questions:
>
> *) Why do you need the 'MallocArray' class? Can't you just use
>    the 'Matrix' or 'ColumnVector' class?
>
> *) You need to do some error checking in the input handling.
>    Basically, you need to some
>
>     if (error_state)
>       {
>         error ("something");
>         return octave_value ();
>       }
>    
>    after each call to '*_value ()' with '*' being e.g. 'double'.
>
> Otherwise I guess things look good :-)
>
> Søren
>
>

/*

Copyright (c) 2008-2009, Evgeni A. Nurminski
Copyright (c) 2008-2009, Pete Gonzalez

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.


Authors:

Evgeni A. Nurminski <nurmi@...>
Institute for Automation and Control Problems
Far East branch of RAS, Vladivostok, Russia

Pete Gonzalez <pgonzalez@...>
Bluel Technologies Corporation

*/

#ifndef CL2BP_STANDALONE  // Define this symbol for usage outside of Octave
#include <octave/oct.h>
#endif

//===========================================================================================================

#ifdef _MSC_VER
#define _CRT_SECURE_NO_WARNINGS
#define _USE_MATH_DEFINES
#include <windows.h>  // OutputDebugStringA()
#endif

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <math.h>
#include <assert.h>
#include <string.h>

//-----------------------------------------------------------------------------------------------------------
// Logging
//-----------------------------------------------------------------------------------------------------------

#ifdef CL2BP_LOGGING // Define this symbol for diagnostic messages
static void debug_print(const char *format, ...) {
  va_list argptr;
  va_start(argptr, format);

  char buf[20*1024];
  if (vsnprintf(buf,20*1024-1,format,argptr) == -1) {
    strcpy(buf, "#ERROR#");
  }
  va_end(argptr);

#ifdef _MSC_VER
  OutputDebugStringA(buf);            // log to the Visual Studio output pane
#else
  octave_stdout << buf;               // log to the Octave pager
#endif
}
#else
#define debug_print(...)  ((void)0)   // don't log anything (for a big speed improvement)
#endif

//-----------------------------------------------------------------------------------------------------------
// Array wrappers
//-----------------------------------------------------------------------------------------------------------

#ifndef CL2BP_STANDALONE
// These extend the native Octave classes with standard C-style operators to improve readability.
template <class T>
class ArrayC : public Array<T> {
public:
  T& operator[](int index) { return Array<T>::xelem(index); }
  T operator[](int index) const { return Array<T>::xelem(index); }
  ArrayC(octave_idx_type n, const T& fill_value) : Array<T>(n,fill_value) { }
  ~ArrayC() { }
};

class ColumnVectorC : public ColumnVector {
public:
  double& operator[](int index) { return ColumnVector::xelem(index); }
  double operator[](int index) const { return ColumnVector::xelem(index); }
  ColumnVectorC(octave_idx_type n, double fill_value) : ColumnVector(n, fill_value) { }
  ColumnVectorC() : ColumnVector() { }
  ~ColumnVectorC() { }
};
#else
// This is a trivial replacement for Octave's Array/ColumnVector classes, for usage in LGPL scenarios
template <class T>
class ArrayC {
  int count;
  T *ptr;
public:
  T& operator[](int index) {
    assert(index >= 0 && index < count);
    return ptr[index];
  }
  T operator[](int index) const {
    assert(index >= 0 && index < count);
    return ptr[index];
  }
  int length() const { return count; }
  void resize(int count_, T fill_value) {
    assert(count_ >= 0 && count_ <= 512*1024*1024);  // verify that the array size is reasonable
    count = count_;
    ptr = (T *)realloc(ptr, count * sizeof(T));
    for (int i=0; i<count; ++i) ptr[i] = fill_value;
  }
  ArrayC(int count_, T fill_value) {
    ptr = 0;
    resize(count_, fill_value);
  }
  ~ArrayC() { free(ptr); }
private:
  ArrayC(const ArrayC& src) { }  // copy constructor is unimplemented and disallowed
};

class ColumnVectorC : public ArrayC<double> {
public:
  ColumnVectorC(int count_, double fill_value) : ArrayC<double>(count_, fill_value) { }
  ColumnVectorC() : ArrayC<double>(0,0.0) { }
  ~ColumnVectorC() { }
};
#define OCTAVE_QUIT ((void)0)
#endif

//-----------------------------------------------------------------------------------------------------------
// Cl2bp Algorithm
//-----------------------------------------------------------------------------------------------------------

#define BIG_NUMBER 100000

//-----------------------------------------------------------------------------------------------------------
static int local_max(const ColumnVectorC& x, int n, ArrayC<int>& ix) {
  int i, mx;

  mx = 0;

  ColumnVectorC xx(n+2, 0.0);

  xx[0] = xx[n + 1] = -BIG_NUMBER;
  for ( i = 1; i <= n; i++)
    xx[i] = x[i - 1];
  for ( i = 0; i < n; i++ ) {
    (((xx[i] <  xx[i + 1]) && (xx[i + 1] >  xx[i + 2])) ||
     ((xx[i] <  xx[i + 1]) && (xx[i + 1] >= xx[i + 2])) ||
     (xx[i] <= xx[i + 1]) && (xx[i + 1] >  xx[i + 2]) )
    && ( ix[mx++] = i );
  }
  return mx;
}

//-----------------------------------------------------------------------------------------------------------
static int solvps(ColumnVectorC& a, ColumnVectorC& b, int n) {
  double t;
  int j, k;
  int a_p;
  int a_q;
  int a_r;
  int a_s;

  OCTAVE_QUIT;  // In empirical tests, 2% to 6% of the execution time was spent in solvps()

  for (j = 0, a_p = 0; j < n; ++j, a_p += n+1) {
    for (a_q = j*n; a_q < a_p; ++a_q)
      a[a_p] -= a[a_q] * a[a_q];
    if (a[a_p] <= 0.)
      return -1;
    a[a_p] = sqrt(a[a_p]);
    for (k = j + 1, a_q = a_p + n; k < n; ++k, a_q += n) {
      for (a_r = j*n, a_s = k*n, t = 0.; a_r < a_p;)
        t += a[a_r++] * a[a_s++];
      a[a_q] -= t;
      a[a_q] /= a[a_p];
    }
  }
  for (j = 0, a_p = 0; j < n; ++j, a_p += n+1) {
    for (k = 0, a_q = j*n; k < j;)
      b[j] -=b [k++] * a[a_q++];
    b[j] /= a[a_p];
  }
  for (j = n - 1, a_p = n*n - 1; j >= 0; --j, a_p -= n + 1) {
    for (k = j + 1, a_q = a_p + n; k < n; a_q += n)
      b[j] -= b[k++]* a[a_q];
    b[j] /= a[a_p];
  }
  return 0;
}

//-----------------------------------------------------------------------------------------------------------
static void rmmult(ColumnVectorC& rm, const ColumnVectorC& a, const ColumnVectorC& b,
  int n, int m, int l) {

  double z;
  int i,j,k;
  int b_p; // index into b
  int a_p; // index into a
  int rm_start = 0;  // index into rm
  int rm_q; // index into rm

  ColumnVectorC q0(m, 0.0);
  for (i = 0; i < l ; ++i, ++rm_start) {
    OCTAVE_QUIT;  // In empirical tests, 87% to 95% of the execution time was spent in rmmult()
    for (k = 0, b_p = i; k < m; b_p += l)
      q0[k++] = b[b_p];
    for (j = 0, a_p = 0, rm_q = rm_start; j < n; ++j, rm_q += l) {
      for (k = 0, z = 0.; k < m;)
        z += a[a_p++] * q0[k++];
      rm[rm_q]=z;
    }
  }
}

//-----------------------------------------------------------------------------------------------------------
static void mattr(ColumnVectorC& a, const ColumnVectorC& b, int m, int n) {
  int i, j;
  int b_p;
  int a_p = 0;
  int b_start = 0;
  for (i = 0; i < n; ++i, ++b_start)
    for ( j = 0, b_p = b_start; j < m; ++j, b_p += n )
      a[a_p++] = b[b_p];
}

//-----------------------------------------------------------------------------------------------------------
static void calcAx(ColumnVectorC& Ax, int m, int L) {
  double r = M_SQRT2, pi = M_PI;
  int i, j;

  Ax.resize((L+1)*(m + 1), 0.0);

  for ( i = 0; i <= L; i++)
    for ( j = 0; j <= m; j++)
      Ax[i*(m+1) + j] = cos(i*j*pi/L);
  for ( i = 0; i < (L+1)*(m+1); i += m + 1 )
    Ax[i] /= r;
}

//-----------------------------------------------------------------------------------------------------------
static double L2error(const ColumnVectorC& x, const ColumnVectorC& w, int L, double w1, double w2) {
  double xx, ww, sumerr, pi = M_PI;
  int i;
  for ( i = 0, sumerr = 0; i < L + 1; i++ ) {
    ww = w[i];
    xx = x[i];
    sumerr += ( ww < w1*pi || ww > w2*pi ) ?  xx*xx : (1 - xx)*(1 - xx);
  }
  return sumerr;
}

//-----------------------------------------------------------------------------------------------------------
static double CHerror(double *wmin, const ColumnVectorC& x, const ColumnVectorC& w,
  int L, double w1, double w2) {

  double xx, ww, err, errmax;
  int i;
  errmax = -1;
  *wmin = 0;
  for ( i = 0; i < L + 1; i++ ) {
    ww = w[i];
    xx = x[i];
    err = (( ww < w1 ) || (ww > w2 )) ?  fabs(xx) : fabs(1 - xx);
    if ( err > errmax ) {
      errmax = err;
      *wmin = ww;
    }
  }
  return errmax;
}

//-----------------------------------------------------------------------------------------------------------
static void Ggen(ColumnVectorC& G, int m, const ColumnVectorC& w, const ArrayC<int>& kmax,
  int l_kmax, const ArrayC<int>& kmin, int l_kmin) {

  G.resize((l_kmax + l_kmin)*(m + 1), 0.0);

  int nsys, i, j;
  double r = M_SQRT2;

  nsys = l_kmax + l_kmin;
  for ( i = 0; i < l_kmax; i++)
    for ( j = 0; j <= m; j++)
      G[i*(m+1) + j] = cos(w[kmax[i]]*j);
  for ( i = l_kmax; i < nsys; i++)
    for ( j = 0; j <= m; j++)
      G[i*(m+1) + j] = -cos(w[kmin[i - l_kmax]]*j);
  for ( i = 0; i < nsys*(m+1); i += m+1 )
    G[i] /= r;
}

//-----------------------------------------------------------------------------------------------------------
// Constrained L2 BandPass FIR filter design
//  h       2*m+1 filter coefficients (output)
//  m       degree of cosine polynomial
//  w1,w2   fist and second band edges
//  up      upper bounds
//  lo      lower bounds
//  L       grid size
//  eps     stopping condition ("SN")
//  mit     maximum number of iterations
// cl2bp() returns true if the solver converged, or false if the maximum number of iterations was reached.
static bool cl2bp(ColumnVectorC& h, int m, double w1, double w2, double up[3], double lo[3], int L,
  double eps, int mit) {

  double r = M_SQRT2;
  double pi = M_PI;
  double wmin, ww, xw;
  int q1, q2, q3, i, iter = 0;
  int off, neg;

  int ik;
  int l_kmax;
  int l_kmin;
  int l_okmax;
  int l_okmin;
  double uvo = 0, lvo = 0;

  double err, diff_up, diff_lo;
  double erru, erro;
  int iup, ilo;

  int nsys, j;

  int imin;
  double umin;

  int k1, k2, ak1, ak2;
  double cerr;

  h.resize(2*m+1, 0.0);

  bool converged = true;

  ColumnVectorC x0(L+1, 0.0);
  ColumnVectorC x1(L+1, 0.0);
  ColumnVectorC xx(L+1, 0.0);
  ColumnVectorC xm1(L+1, 0.0);
  ColumnVectorC w(L+1, 0.0);
  ColumnVectorC u(L+1, 0.0);
  ColumnVectorC l(L+1, 0.0);
  ColumnVectorC a(L+1, 0.0);
  ColumnVectorC c(L+1, 0.0);
  ArrayC<int> kmax(L+1, 0);
  ArrayC<int> kmin(L+1, 0);
  l_kmax  = l_kmin = 0;

  ArrayC<int> okmin(L+1, 0);
  ArrayC<int> okmax(L+1, 0);
  l_okmax = l_okmin = 0;
  ColumnVectorC rhs(m+2, 0.0);
  ColumnVectorC mu(2*(L+1), 0.0);

  for ( i = 0; i <= L; i++ )
    w[i] = i*pi/L;

  ColumnVectorC Z(2*L-1-2*m, 0.0);

  q1 = (int)floor(L*w1/pi);
  q2 = (int)floor(L*(w2-w1)/pi);
  q3 = L + 1 - q1 - q2;

  off = 0;
  for ( i = 0; i < q1; i++) {
    u[i] = up[0];
    l[i] = lo[0];
  }
  off += i;
  for ( i = 0; i < q2; i++) {
    u[off + i] = up[1];
    l[off + i] = lo[1];
  }
  off += i;
  for ( i = 0; i < q3; i++) {
    u[off + i] = up[2];
    l[off + i] = lo[2];
  }


  c[0] = (w2-w1)*r;
  for ( i = 1; i <= m; i++)
    c[i] =  2*(sin(w2*i)-sin(w1*i))/i;
  for ( i = 0; i <= m; i++) {
    c[i] /=  pi;
    a[i] = c[i];
  }
  debug_print("\nInitial approximation. Unconstrained L2 filter coefficients.\n");
  debug_print("=============================================================\n");

  debug_print("\nZero order term %8.3lf\n", a[0]);
  for ( i = 1; i <= m; i++) {
    debug_print("%4d %8.3lf", i, a[i]);
    if (i - 8*(i/8) == 0)
      debug_print("\n");
  }
  debug_print("\n");

  // Precalculate Ax
  ColumnVectorC Ax;
  calcAx(Ax, m, L);

  //calcA(x0, a, m, L);
  rmmult(x0, Ax, a, L + 1, m + 1, 1);

  err = CHerror(&wmin, x0, w, L, w1, w2);
  debug_print("\nChebyshev err %12.4e (%11.5e)  <> L2 err %12.4e\n", err, wmin/pi, L2error(x0, w, L, w1, w2)/(L+1));
  for (iter = 1; ; iter++) {
    debug_print("\n:::::: Iteration %3d :::::::::::\n", iter);

    if ( (uvo > eps*2) || (lvo > eps*2) ) {
      debug_print("\nXXX Take care of old errors: uvo lvo %12.3e %12.3e", uvo, lvo);
      if( k1 >= 0 ) debug_print(" k1 %3d(%d)", k1, okmax[k1]);
      if( k2 >= 0 ) debug_print(" k2 %3d(%d)", k2, okmin[k2]);
      debug_print("\n");

      if ( uvo > lvo ) {
        kmax[l_kmax] = okmax[k1];
        l_kmax++;
        l_okmax--;
        for (i = k1; i < l_okmax; i++)
          okmax[i] = okmax[i + 1];
      } else {
        kmin[l_kmin] = okmin[k2];
        l_kmin++;
        l_okmin--;
        for (i = k2; i < l_okmin; i++)
          okmin[i] = okmin[i + 1];
      }
      nsys = l_kmax + l_kmin;

      /*
        for (i = 0; i < l_kmax; i++)
          debug_print("i %3d kmax %3d mu %12.4e\n",
             i, kmax[i], mu[i]);
        debug_print("\n");
        for (i = 0; i < l_kmin; i++)
          debug_print("i %3d kmin %3d mu %12.4e\n",
             i, kmin[i], mu[i + l_kmax]);
        debug_print("\n\n");
      */
    } else {

      //calcA(x1, a, m, L);
      rmmult(x1, Ax, a, L + 1, m + 1, 1);


      for ( i = 0; i < l_kmax; i++)
        okmax[i] = kmax[i];
      for ( i = 0; i < l_kmin; i++)
        okmin[i] = kmin[i];
      l_okmax = l_kmax;
      l_okmin = l_kmin;

      l_kmax = local_max(x1, L + 1, kmax);


      for ( i = 0; i < L + 1; i++)
        xm1[i] = -x1[i];
      l_kmin = local_max(xm1, L + 1, kmin);

      debug_print("\nSignificant deviations from the ideal filter. Levels:");
      debug_print(" %8.2e %8.2e %8.2e (lo)", lo[0], lo[1], lo[2]);
      debug_print(" %8.2e %8.2e %8.2e (up)", up[0], up[1], up[2]);
      debug_print("\n");

      for ( i = 0, ik = 0; i < l_kmax; i++) {
        j = kmax[i];
        if ( x1[j] > u[j] - 10*eps ) {
          kmax[ik] = j;
          ik++;
        }
      }
      l_kmax = ik;

      debug_print("overshoots = %d\n", l_kmax);
      for ( i = 0; i < l_kmax; i++) {
        j = kmax[i];
        ww = w[j];
        xw = x1[j];
        err = (w1 < ww && ww < w2) ? xw - 1 : xw;
        debug_print(" i = %3d kmax = %3d xw = %12.5e err = %10.3e u = %12.4e max = %9.2e\n",
               i, j, xw, err, u[j], xw - u[j]);
      }

      for ( i = 0, ik = 0; i < l_kmin; i++) {
        j = kmin[i];
        if ( x1[j] < l[j] + 10*eps ) {
          kmin[ik] = j;
          ik++;
        }
      }
      l_kmin = ik;

      debug_print("undershoots = %d\n", l_kmin);
      for ( i = 0; i < l_kmin; i++) {
        j = kmin[i];
        ww = w[j];
        xw = -xm1[j];
        err =(w1 < ww && ww < w2) ? xw - 1 : xw;
        debug_print(" i = %3d kmin = %3d xw = %12.5e err = %10.3e l = %12.4e min = %9.2e\n",
               i, j, xw, err, l[j], xw - l[j]);
      }

      err = erru = erro = diff_up = diff_lo = 0;
      iup = ilo = 0;
      for ( i = 0; i < l_kmax; i++) {
        ik = kmax[i];
        diff_up = x1[ik] - u[ik];
        if ( diff_up > erru ) {
          erru = diff_up;
          iup = i;
        }
      }
      for ( i = 0; i < l_kmin; i++) {
        ik = kmin[i];
        diff_lo = l[ik] - x1[ik];
        if ( diff_lo > erro ) {
          erro = diff_lo;
          ilo = i;
        }
      }
      err = erro > erru ? erro : erru;
      debug_print("erru = %14.6e iup = %4d kmax = %4d ", erru, iup, kmax[iup]);
      debug_print("erro = %14.6e ilo = %4d kmin = %4d\n", erro, ilo, kmin[ilo]);

      if ( err < eps )
        break;
    }


    nsys = l_kmax + l_kmin;

    ColumnVectorC G(nsys*(m + 1), 0.0);
    ColumnVectorC GT(nsys*(m + 1), 0.0);
    ColumnVectorC GG(nsys*nsys, 0.0);

    for ( i = 0; i < l_kmax; i++)
      for ( j = 0; j <= m; j++)
        G[i*(m+1) + j] = cos(w[kmax[i]]*j);
    for ( i = l_kmax; i < nsys; i++)
      for ( j = 0; j <= m; j++)
        G[i*(m+1) + j] = -cos(w[kmin[i - l_kmax]]*j);
    for ( i = 0; i < nsys*(m+1); i += m+1 )
      G[i] /= r;
    mattr(GT, G, nsys, m + 1);
    rmmult(GG, G, GT, nsys, m + 1, nsys);


    rmmult(rhs, G, c, nsys, m + 1, 1);
    for ( i = 0; i < nsys; i++)
      if ( i < l_kmax )
        rhs[i] -= u[kmax[i]];
      else
        rhs[i] += l[kmin[i - l_kmax]];

    for ( i = 0; i < nsys; ++i)
      mu[i] = rhs[i];

    solvps(GG, mu, nsys);
    debug_print("\nXXX KT-system with %d equations resolved.\n", nsys);
    for ( i = 0, neg = 0; i < nsys; i++)
      if ( mu[i] < 0 ) neg++;
    debug_print("\nTotal number of negative multipliers = %3d\n\n", neg);
    while (neg) {


      umin = BIG_NUMBER;
      for ( i = 0, j = 0; i < nsys; i++) {
        if ( mu[i] >= 0 ) continue;
        j++;
        if ( mu[i] < umin ) {
          imin = i;
          umin = mu[i];
        }
      }

      if ( umin > 0 )
        break;
      debug_print(" neg = %3d of %3d imin = %3d mu-min = %12.4e\n", j, nsys, imin, umin);

      for ( i = imin; i < nsys - 1; i++) {
        rhs[i] = rhs[i + 1];
        for ( j = 0; j <= m; j++)
          G[i*(m + 1) + j] =  G[(i + 1)*(m + 1) + j];
      }

      if ( imin < l_kmax ) {
        for ( i = imin; i < l_kmax - 1; i++)
          kmax[i] = kmax[i + 1];
        l_kmax--;
      } else {
        for ( i = imin; i < nsys - 1; i++)
          kmin[i - l_kmax] = kmin[i - l_kmax + 1];
        l_kmin--;
      }

      --nsys;

      mattr(GT, G, nsys, m + 1);
      rmmult(GG, G, GT, nsys, m + 1, nsys);
      for ( i = 0; i < nsys; ++i)
        mu[i] = rhs[i];
      solvps(GG, mu, nsys);
    }

    ColumnVectorC da(m + 1, 0.0);
    ColumnVectorC zd(nsys, 0.0);

    rmmult(da, GT, mu, m + 1, nsys, 1);
    for ( i = 0; i <= m; i++)
      a[i] = c[i] - da[i];
    rmmult(zd, G, a, nsys, m + 1, 1);

    ColumnVectorC zz(l_okmax + l_okmin, 0.0);
    Ggen(G, m, w, okmax, l_okmax, okmin, l_okmin);
    rmmult(zz, G, a, l_okmax + l_okmin, m + 1, 1);
    uvo = lvo = eps;
    k1 = k2 = -1;
    ak1 = ak2 = -1;
    if (l_okmax + l_okmin > 0)
      debug_print("\nErrors on the previous set of freqs\n\n");
    for (i = 0; i < l_okmax; i++) {
      j = okmax[i];
      cerr = zz[i] - u[j];
      debug_print(" i %2d j %3d u %12.4e Ga %12.4e cerr %12.4e\n",
             i, j, u[j], zz[i], cerr);
      if ( cerr > uvo ) {
        uvo = cerr;
        k1 = i;
        ak1 = j;
      }
    }
    cerr = 0;
    for (i = 0; i < l_okmin; i++) {
      j = okmin[i];
      cerr = l[j] + zz[i + l_okmax];
      debug_print(" i %2d j %3d l %12.4e Ga %12.4e cerr %12.4e\n",
             i, j, l[j], zz[i + l_okmax], cerr);
      if ( cerr > lvo ) {
        lvo = cerr;
        k2 = i, ak2 = j;
      }
    }
    if ( l_okmax + l_okmin > 0 ) {
      debug_print("\n uvo = %12.4e k1 = %4d (%3d) ", uvo, k1, ak1);
      debug_print("  lvo = %12.4e k2 = %4d (%3d) ", lvo, k2, ak2);
      debug_print(" maxerr = %12.4e\n", uvo > lvo ? uvo : lvo);
    }

    debug_print("\nConstrained L2 band filter coefficients.\n");
    debug_print("=====================================\n");

#ifdef CL2BP_LOGGING
    debug_print("\nZero order term %8.3lf\n", a[0]);
    for ( i = 1; i <= m; i++) {
      debug_print("%4d %8.3lf", i, a[i]);
      if (i - 8*(i/8) == 0)
        debug_print("\n");
    }
    debug_print("\n");
#endif

    //calcA(xx, a, m, L);
    rmmult(xx, Ax, a, L + 1, m + 1, 1);

    if ( iter >= mit ) {
      debug_print("Maximum iterations reached\n");
      converged = false;
    }
  }
  for (i = 0; i < m; i++) {
    h[i] = a[m - i]/2;
    h[m + i + 1] = a[i + 1]/2;
  }
  h[m] = a[0]*r/2;

  return converged;
}


//===========================================================================================================
#ifndef CL2BP_STANDALONE

/* == Octave interface starts here ====================================== */
DEFUN_DLD (cl2bp, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {@var{h} =} cl2bp (@var{m}, @var{w1}, @var{w2}, @var{up}, @var{lo}\
 [, @var{gridsize}])\n\
\n\
Constrained L2 bandpass FIR filter design.  This is a fast implementation of the algorithm\n\
cited below.  Compared to @dfn{remez}, it offers implicit specification of transition bands,\n\
a higher likelihood of convergence, and an error criterion combining features of both L2 and\n\
Chebyshev approaches.@*\n\
Inputs:@*\n\
@var{m}: degree of cosine polynomial, i.e. the number of output coefficients will be @var{m}*2+1@*\n\
@var{w1}, @var{w2}: bandpass filter cutoffs in the range 0 <= @var{w1} < @var{w2} <= pi,\n\
where pi is the Nyquist frequency@*\n\
@var{up}: vector of 3 upper bounds for [stopband1, passband, stopband2]@*\n\
@var{lo}: vector of 3 lower bounds for [stopband1, passband, stopband2]@*\n\
@var{gridsize}: search grid size; larger values may improve accuracy,\n\
but greatly increase calculation time.@*\n\
Output:@*\n\
A vector of @var{m}*2+1 FIR coefficients, or an empty value if the solver failed to converge.@*\n\
Example:\n\
@example\n\
@code{h = cl2bp(30, 0.3*pi, 0.6*pi, [0.02, 1.02, 0.02], [-0.02, 0.98, -0.02], 2^11);}\n\
@end example\n\
Original Paper:  I. W. Selesnick, M. Lang, and C. S. Burrus.  A modified algorithm for\n\
constrained least square design of multiband FIR filters without specified transition bands.\n\
IEEE Trans. on Signal Processing, 46(2):497-501, February 1998.\n\
@end deftypefn\n\
@seealso{remez}\n")
{
  octave_value_list retval;
  int i;

  int nargin = args.length();
  if (nargin < 5 || nargin > 6) {
    print_usage ();
    return retval;
  }

  int m = args(0).int_value(true);
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (0));
    return retval;
  }
  if (m < 2 || m > 5000) {
    error("cl2bp: The m count must be between 2 and 5000");
    return retval;
  }

  double w1 = args(1).double_value();
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (1));
    return retval;
  }
  double w2 = args(2).double_value();
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (2));
    return retval;
  }
  if (w1 < 0 || w1 > w2 || w2 > 2*M_PI) {
    error("cl2bp: The cutoffs must be in the range 0 < w1 < w2 < 2*pi");
    return retval;
  }

  ColumnVector up_vector(args(3).vector_value());
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (3));
    return retval;
  }
  ColumnVector lo_vector(args(4).vector_value());
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (4));
    return retval;
  }
  if (up_vector.length() != 3 || lo_vector.length() != 3) {
    error("cl2bp: The up and lo vectors must contain 3 values");
    return retval;
  }

  double up[3];
  double lo[3];
  for (int i=0; i<3; ++i) {
    up[i] = up_vector(i);
    lo[i] = lo_vector(i);
  }

  int L = args(5).int_value(true);
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (5));
    return retval;
  }
  if (L > 1000000) {
    error("cl2bp: The \"gridsize\" parameter cannot exceed 1000000");
    return retval;
  }
  // Z is allocated as Z(2*L-1-2*m);
  if (L <= m) {
    error("cl2bp: The \"gridsize\" parameter must be larger than \"m\"");
    return retval;
  }

  ColumnVectorC h;
  cl2bp(h, m, w1, w2, up, lo, L, 1.e-5, 20);

  return octave_value(h);
}
#endif

/*
%!test
%! b = [
%!    0.0000000000000000
%!    0.0563980420304213
%!   -0.0000000000000000
%!   -0.0119990278695041
%!   -0.0000000000000001
%!   -0.3016146759510104
%!    0.0000000000000001
%!    0.5244313235801866
%!    0.0000000000000001
%!   -0.3016146759510104
%!   -0.0000000000000001
%!   -0.0119990278695041
%!   -0.0000000000000000
%!    0.0563980420304213
%!    0.0000000000000000];
%! assert(cl2bp(7, 0.25*pi, 0.75*pi, [0.01, 1.04, 0.01], [-0.01, 0.96, -0.01], 2^11), b, 1e-14);

*/

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Jaroslav Hajek-2]

On Fri, Oct 2, 2009 at 2:49 PM, Pete Gonzalez <pgonzalez@...> wrote:
> Hi Søren,
>

Hi,

> Attached is an updated version of cl2bp.cc with the following changes:
>
> - As requested, I added error_state checks for the input parameters. BTW I
> noticed that a lot of shipping Octave code is missing these checks; it
> should probably be audited at some point.
>
> - The MallocArray class you asked about was intended to support LGPL
> scenarios, i.e. where people can't use the GPL'ed oct.h.  That code is now
> isolated with a CL2BP_STANDALONE switch.  In a normal build, the native
> ColumnVector/Array classes are now used, but with C-style indexing
> operators.  (I didn't convert x[i] to x(i) everywhere because I think it
> would hinder readability.)
>

First the legal thing:
I'm afraid GPL can't be worked around as easily as you think (if it
could, dozens of companies would instantly do it). If you use Octave's
C++ API within your code, your code becomes a derived work of Octave
and hence its license needs to comply with GPL, use of LGPL is not
possible. It doesn't matter that the Octave-derived part can be
switched off by a directive. Even if it was commented out, GPL would
still apply.
In the present form, it can only be included in OctaveForge if you
accept changing the license to GPL.
If you really want a LGPL part, a possible way out is to isolate the
library from the Octave interface and distribute them separately
(In OctaveForge, they can be put together, but the license will become GPL).

Then the technical thing:
is rmmult really doing a matrix-matrix multiply? If yes, what are the
typical dimensions? For anything bigger than small (left vague), it
will probably be worth to defer this operation to BLAS to gain
performance.

hth

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Pete Gonzalez]

Jaroslav Hajek wrote:

> First the legal thing:
> I'm afraid GPL can't be worked around as easily as you think (if it
> could, dozens of companies would instantly do it). If you use Octave's
> C++ API within your code, your code becomes a derived work of Octave
> and hence its license needs to comply with GPL, use of LGPL is not
> possible. It doesn't matter that the Octave-derived part can be
> switched off by a directive. Even if it was commented out, GPL would
> still apply.
> In the present form, it can only be included in OctaveForge if you
> accept changing the license to GPL.
> If you really want a LGPL part, a possible way out is to isolate the
> library from the Octave interface and distribute them separately
> (In OctaveForge, they can be put together, but the license will become GPL).

What you are saying is that the files distributed with Octave become
GPL, regardless of their unmodified copyright banners displaying other
licenses such as LGPL, BSD, public domain, etc.  Having chosen LGPL, we
are of course implicitly granting Octave this right.

However, as the copyright holders, we can separately contribute our
original code (without GPL modifications) to other open source projects
under LGPL, right?  The LGPL banner/switch was intended to advertise
that this other licensing option exists, and to make it easier to track
changes between the two (e.g. if someone wanted to contribute the same
patch to both branches).

I can see that the CL2BP_STANDALONE directive may be intertwining those
branches a little too intimately.  If you prefer, I can separate
cl2bp.cc into two files -- one with an LGPL license doesn't reference
Octave at all, and one with an explicit GPL license that contains the
Octave wrapper.  Would that satisfy your concerns?

> Then the technical thing:
> is rmmult really doing a matrix-matrix multiply? If yes, what are the
> typical dimensions? For anything bigger than small (left vague), it
> will probably be worth to defer this operation to BLAS to gain
> performance.

You are right that an optimized matrix library might improve
performance.  However, the library would need to be compatible with the
LGPL license, and it would need to be relatively small and portable.  (I
think BLAS is a little "baroque" for our needs.)  Also, this change is
nontrivial programming work that we probably won't be able to do right
away.  But I'm open to suggestions.

Thanks for your comments.

Cheers,
-Pete

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Jaroslav Hajek-2]

On Fri, Oct 2, 2009 at 8:18 PM, Pete Gonzalez <pgonzalez@...> wrote:

>
>
> Jaroslav Hajek wrote:
>>
>> First the legal thing:
>> I'm afraid GPL can't be worked around as easily as you think (if it
>> could, dozens of companies would instantly do it). If you use Octave's
>> C++ API within your code, your code becomes a derived work of Octave
>> and hence its license needs to comply with GPL, use of LGPL is not
>> possible. It doesn't matter that the Octave-derived part can be
>> switched off by a directive. Even if it was commented out, GPL would
>> still apply.
>> In the present form, it can only be included in OctaveForge if you
>> accept changing the license to GPL.
>> If you really want a LGPL part, a possible way out is to isolate the
>> library from the Octave interface and distribute them separately
>> (In OctaveForge, they can be put together, but the license will become
>> GPL).
>
> What you are saying is that the files distributed with Octave become GPL,
> regardless of their unmodified copyright banners displaying other licenses
> such as LGPL, BSD, public domain, etc.  Having chosen LGPL, we are of course
> implicitly granting Octave this right.
>

It's not "of course" - you, as the author, choose licenses for your
work. You *may* choose not to comply with GPL, but then you simply can
not distribute the source.

> However, as the copyright holders, we can separately contribute our original
> code (without GPL modifications) to other open source projects under LGPL,
> right?  The LGPL banner/switch was intended to advertise that this other
> licensing option exists, and to make it easier to track changes between the
> two (e.g. if someone wanted to contribute the same patch to both branches).
>

No, it doesn't work this way. Once you create a derivative work of
Octave (which your source is) and release it under GPL (which you need
if you want to distribute it), you can't provide an option to
"downgrade" the license. You can, of course, take the source you
posted (because it is your work), strip away the Octave stuff, and
then license it under LGPL. Users, however, cannot "unGPL" the code -
only you can. This is an example of the "viral" nature of GPL that is
often referred to.

> I can see that the CL2BP_STANDALONE directive may be intertwining those
> branches a little too intimately.  If you prefer, I can separate cl2bp.cc
> into two files -- one with an LGPL license doesn't reference Octave at all,
> and one with an explicit GPL license that contains the Octave wrapper.
>  Would that satisfy your concerns?

I think that would be OK. However, if it becomes part of OctaveForge,
it will be covered by GPL. Anyone getting the file from OctaveForge
will be bound by GPL. Anyone getting the same source from you directly
will be allowed to use LGPL.

>> Then the technical thing:
>> is rmmult really doing a matrix-matrix multiply? If yes, what are the
>> typical dimensions? For anything bigger than small (left vague), it
>> will probably be worth to defer this operation to BLAS to gain
>> performance.
>
> You are right that an optimized matrix library might improve performance.
>  However, the library would need to be compatible with the LGPL license, and
> it would need to be relatively small and portable.  (I think BLAS is a
> little "baroque" for our needs.)  Also, this change is nontrivial
> programming work that we probably won't be able to do right away.  But I'm
> open to suggestions.
>

As I said, it depends on the dimensions. If the dimensions are, say,
10x10, you'll gain little or nothing. If it's 100x100 or 1000x1000,
you'll gain a lot. BLAS is today more of a standard interface than a
library, so there's no licensing problem - you can use anything. If
you use Octave's classes - say, Matrix and ColumnVector in your file,
then Matrix*ColumnVector operation is already linked to BLAS.

hth

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Pete Gonzalez]

Jaroslav Hajek wrote:
>> What you are saying is that the files distributed with Octave become GPL,
>> regardless of their unmodified copyright banners displaying other licenses
>> such as LGPL, BSD, public domain, etc.  Having chosen LGPL, we are of course
>> implicitly granting Octave this right.
>
> It's not "of course" - you, as the author, choose licenses for your
> work. You *may* choose not to comply with GPL, but then you simply can
> not distribute the source.

My point was that if we choose to reject GPL, then we must also reject
LGPL.  By choosing LGPL, we implicitly are giving the Octave maintainers
(and anyone else) the option to repackage our code as GPL.

>> I can see that the CL2BP_STANDALONE directive may be intertwining those
>> branches a little too intimately.  If you prefer, I can separate cl2bp.cc
>> into two files -- one with an LGPL license doesn't reference Octave at all,
>> and one with an explicit GPL license that contains the Octave wrapper.
>>  Would that satisfy your concerns?
>
> I think that would be OK. However, if it becomes part of OctaveForge,
> it will be covered by GPL. Anyone getting the file from OctaveForge
> will be bound by GPL. Anyone getting the same source from you directly
> will be allowed to use LGPL.

In other words, the GPL and LGPL branches can evolve in parallel.
Patches cannot be copied verbatim from the GPL branch into the LGPL
branch, but the author of a patch is allowed to apply their patch to
both branches.  And merging from the LGPL branch to the GPL branch is
always allowed.

But this assumes that the two branches don't diverge too much, which I
think is a valid reason for preserving the LGPL form in Octave, even
though the license has effectively become GPL.

Attached are three new files.  The files cl2bp_lib.cc and .h are based
on our original prototype before Octave integration.  The license banner
is LGPL, and no Octave header files are used.  (In particular, it uses
the original MallocArray class instead of Octave's ColumnVector.)  The
Octave wrapper is in cl2bp.cc, which has the full GPL license banner.

Let me know if this looks good to you.

>> You are right that an optimized matrix library might improve performance.
>>  However, the library would need to be compatible with the LGPL license, and
>> it would need to be relatively small and portable.  (I think BLAS is a
>> little "baroque" for our needs.)  Also, this change is nontrivial
>> programming work that we probably won't be able to do right away.  But I'm
>> open to suggestions.
>
> As I said, it depends on the dimensions. If the dimensions are, say,
> 10x10, you'll gain little or nothing. If it's 100x100 or 1000x1000,
> you'll gain a lot. BLAS is today more of a standard interface than a
> library, so there's no licensing problem - you can use anything. If
> you use Octave's classes - say, Matrix and ColumnVector in your file,
> then Matrix*ColumnVector operation is already linked to BLAS.

Yes, but Matrix and ColumnVector are GPL, not LGPL.  ;-)

Cheers,
-Pete

/*

Copyright (c) 2008-2009, Evgeni A. Nurminski
Copyright (c) 2008-2009, Pete Gonzalez

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 software; see the file COPYING.  If not, see
<http://www.gnu.org/licenses/>.


Authors:

Evgeni A. Nurminski <nurmi@...>
Institute for Automation and Control Problems
Far East branch of RAS, Vladivostok, Russia

Pete Gonzalez <pgonzalez@...>
Bluel Technologies Corporation

*/

#ifndef CL2BP_STANDALONE  // Define this symbol for usage outside of Octave
#include <octave/oct.h>
#endif

#include "cl2bp_lib.h"

static void cancel_handler(void *) {
  OCTAVE_QUIT;
}

// When CL2BP_LOGGING is defined, this will direct messages to the Octave pager
void cl2bp_log(const char *message) {
  octave_stdout << message;
}


DEFUN_DLD (cl2bp, args, ,
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {@var{h} =} cl2bp (@var{m}, @var{w1}, @var{w2}, @var{up}, @var{lo}\
 [, @var{gridsize}])\n\
\n\
Constrained L2 bandpass FIR filter design.  This is a fast implementation of the algorithm\n\
cited below.  Compared to @dfn{remez}, it offers implicit specification of transition bands,\n\
a higher likelihood of convergence, and an error criterion combining features of both L2 and\n\
Chebyshev approaches.@*\n\
Inputs:@*\n\
@var{m}: degree of cosine polynomial, i.e. the number of output coefficients will be @var{m}*2+1@*\n\
@var{w1}, @var{w2}: bandpass filter cutoffs in the range 0 <= @var{w1} < @var{w2} <= pi,\n\
where pi is the Nyquist frequency@*\n\
@var{up}: vector of 3 upper bounds for [stopband1, passband, stopband2]@*\n\
@var{lo}: vector of 3 lower bounds for [stopband1, passband, stopband2]@*\n\
@var{gridsize}: search grid size; larger values may improve accuracy,\n\
but greatly increase calculation time.@*\n\
Output:@*\n\
A vector of @var{m}*2+1 FIR coefficients, or an empty value if the solver failed to converge.@*\n\
Example:\n\
@example\n\
@code{h = cl2bp(30, 0.3*pi, 0.6*pi, [0.02, 1.02, 0.02], [-0.02, 0.98, -0.02], 2^11);}\n\
@end example\n\
Original Paper:  I. W. Selesnick, M. Lang, and C. S. Burrus.  A modified algorithm for\n\
constrained least square design of multiband FIR filters without specified transition bands.\n\
IEEE Trans. on Signal Processing, 46(2):497-501, February 1998.\n\
@end deftypefn\n\
@seealso{remez}\n")
{
  octave_value_list retval;
  int i;

  int nargin = args.length();
  if (nargin < 5 || nargin > 6) {
    print_usage ();
    return retval;
  }

  int m = args(0).int_value(true);
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (0));
    return retval;
  }
  if (m < 2 || m > 5000) {
    error("cl2bp: The m count must be between 2 and 5000");
    return retval;
  }

  double w1 = args(1).double_value();
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (1));
    return retval;
  }
  double w2 = args(2).double_value();
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (2));
    return retval;
  }
  if (w1 < 0 || w1 > w2 || w2 > 2*M_PI) {
    error("cl2bp: The cutoffs must be in the range 0 < w1 < w2 < 2*pi");
    return retval;
  }

  ColumnVector up_vector(args(3).vector_value());
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (3));
    return retval;
  }
  ColumnVector lo_vector(args(4).vector_value());
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (4));
    return retval;
  }
  if (up_vector.length() != 3 || lo_vector.length() != 3) {
    error("cl2bp: The up and lo vectors must contain 3 values");
    return retval;
  }

  double up[3];
  double lo[3];
  for (int i=0; i<3; ++i) {
    up[i] = up_vector(i);
    lo[i] = lo_vector(i);
  }

  int L = args(5).int_value(true);
  if (error_state) {
    gripe_wrong_type_arg ("cl2bp", args (5));
    return retval;
  }
  if (L > 1000000) {
    error("cl2bp: The \"gridsize\" parameter cannot exceed 1000000");
    return retval;
  }
  // Z is allocated as Z(2*L-1-2*m);
  if (L <= m) {
    error("cl2bp: The \"gridsize\" parameter must be larger than \"m\"");
    return retval;
  }

  MallocArray<double> h;
  cl2bp(h, m, w1, w2, up, lo, L, 1.e-5, 20, cancel_handler);

  ColumnVector h_vector(h.get_length());

  for (int i=0; i<h.get_length(); ++i)
    h_vector(i) = h[i];

  return octave_value(h_vector);
}

/*
%!test
%! b = [
%!    0.0000000000000000
%!    0.0563980420304213
%!   -0.0000000000000000
%!   -0.0119990278695041
%!   -0.0000000000000001
%!   -0.3016146759510104
%!    0.0000000000000001
%!    0.5244313235801866
%!    0.0000000000000001
%!   -0.3016146759510104
%!   -0.0000000000000001
%!   -0.0119990278695041
%!   -0.0000000000000000
%!    0.0563980420304213
%!    0.0000000000000000];
%! assert(cl2bp(7, 0.25*pi, 0.75*pi, [0.01, 1.04, 0.01], [-0.01, 0.96, -0.01], 2^11), b, 1e-14);

*/

/*

Copyright (c) 2008-2009, Evgeni A. Nurminski
Copyright (c) 2008-2009, Pete Gonzalez

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.


Authors:

Evgeni A. Nurminski <nurmi@...>
Institute for Automation and Control Problems
Far East branch of RAS, Vladivostok, Russia

Pete Gonzalez <pgonzalez@...>
Bluel Technologies Corporation

*/

#ifdef _MSC_VER
#define _CRT_SECURE_NO_WARNINGS  // disable spurious warnings
#define _USE_MATH_DEFINES // for math.h
#endif

#include "cl2bp_lib.h"

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <math.h>
#include <string.h>

#define BIG_NUMBER 100000

//-----------------------------------------------------------------------------------------------------------
#ifdef CL2BP_LOGGING
static void log_printf(const char *format, ...) {
  va_list argptr;
  va_start(argptr, format);

  char buf[20*1024];
  if (vsnprintf(buf,20*1024-1,format,argptr) == -1) {
    strcpy(buf, "#ERROR#");
  }
  va_end(argptr);

  cl2bp_log(buf);
}
#else
#define log_printf(...)  ((void)0)   // don't log anything (for a big speed improvement)
#endif

//-----------------------------------------------------------------------------------------------------------
static int local_max(const MallocArray<double>& x, int n, MallocArray<int>& ix) {
  int i, mx;

  mx = 0;

  MallocArray<double> xx(n+2);

  xx[0] = xx[n + 1] = -BIG_NUMBER;
  for ( i = 1; i <= n; i++)
    xx[i] = x[i - 1];
  for ( i = 0; i < n; i++ ) {
    (((xx[i] <  xx[i + 1]) && (xx[i + 1] >  xx[i + 2])) ||
     ((xx[i] <  xx[i + 1]) && (xx[i + 1] >= xx[i + 2])) ||
     (xx[i] <= xx[i + 1]) && (xx[i + 1] >  xx[i + 2]) )
    && ( ix[mx++] = i );
  }
  return mx;
}

//-----------------------------------------------------------------------------------------------------------
static int solvps(MallocArray<double>& a, MallocArray<double>& b, int n,
  void (*cancel_handler)(void *), void *cancel_state) {

  double t;
  int j, k;
  int a_p;
  int a_q;
  int a_r;
  int a_s;

  // In empirical tests, 2% to 6% of the execution time was spent in solvps()
  if (cancel_handler) cancel_handler(cancel_state);

  for (j = 0, a_p = 0; j < n; ++j, a_p += n+1) {
    for (a_q = j*n; a_q < a_p; ++a_q)
      a[a_p] -= a[a_q] * a[a_q];
    if (a[a_p] <= 0.)
      return -1;
    a[a_p] = sqrt(a[a_p]);
    for (k = j + 1, a_q = a_p + n; k < n; ++k, a_q += n) {
      for (a_r = j*n, a_s = k*n, t = 0.; a_r < a_p;)
        t += a[a_r++] * a[a_s++];
      a[a_q] -= t;
      a[a_q] /= a[a_p];
    }
  }
  for (j = 0, a_p = 0; j < n; ++j, a_p += n+1) {
    for (k = 0, a_q = j*n; k < j;)
      b[j] -=b [k++] * a[a_q++];
    b[j] /= a[a_p];
  }
  for (j = n - 1, a_p = n*n - 1; j >= 0; --j, a_p -= n + 1) {
    for (k = j + 1, a_q = a_p + n; k < n; a_q += n)
      b[j] -= b[k++]* a[a_q];
    b[j] /= a[a_p];
  }
  return 0;
}

//-----------------------------------------------------------------------------------------------------------
static void rmmult(MallocArray<double>& rm, const MallocArray<double>& a, const MallocArray<double>& b,
  int n, int m, int l,
  void (*cancel_handler)(void *), void *cancel_state) {

  double z;
  int i,j,k;
  int b_p; // index into b
  int a_p; // index into a
  int rm_start = 0;  // index into rm
  int rm_q; // index into rm

  MallocArray<double> q0(m);
  for (i = 0; i < l ; ++i, ++rm_start) {
    // In empirical tests, 87% to 95% of the execution time was spent in rmmult()
    if (cancel_handler) cancel_handler(cancel_state);
    for (k = 0, b_p = i; k < m; b_p += l)
      q0[k++] = b[b_p];
    for (j = 0, a_p = 0, rm_q = rm_start; j < n; ++j, rm_q += l) {
      for (k = 0, z = 0.; k < m;)
        z += a[a_p++] * q0[k++];
      rm[rm_q]=z;
    }
  }
}

//-----------------------------------------------------------------------------------------------------------
static void mattr(MallocArray<double>& a, const MallocArray<double>& b, int m, int n) {
  int i, j;
  int b_p;
  int a_p = 0;
  int b_start = 0;
  for (i = 0; i < n; ++i, ++b_start)
    for ( j = 0, b_p = b_start; j < m; ++j, b_p += n )
      a[a_p++] = b[b_p];
}

//-----------------------------------------------------------------------------------------------------------
static void calcAx(MallocArray<double>& Ax, int m, int L) {
  double r = M_SQRT2, pi = M_PI;
  int i, j;

  Ax.resize((L+1)*(m + 1));

  for ( i = 0; i <= L; i++)
    for ( j = 0; j <= m; j++)
      Ax[i*(m+1) + j] = cos(i*j*pi/L);
  for ( i = 0; i < (L+1)*(m+1); i += m + 1 )
    Ax[i] /= r;
}

//-----------------------------------------------------------------------------------------------------------
static double L2error(const MallocArray<double>& x, const MallocArray<double>& w, int L, double w1, double w2) {
  double xx, ww, sumerr, pi = M_PI;
  int i;
  for ( i = 0, sumerr = 0; i < L + 1; i++ ) {
    ww = w[i];
    xx = x[i];
    sumerr += ( ww < w1*pi || ww > w2*pi ) ?  xx*xx : (1 - xx)*(1 - xx);
  }
  return sumerr;
}

//-----------------------------------------------------------------------------------------------------------
static double CHerror(double *wmin, const MallocArray<double>& x, const MallocArray<double>& w,
  int L, double w1, double w2) {

  double xx, ww, err, errmax;
  int i;
  errmax = -1;
  *wmin = 0;
  for ( i = 0; i < L + 1; i++ ) {
    ww = w[i];
    xx = x[i];
    err = (( ww < w1 ) || (ww > w2 )) ?  fabs(xx) : fabs(1 - xx);
    if ( err > errmax ) {
      errmax = err;
      *wmin = ww;
    }
  }
  return errmax;
}

//-----------------------------------------------------------------------------------------------------------
static void Ggen(MallocArray<double>& G, int m, const MallocArray<double>& w, const MallocArray<int>& kmax,
  int l_kmax, const MallocArray<int>& kmin, int l_kmin) {

  G.resize((l_kmax + l_kmin)*(m + 1));

  int nsys, i, j;
  double r = M_SQRT2;

  nsys = l_kmax + l_kmin;
  for ( i = 0; i < l_kmax; i++)
    for ( j = 0; j <= m; j++)
      G[i*(m+1) + j] = cos(w[kmax[i]]*j);
  for ( i = l_kmax; i < nsys; i++)
    for ( j = 0; j <= m; j++)
      G[i*(m+1) + j] = -cos(w[kmin[i - l_kmax]]*j);
  for ( i = 0; i < nsys*(m+1); i += m+1 )
    G[i] /= r;
}

//-----------------------------------------------------------------------------------------------------------
bool cl2bp(MallocArray<double>& h, int m, double w1, double w2, double up[3], double lo[3], int L,
  double eps, int mit, void (*cancel_handler)(void *), void *cancel_state) {

  double r = M_SQRT2;
  double pi = M_PI;
  double wmin, ww, xw;
  int q1, q2, q3, i, iter = 0;
  int off, neg;

  int ik;
  int l_kmax;
  int l_kmin;
  int l_okmax;
  int l_okmin;
  double uvo = 0, lvo = 0;

  double err, diff_up, diff_lo;
  double erru, erro;
  int iup, ilo;

  int nsys, j;

  int imin;
  double umin;

  int k1, k2, ak1, ak2;
  double cerr;

  h.resize(2*m+1);

  bool converged = true;

  MallocArray<double> x0(L+1);
  MallocArray<double> x1(L+1);
  MallocArray<double> xx(L+1);
  MallocArray<double> xm1(L+1);
  MallocArray<double> w(L+1);
  MallocArray<double> u(L+1);
  MallocArray<double> l(L+1);
  MallocArray<double> a(L+1);
  MallocArray<double> c(L+1);
  MallocArray<int> kmax(L+1);
  MallocArray<int> kmin(L+1);
  l_kmax  = l_kmin = 0;

  MallocArray<int> okmin(L+1);
  MallocArray<int> okmax(L+1);
  l_okmax = l_okmin = 0;
  MallocArray<double> rhs(m+2);
  MallocArray<double> mu(2*(L+1));

  for ( i = 0; i <= L; i++ )
    w[i] = i*pi/L;

  MallocArray<double> Z(2*L-1-2*m);

  q1 = (int)floor(L*w1/pi);
  q2 = (int)floor(L*(w2-w1)/pi);
  q3 = L + 1 - q1 - q2;

  off = 0;
  for ( i = 0; i < q1; i++) {
    u[i] = up[0];
    l[i] = lo[0];
  }
  off += i;
  for ( i = 0; i < q2; i++) {
    u[off + i] = up[1];
    l[off + i] = lo[1];
  }
  off += i;
  for ( i = 0; i < q3; i++) {
    u[off + i] = up[2];
    l[off + i] = lo[2];
  }


  c[0] = (w2-w1)*r;
  for ( i = 1; i <= m; i++)
    c[i] =  2*(sin(w2*i)-sin(w1*i))/i;
  for ( i = 0; i <= m; i++) {
    c[i] /=  pi;
    a[i] = c[i];
  }
  log_printf("\nInitial approximation. Unconstrained L2 filter coefficients.\n");
  log_printf("=============================================================\n");

  log_printf("\nZero order term %8.3lf\n", a[0]);
  for ( i = 1; i <= m; i++) {
    log_printf("%4d %8.3lf", i, a[i]);
    if (i - 8*(i/8) == 0)
      log_printf("\n");
  }
  log_printf("\n");

  // Precalculate Ax
  MallocArray<double> Ax;
  calcAx(Ax, m, L);

  //calcA(x0, a, m, L);
  rmmult(x0, Ax, a, L + 1, m + 1, 1, cancel_handler, cancel_state);

  err = CHerror(&wmin, x0, w, L, w1, w2);
  log_printf("\nChebyshev err %12.4e (%11.5e)  <> L2 err %12.4e\n", err, wmin/pi, L2error(x0, w, L, w1, w2)/(L+1));
  for (iter = 1; ; iter++) {
    log_printf("\n:::::: Iteration %3d :::::::::::\n", iter);

    if ( (uvo > eps*2) || (lvo > eps*2) ) {
      log_printf("\nXXX Take care of old errors: uvo lvo %12.3e %12.3e", uvo, lvo);
      if( k1 >= 0 ) log_printf(" k1 %3d(%d)", k1, okmax[k1]);
      if( k2 >= 0 ) log_printf(" k2 %3d(%d)", k2, okmin[k2]);
      log_printf("\n");

      if ( uvo > lvo ) {
        kmax[l_kmax] = okmax[k1];
        l_kmax++;
        l_okmax--;
        for (i = k1; i < l_okmax; i++)
          okmax[i] = okmax[i + 1];
      } else {
        kmin[l_kmin] = okmin[k2];
        l_kmin++;
        l_okmin--;
        for (i = k2; i < l_okmin; i++)
          okmin[i] = okmin[i + 1];
      }
      nsys = l_kmax + l_kmin;

      /*
        for (i = 0; i < l_kmax; i++)
          log_printf("i %3d kmax %3d mu %12.4e\n",
             i, kmax[i], mu[i]);
        log_printf("\n");
        for (i = 0; i < l_kmin; i++)
          log_printf("i %3d kmin %3d mu %12.4e\n",
             i, kmin[i], mu[i + l_kmax]);
        log_printf("\n\n");
      */
    } else {

      //calcA(x1, a, m, L);
      rmmult(x1, Ax, a, L + 1, m + 1, 1, cancel_handler, cancel_state);


      for ( i = 0; i < l_kmax; i++)
        okmax[i] = kmax[i];
      for ( i = 0; i < l_kmin; i++)
        okmin[i] = kmin[i];
      l_okmax = l_kmax;
      l_okmin = l_kmin;

      l_kmax = local_max(x1, L + 1, kmax);


      for ( i = 0; i < L + 1; i++)
        xm1[i] = -x1[i];
      l_kmin = local_max(xm1, L + 1, kmin);

      log_printf("\nSignificant deviations from the ideal filter. Levels:");
      log_printf(" %8.2e %8.2e %8.2e (lo)", lo[0], lo[1], lo[2]);
      log_printf(" %8.2e %8.2e %8.2e (up)", up[0], up[1], up[2]);
      log_printf("\n");

      for ( i = 0, ik = 0; i < l_kmax; i++) {
        j = kmax[i];
        if ( x1[j] > u[j] - 10*eps ) {
          kmax[ik] = j;
          ik++;
        }
      }
      l_kmax = ik;

      log_printf("overshoots = %d\n", l_kmax);
      for ( i = 0; i < l_kmax; i++) {
        j = kmax[i];
        ww = w[j];
        xw = x1[j];
        err = (w1 < ww && ww < w2) ? xw - 1 : xw;
        log_printf(" i = %3d kmax = %3d xw = %12.5e err = %10.3e u = %12.4e max = %9.2e\n",
               i, j, xw, err, u[j], xw - u[j]);
      }

      for ( i = 0, ik = 0; i < l_kmin; i++) {
        j = kmin[i];
        if ( x1[j] < l[j] + 10*eps ) {
          kmin[ik] = j;
          ik++;
        }
      }
      l_kmin = ik;

      log_printf("undershoots = %d\n", l_kmin);
      for ( i = 0; i < l_kmin; i++) {
        j = kmin[i];
        ww = w[j];
        xw = -xm1[j];
        err =(w1 < ww && ww < w2) ? xw - 1 : xw;
        log_printf(" i = %3d kmin = %3d xw = %12.5e err = %10.3e l = %12.4e min = %9.2e\n",
               i, j, xw, err, l[j], xw - l[j]);
      }

      err = erru = erro = diff_up = diff_lo = 0;
      iup = ilo = 0;
      for ( i = 0; i < l_kmax; i++) {
        ik = kmax[i];
        diff_up = x1[ik] - u[ik];
        if ( diff_up > erru ) {
          erru = diff_up;
          iup = i;
        }
      }
      for ( i = 0; i < l_kmin; i++) {
        ik = kmin[i];
        diff_lo = l[ik] - x1[ik];
        if ( diff_lo > erro ) {
          erro = diff_lo;
          ilo = i;
        }
      }
      err = erro > erru ? erro : erru;
      log_printf("erru = %14.6e iup = %4d kmax = %4d ", erru, iup, kmax[iup]);
      log_printf("erro = %14.6e ilo = %4d kmin = %4d\n", erro, ilo, kmin[ilo]);

      if ( err < eps )
        break;
    }


    nsys = l_kmax + l_kmin;

    MallocArray<double> G(nsys*(m + 1));
    MallocArray<double> GT(nsys*(m + 1));
    MallocArray<double> GG(nsys*nsys);

    for ( i = 0; i < l_kmax; i++)
      for ( j = 0; j <= m; j++)
        G[i*(m+1) + j] = cos(w[kmax[i]]*j);
    for ( i = l_kmax; i < nsys; i++)
      for ( j = 0; j <= m; j++)
        G[i*(m+1) + j] = -cos(w[kmin[i - l_kmax]]*j);
    for ( i = 0; i < nsys*(m+1); i += m+1 )
      G[i] /= r;
    mattr(GT, G, nsys, m + 1);
    rmmult(GG, G, GT, nsys, m + 1, nsys, cancel_handler, cancel_state);


    rmmult(rhs, G, c, nsys, m + 1, 1, cancel_handler, cancel_state);
    for ( i = 0; i < nsys; i++)
      if ( i < l_kmax )
        rhs[i] -= u[kmax[i]];
      else
        rhs[i] += l[kmin[i - l_kmax]];

    for ( i = 0; i < nsys; ++i)
      mu[i] = rhs[i];

    solvps(GG, mu, nsys, cancel_handler, cancel_state);
    log_printf("\nXXX KT-system with %d equations resolved.\n", nsys);
    for ( i = 0, neg = 0; i < nsys; i++)
      if ( mu[i] < 0 ) neg++;
    log_printf("\nTotal number of negative multipliers = %3d\n\n", neg);
    while (neg) {


      umin = BIG_NUMBER;
      for ( i = 0, j = 0; i < nsys; i++) {
        if ( mu[i] >= 0 ) continue;
        j++;
        if ( mu[i] < umin ) {
          imin = i;
          umin = mu[i];
        }
      }

      if ( umin > 0 )
        break;
      log_printf(" neg = %3d of %3d imin = %3d mu-min = %12.4e\n", j, nsys, imin, umin);

      for ( i = imin; i < nsys - 1; i++) {
        rhs[i] = rhs[i + 1];
        for ( j = 0; j <= m; j++)
          G[i*(m + 1) + j] =  G[(i + 1)*(m + 1) + j];
      }

      if ( imin < l_kmax ) {
        for ( i = imin; i < l_kmax - 1; i++)
          kmax[i] = kmax[i + 1];
        l_kmax--;
      } else {
        for ( i = imin; i < nsys - 1; i++)
          kmin[i - l_kmax] = kmin[i - l_kmax + 1];
        l_kmin--;
      }

      --nsys;

      mattr(GT, G, nsys, m + 1);
      rmmult(GG, G, GT, nsys, m + 1, nsys, cancel_handler, cancel_state);
      for ( i = 0; i < nsys; ++i)
        mu[i] = rhs[i];
      solvps(GG, mu, nsys, cancel_handler, cancel_state);
    }

    MallocArray<double> da(m + 1);
    MallocArray<double> zd(nsys);

    rmmult(da, GT, mu, m + 1, nsys, 1, cancel_handler, cancel_state);
    for ( i = 0; i <= m; i++)
      a[i] = c[i] - da[i];
    rmmult(zd, G, a, nsys, m + 1, 1, cancel_handler, cancel_state);

    MallocArray<double> zz(l_okmax + l_okmin);
    Ggen(G, m, w, okmax, l_okmax, okmin, l_okmin);
    rmmult(zz, G, a, l_okmax + l_okmin, m + 1, 1, cancel_handler, cancel_state);
    uvo = lvo = eps;
    k1 = k2 = -1;
    ak1 = ak2 = -1;
    if (l_okmax + l_okmin > 0)
      log_printf("\nErrors on the previous set of freqs\n\n");
    for (i = 0; i < l_okmax; i++) {
      j = okmax[i];
      cerr = zz[i] - u[j];
      log_printf(" i %2d j %3d u %12.4e Ga %12.4e cerr %12.4e\n",
             i, j, u[j], zz[i], cerr);
      if ( cerr > uvo ) {
        uvo = cerr;
        k1 = i;
        ak1 = j;
      }
    }
    cerr = 0;
    for (i = 0; i < l_okmin; i++) {
      j = okmin[i];
      cerr = l[j] + zz[i + l_okmax];
      log_printf(" i %2d j %3d l %12.4e Ga %12.4e cerr %12.4e\n",
             i, j, l[j], zz[i + l_okmax], cerr);
      if ( cerr > lvo ) {
        lvo = cerr;
        k2 = i, ak2 = j;
      }
    }
    if ( l_okmax + l_okmin > 0 ) {
      log_printf("\n uvo = %12.4e k1 = %4d (%3d) ", uvo, k1, ak1);
      log_printf("  lvo = %12.4e k2 = %4d (%3d) ", lvo, k2, ak2);
      log_printf(" maxerr = %12.4e\n", uvo > lvo ? uvo : lvo);
    }

    log_printf("\nConstrained L2 band filter coefficients.\n");
    log_printf("=====================================\n");

#ifdef CL2BP_LOGGING
    log_printf("\nZero order term %8.3lf\n", a[0]);
    for ( i = 1; i <= m; i++) {
      log_printf("%4d %8.3lf", i, a[i]);
      if (i - 8*(i/8) == 0)
        log_printf("\n");
    }
    log_printf("\n");
#endif

    //calcA(xx, a, m, L);
    rmmult(xx, Ax, a, L + 1, m + 1, 1, cancel_handler, cancel_state);

    if ( iter >= mit ) {
      log_printf("Maximum iterations reached\n");
      converged = false;
    }
  }
  for (i = 0; i < m; i++) {
    h[i] = a[m - i]/2;
    h[m + i + 1] = a[i + 1]/2;
  }
  h[m] = a[0]*r/2;

  return converged;
}

/*

Copyright (c) 2008-2009, Evgeni A. Nurminski
Copyright (c) 2008-2009, Pete Gonzalez

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.


Authors:

Evgeni A. Nurminski <nurmi@...>
Institute for Automation and Control Problems
Far East branch of RAS, Vladivostok, Russia

Pete Gonzalez <pgonzalez@...>
Bluel Technologies Corporation

*/

#ifndef CL2BP_H
#define CL2BP_H

#include <assert.h>

//-----------------------------------------------------------------------------------------------------------
// If you want to debug the cl2bp algorithm, define the CL2BP_LOGGING symbol and provide an
// implementation of cl2bp_log().
#ifdef CL2BP_LOGGING
extern void cl2bp_log(const char *message);
#endif

//-----------------------------------------------------------------------------------------------------------
// This is a simple helper class that performs bounds-checking for debug builds, but reduces to an unchecked
// malloc() array for release builds.
template <class T>
class MallocArray {
  int length;
  T *ptr;
public:
  T& operator[](int index) {
    assert(index >= 0 && index < length);
    return ptr[index];
  }
  T operator[](int index) const {
    assert(index >= 0 && index < length);
    return ptr[index];
  }

  int get_length() const { return length; }

  void resize(int length_) {
    assert(length_ >= 0 && length_ <= 512*1024*1024);  // verify that the array size is reasonable
    length = length_;
    ptr = (T *)realloc(ptr, length * sizeof(T));
    memset(ptr, 0, length * sizeof(T));
  }

  MallocArray(int length_=0) {
    ptr = 0;
    length = 0;
    if (length_) resize(length_);
  }
  ~MallocArray() { free(ptr); }
private:
  MallocArray(const MallocArray& src) { }  // copy constructor is unimplemented and disallowed
};

//-----------------------------------------------------------------------------------------------------------
// Constrained L2 BandPass FIR filter design
//  h       2*m+1 filter coefficients (output)
//  m       degree of cosine polynomial
//  w1,w2   fist and second band edges
//  up      upper bounds
//  lo      lower bounds
//  L       grid size
//  eps     stopping condition ("SN")
//  mit     maximum number of iterations
//
// cl2bp() returns true if the solver converged, or false if the maximum number of iterations was reached.
// If provided, the cancel_handler function pointer will be called periodically inside long-running loops
// giving an opportunity to abort.  The cancel_state is a user-defined object, e.g. a pointer to a cancel
// button or other object.
//
// Example usage:
//   MallocArray<double> coefficients;
//   double up[3] = { 0.02, 1.02, 0.02 };
//   double lo[3] = { -0.02, 0.98, -0.02 };
//   cl2bp(coefficients, 30, 0.3*M_PI, 0.6*M_PI, up, lo, 1<<11, 1.e-5, 20);
//
// The algorithm is described in this paper:
//   I. W. Selesnick, M. Lang, and C. S. Burrus.  A modified algorithm for constrained least square
//   design of multiband FIR filters without specified transition bands.  IEEE Trans. on Signal
//   Processing, 46(2):497-501, February 1998.
bool cl2bp(MallocArray<double>& h, int m, double w1, double w2, double up[3], double lo[3], int L,
  double eps, int mit, void (*cancel_handler)(void *)=0, void *cancel_state=0);

#endif

------------------------------------------------------------------------------
Come build with us! The BlackBerry® Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9-12, 2009. Register now!
http://p.sf.net/sfu/devconf
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Søren Hauberg]

søn, 04 10 2009 kl. 05:18 -0400, skrev Pete Gonzalez:
> Attached are three new files.  The files cl2bp_lib.cc and .h are based
> on our original prototype before Octave integration.  The license banner
> is LGPL, and no Octave header files are used.  (In particular, it uses
> the original MallocArray class instead of Octave's ColumnVector.)  The
> Octave wrapper is in cl2bp.cc, which has the full GPL license banner.
>
> Let me know if this looks good to you.

It seems nobody objects to your suggested code, so I guess it is fine
for inclusion. Could you send a patch that integrates the code with the
package, i.e. it should also make needed changes to 'Makefile' and
'INDEX' ?

Søren


------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Pete Gonzalez]

Søren Hauberg wrote:
 > It seems nobody objects to your suggested code, so I guess it is fine
 > for inclusion. Could you send a patch that integrates the code with the
 > package, i.e. it should also make needed changes to 'Makefile' and
 > 'INDEX' ?

Attached is a patch against the svn HEAD revision for this directory:

https://octave.svn.sourceforge.net/svnroot/octave/trunk/octave-forge/main/signal

I used "diff --new-file" to represent the contributed files.

Cheers,
-Pete



diff --new-file --recursive --unified=3 -w signal.orig/src/cl2bp.cc signal/src/cl2bp.cc
--- signal.orig/src/cl2bp.cc 1969-12-31 19:00:00.000000000 -0500
+++ signal/src/cl2bp.cc 2009-10-17 18:19:36.000000000 -0400
@@ -0,0 +1,168 @@
+/*
+
+Copyright (c) 2008-2009, Evgeni A. Nurminski
+Copyright (c) 2008-2009, Pete Gonzalez
+
+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 software; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+
+Authors:
+
+Evgeni A. Nurminski <nurmi@...>
+Institute for Automation and Control Problems
+Far East branch of RAS, Vladivostok, Russia
+
+Pete Gonzalez <pgonzalez@...>
+Bluel Technologies Corporation
+
+*/
+
+#include <octave/oct.h>
+
+#include "cl2bp_lib.h"
+
+static void cancel_handler(void *) {
+  OCTAVE_QUIT;
+}
+
+// When CL2BP_LOGGING is defined, this will direct messages to the Octave pager
+void cl2bp_log(const char *message) {
+  octave_stdout << message;
+}
+
+
+DEFUN_DLD (cl2bp, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{h} =} cl2bp (@var{m}, @var{w1}, @var{w2}, @var{up}, @var{lo}\
+ [, @var{gridsize}])\n\
+\n\
+Constrained L2 bandpass FIR filter design.  This is a fast implementation of the algorithm\n\
+cited below.  Compared to @dfn{remez}, it offers implicit specification of transition bands,\n\
+a higher likelihood of convergence, and an error criterion combining features of both L2 and\n\
+Chebyshev approaches.@*\n\
+Inputs:@*\n\
+@var{m}: degree of cosine polynomial, i.e. the number of output coefficients will be @var{m}*2+1@*\n\
+@var{w1}, @var{w2}: bandpass filter cutoffs in the range 0 <= @var{w1} < @var{w2} <= pi,\n\
+where pi is the Nyquist frequency@*\n\
+@var{up}: vector of 3 upper bounds for [stopband1, passband, stopband2]@*\n\
+@var{lo}: vector of 3 lower bounds for [stopband1, passband, stopband2]@*\n\
+@var{gridsize}: search grid size; larger values may improve accuracy,\n\
+but greatly increase calculation time.@*\n\
+Output:@*\n\
+A vector of @var{m}*2+1 FIR coefficients, or an empty value if the solver failed to converge.@*\n\
+Example:\n\
+@example\n\
+@code{h = cl2bp(30, 0.3*pi, 0.6*pi, [0.02, 1.02, 0.02], [-0.02, 0.98, -0.02], 2^11);}\n\
+@end example\n\
+Original Paper:  I. W. Selesnick, M. Lang, and C. S. Burrus.  A modified algorithm for\n\
+constrained least square design of multiband FIR filters without specified transition bands.\n\
+IEEE Trans. on Signal Processing, 46(2):497-501, February 1998.\n\
+@end deftypefn\n\
+@seealso{remez}\n")
+{
+  octave_value_list retval;
+  int i;
+
+  int nargin = args.length();
+  if (nargin < 5 || nargin > 6) {
+    print_usage ();
+    return retval;
+  }
+
+  int m = args(0).int_value(true);
+  if (error_state) {
+    gripe_wrong_type_arg("cl2bp", args (0));
+    return retval;
+  }
+  double w1 = args(1).double_value();
+  if (error_state) {
+    gripe_wrong_type_arg("cl2bp", args (1));
+    return retval;
+  }
+  double w2 = args(2).double_value();
+  if (error_state) {
+    gripe_wrong_type_arg("cl2bp", args (2));
+    return retval;
+  }
+  ColumnVector up_vector(args(3).vector_value());
+  if (error_state) {
+    gripe_wrong_type_arg("cl2bp", args (3));
+    return retval;
+  }
+  ColumnVector lo_vector(args(4).vector_value());
+  if (error_state) {
+    gripe_wrong_type_arg("cl2bp", args (4));
+    return retval;
+  }
+  if (up_vector.length() != 3 || lo_vector.length() != 3) {
+    error("cl2bp: The up and lo vectors must contain 3 values");
+    return retval;
+  }
+
+  double up[3];
+  double lo[3];
+  for (int i=0; i<3; ++i) {
+    up[i] = up_vector(i);
+    lo[i] = lo_vector(i);
+  }
+
+  int L = args(5).int_value(true);
+  if (error_state) {
+    gripe_wrong_type_arg("cl2bp", args (5));
+    return retval;
+  }
+  if (L > 1000000) {
+    error("cl2bp: The \"gridsize\" parameter cannot exceed 1000000");
+    return retval;
+  }
+  
+  MallocArray<double> h;
+  try {
+    cl2bp(h, m, w1, w2, up, lo, L, 1.e-5, 20, cancel_handler);
+  }
+  catch (std::exception &ex) {
+    error(ex.what());
+    return retval;
+  }
+
+  ColumnVector h_vector(h.get_length());
+
+  for (int i=0; i<h.get_length(); ++i)
+    h_vector(i) = h[i];
+
+  return octave_value(h_vector);
+}
+
+/*
+%!test
+%! b = [
+%!    0.0000000000000000
+%!    0.0563980420304213
+%!   -0.0000000000000000
+%!   -0.0119990278695041
+%!   -0.0000000000000001
+%!   -0.3016146759510104
+%!    0.0000000000000001
+%!    0.5244313235801866
+%!    0.0000000000000001
+%!   -0.3016146759510104
+%!   -0.0000000000000001
+%!   -0.0119990278695041
+%!   -0.0000000000000000
+%!    0.0563980420304213
+%!    0.0000000000000000];
+%! assert(cl2bp(7, 0.25*pi, 0.75*pi, [0.01, 1.04, 0.01], [-0.01, 0.96, -0.01], 2^11), b, 1e-14);
+
+*/
diff --new-file --recursive --unified=3 -w signal.orig/src/cl2bp_lib.cc signal/src/cl2bp_lib.cc
--- signal.orig/src/cl2bp_lib.cc 1969-12-31 19:00:00.000000000 -0500
+++ signal/src/cl2bp_lib.cc 2009-10-17 18:19:40.000000000 -0400
@@ -0,0 +1,616 @@
+/*
+
+Copyright (c) 2008-2009, Evgeni A. Nurminski
+Copyright (c) 2008-2009, Pete Gonzalez
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+Authors:
+
+Evgeni A. Nurminski <nurmi@...>
+Institute for Automation and Control Problems
+Far East branch of RAS, Vladivostok, Russia
+
+Pete Gonzalez <pgonzalez@...>
+Bluel Technologies Corporation
+
+*/
+
+#ifdef _MSC_VER
+#define _CRT_SECURE_NO_WARNINGS  // disable spurious warnings
+#define _USE_MATH_DEFINES // for math.h
+#endif
+
+#include "cl2bp_lib.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <math.h>
+#include <string.h>
+
+#include <stdexcept>
+
+#define BIG_NUMBER 100000
+
+//-----------------------------------------------------------------------------------------------------------
+#ifdef CL2BP_LOGGING
+static void log_printf(const char *format, ...) {
+  va_list argptr;
+  va_start(argptr, format);
+
+  char buf[20*1024];
+  if (vsnprintf(buf,20*1024-1,format,argptr) == -1) {
+    strcpy(buf, "#ERROR#");
+  }
+  va_end(argptr);
+
+  cl2bp_log(buf);
+}
+#else
+#define log_printf(...)  ((void)0)   // don't log anything (for a big speed improvement)
+#endif
+
+//-----------------------------------------------------------------------------------------------------------
+static int local_max(const MallocArray<double>& x, int n, MallocArray<int>& ix) {
+  int i, mx;
+
+  mx = 0;
+
+  MallocArray<double> xx(n+2);
+
+  xx[0] = xx[n + 1] = -BIG_NUMBER;
+  for ( i = 1; i <= n; i++)
+    xx[i] = x[i - 1];
+  for ( i = 0; i < n; i++ ) {
+    (((xx[i] <  xx[i + 1]) && (xx[i + 1] >  xx[i + 2])) ||
+     ((xx[i] <  xx[i + 1]) && (xx[i + 1] >= xx[i + 2])) ||
+     (xx[i] <= xx[i + 1]) && (xx[i + 1] >  xx[i + 2]) )
+    && ( ix[mx++] = i );
+  }
+  return mx;
+}
+
+//-----------------------------------------------------------------------------------------------------------
+static int solvps(MallocArray<double>& a, MallocArray<double>& b, int n,
+  void (*cancel_handler)(void *), void *cancel_state) {
+
+  double t;
+  int j, k;
+  int a_p;
+  int a_q;
+  int a_r;
+  int a_s;
+
+  // In empirical tests, 2% to 6% of the execution time was spent in solvps()
+  if (cancel_handler) cancel_handler(cancel_state);
+
+  for (j = 0, a_p = 0; j < n; ++j, a_p += n+1) {
+    for (a_q = j*n; a_q < a_p; ++a_q)
+      a[a_p] -= a[a_q] * a[a_q];
+    if (a[a_p] <= 0.)
+      return -1;
+    a[a_p] = sqrt(a[a_p]);
+    for (k = j + 1, a_q = a_p + n; k < n; ++k, a_q += n) {
+      for (a_r = j*n, a_s = k*n, t = 0.; a_r < a_p;)
+        t += a[a_r++] * a[a_s++];
+      a[a_q] -= t;
+      a[a_q] /= a[a_p];
+    }
+  }
+  for (j = 0, a_p = 0; j < n; ++j, a_p += n+1) {
+    for (k = 0, a_q = j*n; k < j;)
+      b[j] -=b [k++] * a[a_q++];
+    b[j] /= a[a_p];
+  }
+  for (j = n - 1, a_p = n*n - 1; j >= 0; --j, a_p -= n + 1) {
+    for (k = j + 1, a_q = a_p + n; k < n; a_q += n)
+      b[j] -= b[k++]* a[a_q];
+    b[j] /= a[a_p];
+  }
+  return 0;
+}
+
+//-----------------------------------------------------------------------------------------------------------
+static void rmmult(MallocArray<double>& rm, const MallocArray<double>& a, const MallocArray<double>& b,
+  int n, int m, int l,
+  void (*cancel_handler)(void *), void *cancel_state) {
+
+  double z;
+  int i,j,k;
+  int b_p; // index into b
+  int a_p; // index into a
+  int rm_start = 0;  // index into rm
+  int rm_q; // index into rm
+
+  MallocArray<double> q0(m);
+  for (i = 0; i < l ; ++i, ++rm_start) {
+    // In empirical tests, 87% to 95% of the execution time was spent in rmmult()
+    if (cancel_handler) cancel_handler(cancel_state);
+    for (k = 0, b_p = i; k < m; b_p += l)
+      q0[k++] = b[b_p];
+    for (j = 0, a_p = 0, rm_q = rm_start; j < n; ++j, rm_q += l) {
+      for (k = 0, z = 0.; k < m;)
+        z += a[a_p++] * q0[k++];
+      rm[rm_q]=z;
+    }
+  }
+}
+
+//-----------------------------------------------------------------------------------------------------------
+static void mattr(MallocArray<double>& a, const MallocArray<double>& b, int m, int n) {
+  int i, j;
+  int b_p;
+  int a_p = 0;
+  int b_start = 0;
+  for (i = 0; i < n; ++i, ++b_start)
+    for ( j = 0, b_p = b_start; j < m; ++j, b_p += n )
+      a[a_p++] = b[b_p];
+}
+
+//-----------------------------------------------------------------------------------------------------------
+static void calcAx(MallocArray<double>& Ax, int m, int L) {
+  double r = M_SQRT2, pi = M_PI;
+  int i, j;
+
+  Ax.resize((L+1)*(m + 1));
+
+  for ( i = 0; i <= L; i++)
+    for ( j = 0; j <= m; j++)
+      Ax[i*(m+1) + j] = cos(i*j*pi/L);
+  for ( i = 0; i < (L+1)*(m+1); i += m + 1 )
+    Ax[i] /= r;
+}
+
+//-----------------------------------------------------------------------------------------------------------
+static double L2error(const MallocArray<double>& x, const MallocArray<double>& w, int L, double w1, double w2) {
+  double xx, ww, sumerr, pi = M_PI;
+  int i;
+  for ( i = 0, sumerr = 0; i < L + 1; i++ ) {
+    ww = w[i];
+    xx = x[i];
+    sumerr += ( ww < w1*pi || ww > w2*pi ) ?  xx*xx : (1 - xx)*(1 - xx);
+  }
+  return sumerr;
+}
+
+//-----------------------------------------------------------------------------------------------------------
+static double CHerror(double *wmin, const MallocArray<double>& x, const MallocArray<double>& w,
+  int L, double w1, double w2) {
+
+  double xx, ww, err, errmax;
+  int i;
+  errmax = -1;
+  *wmin = 0;
+  for ( i = 0; i < L + 1; i++ ) {
+    ww = w[i];
+    xx = x[i];
+    err = (( ww < w1 ) || (ww > w2 )) ?  fabs(xx) : fabs(1 - xx);
+    if ( err > errmax ) {
+      errmax = err;
+      *wmin = ww;
+    }
+  }
+  return errmax;
+}
+
+//-----------------------------------------------------------------------------------------------------------
+static void Ggen(MallocArray<double>& G, int m, const MallocArray<double>& w, const MallocArray<int>& kmax,
+  int l_kmax, const MallocArray<int>& kmin, int l_kmin) {
+
+  G.resize((l_kmax + l_kmin)*(m + 1));
+
+  int nsys, i, j;
+  double r = M_SQRT2;
+
+  nsys = l_kmax + l_kmin;
+  for ( i = 0; i < l_kmax; i++)
+    for ( j = 0; j <= m; j++)
+      G[i*(m+1) + j] = cos(w[kmax[i]]*j);
+  for ( i = l_kmax; i < nsys; i++)
+    for ( j = 0; j <= m; j++)
+      G[i*(m+1) + j] = -cos(w[kmin[i - l_kmax]]*j);
+  for ( i = 0; i < nsys*(m+1); i += m+1 )
+    G[i] /= r;
+}
+
+//-----------------------------------------------------------------------------------------------------------
+bool cl2bp(MallocArray<double>& h, int m, double w1, double w2, double up[3], double lo[3], int L,
+  double eps, int mit, void (*cancel_handler)(void *), void *cancel_state) {
+
+  // Ensure sane values for input parameters
+  if (m < 2 || m > 5000)
+    throw std::invalid_argument("cl2bp: The m count must be between 2 and 5000");
+
+  if (w1 < 0 || w1 > w2 || w2 > 2*M_PI)
+    throw std::invalid_argument("cl2bp: The cutoffs must be in the range 0 < w1 < w2 < 2*pi");
+
+  // Z is allocated as Z(2*L-1-2*m)
+  if (L <= m)
+    throw std::invalid_argument("cl2bp: The \"gridsize\" parameter must be larger than \"m\"");
+
+  double r = M_SQRT2;
+  double pi = M_PI;
+  double wmin, ww, xw;
+  int q1, q2, q3, i, iter = 0;
+  int off, neg;
+
+  int ik;
+  int l_kmax;
+  int l_kmin;
+  int l_okmax;
+  int l_okmin;
+  double uvo = 0, lvo = 0;
+
+  double err, diff_up, diff_lo;
+  double erru, erro;
+  int iup, ilo;
+
+  int nsys, j;
+
+  int imin;
+  double umin;
+
+  int k1, k2, ak1, ak2;
+  double cerr;
+
+  h.resize(2*m+1);
+
+  bool converged = true;
+
+  MallocArray<double> x0(L+1);
+  MallocArray<double> x1(L+1);
+  MallocArray<double> xx(L+1);
+  MallocArray<double> xm1(L+1);
+  MallocArray<double> w(L+1);
+  MallocArray<double> u(L+1);
+  MallocArray<double> l(L+1);
+  MallocArray<double> a(L+1);
+  MallocArray<double> c(L+1);
+  MallocArray<int> kmax(L+1);
+  MallocArray<int> kmin(L+1);
+  l_kmax  = l_kmin = 0;
+
+  MallocArray<int> okmin(L+1);
+  MallocArray<int> okmax(L+1);
+  l_okmax = l_okmin = 0;
+  MallocArray<double> rhs(m+2);
+  MallocArray<double> mu(2*(L+1));
+
+  for ( i = 0; i <= L; i++ )
+    w[i] = i*pi/L;
+
+  MallocArray<double> Z(2*L-1-2*m);
+
+  q1 = (int)floor(L*w1/pi);
+  q2 = (int)floor(L*(w2-w1)/pi);
+  q3 = L + 1 - q1 - q2;
+
+  off = 0;
+  for ( i = 0; i < q1; i++) {
+    u[i] = up[0];
+    l[i] = lo[0];
+  }
+  off += i;
+  for ( i = 0; i < q2; i++) {
+    u[off + i] = up[1];
+    l[off + i] = lo[1];
+  }
+  off += i;
+  for ( i = 0; i < q3; i++) {
+    u[off + i] = up[2];
+    l[off + i] = lo[2];
+  }
+
+
+  c[0] = (w2-w1)*r;
+  for ( i = 1; i <= m; i++)
+    c[i] =  2*(sin(w2*i)-sin(w1*i))/i;
+  for ( i = 0; i <= m; i++) {
+    c[i] /=  pi;
+    a[i] = c[i];
+  }
+  log_printf("\nInitial approximation. Unconstrained L2 filter coefficients.\n");
+  log_printf("=============================================================\n");
+
+  log_printf("\nZero order term %8.3lf\n", a[0]);
+  for ( i = 1; i <= m; i++) {
+    log_printf("%4d %8.3lf", i, a[i]);
+    if (i - 8*(i/8) == 0)
+      log_printf("\n");
+  }
+  log_printf("\n");
+
+  // Precalculate Ax
+  MallocArray<double> Ax;
+  calcAx(Ax, m, L);
+
+  //calcA(x0, a, m, L);
+  rmmult(x0, Ax, a, L + 1, m + 1, 1, cancel_handler, cancel_state);
+
+  err = CHerror(&wmin, x0, w, L, w1, w2);
+  log_printf("\nChebyshev err %12.4e (%11.5e)  <> L2 err %12.4e\n", err, wmin/pi, L2error(x0, w, L, w1, w2)/(L+1));
+  for (iter = 1; ; iter++) {
+    log_printf("\n:::::: Iteration %3d :::::::::::\n", iter);
+
+    if ( (uvo > eps*2) || (lvo > eps*2) ) {
+      log_printf("\nXXX Take care of old errors: uvo lvo %12.3e %12.3e", uvo, lvo);
+      if( k1 >= 0 ) log_printf(" k1 %3d(%d)", k1, okmax[k1]);
+      if( k2 >= 0 ) log_printf(" k2 %3d(%d)", k2, okmin[k2]);
+      log_printf("\n");
+
+      if ( uvo > lvo ) {
+        kmax[l_kmax] = okmax[k1];
+        l_kmax++;
+        l_okmax--;
+        for (i = k1; i < l_okmax; i++)
+          okmax[i] = okmax[i + 1];
+      } else {
+        kmin[l_kmin] = okmin[k2];
+        l_kmin++;
+        l_okmin--;
+        for (i = k2; i < l_okmin; i++)
+          okmin[i] = okmin[i + 1];
+      }
+      nsys = l_kmax + l_kmin;
+
+      /*
+        for (i = 0; i < l_kmax; i++)
+          log_printf("i %3d kmax %3d mu %12.4e\n",
+             i, kmax[i], mu[i]);
+        log_printf("\n");
+        for (i = 0; i < l_kmin; i++)
+          log_printf("i %3d kmin %3d mu %12.4e\n",
+             i, kmin[i], mu[i + l_kmax]);
+        log_printf("\n\n");
+      */
+    } else {
+
+      //calcA(x1, a, m, L);
+      rmmult(x1, Ax, a, L + 1, m + 1, 1, cancel_handler, cancel_state);
+
+
+      for ( i = 0; i < l_kmax; i++)
+        okmax[i] = kmax[i];
+      for ( i = 0; i < l_kmin; i++)
+        okmin[i] = kmin[i];
+      l_okmax = l_kmax;
+      l_okmin = l_kmin;
+
+      l_kmax = local_max(x1, L + 1, kmax);
+
+
+      for ( i = 0; i < L + 1; i++)
+        xm1[i] = -x1[i];
+      l_kmin = local_max(xm1, L + 1, kmin);
+
+      log_printf("\nSignificant deviations from the ideal filter. Levels:");
+      log_printf(" %8.2e %8.2e %8.2e (lo)", lo[0], lo[1], lo[2]);
+      log_printf(" %8.2e %8.2e %8.2e (up)", up[0], up[1], up[2]);
+      log_printf("\n");
+
+      for ( i = 0, ik = 0; i < l_kmax; i++) {
+        j = kmax[i];
+        if ( x1[j] > u[j] - 10*eps ) {
+          kmax[ik] = j;
+          ik++;
+        }
+      }
+      l_kmax = ik;
+
+      log_printf("overshoots = %d\n", l_kmax);
+      for ( i = 0; i < l_kmax; i++) {
+        j = kmax[i];
+        ww = w[j];
+        xw = x1[j];
+        err = (w1 < ww && ww < w2) ? xw - 1 : xw;
+        log_printf(" i = %3d kmax = %3d xw = %12.5e err = %10.3e u = %12.4e max = %9.2e\n",
+               i, j, xw, err, u[j], xw - u[j]);
+      }
+
+      for ( i = 0, ik = 0; i < l_kmin; i++) {
+        j = kmin[i];
+        if ( x1[j] < l[j] + 10*eps ) {
+          kmin[ik] = j;
+          ik++;
+        }
+      }
+      l_kmin = ik;
+
+      log_printf("undershoots = %d\n", l_kmin);
+      for ( i = 0; i < l_kmin; i++) {
+        j = kmin[i];
+        ww = w[j];
+        xw = -xm1[j];
+        err =(w1 < ww && ww < w2) ? xw - 1 : xw;
+        log_printf(" i = %3d kmin = %3d xw = %12.5e err = %10.3e l = %12.4e min = %9.2e\n",
+               i, j, xw, err, l[j], xw - l[j]);
+      }
+
+      err = erru = erro = diff_up = diff_lo = 0;
+      iup = ilo = 0;
+      for ( i = 0; i < l_kmax; i++) {
+        ik = kmax[i];
+        diff_up = x1[ik] - u[ik];
+        if ( diff_up > erru ) {
+          erru = diff_up;
+          iup = i;
+        }
+      }
+      for ( i = 0; i < l_kmin; i++) {
+        ik = kmin[i];
+        diff_lo = l[ik] - x1[ik];
+        if ( diff_lo > erro ) {
+          erro = diff_lo;
+          ilo = i;
+        }
+      }
+      err = erro > erru ? erro : erru;
+      log_printf("erru = %14.6e iup = %4d kmax = %4d ", erru, iup, kmax[iup]);
+      log_printf("erro = %14.6e ilo = %4d kmin = %4d\n", erro, ilo, kmin[ilo]);
+
+      if ( err < eps )
+        break;
+    }
+
+
+    nsys = l_kmax + l_kmin;
+
+    MallocArray<double> G(nsys*(m + 1));
+    MallocArray<double> GT(nsys*(m + 1));
+    MallocArray<double> GG(nsys*nsys);
+
+    for ( i = 0; i < l_kmax; i++)
+      for ( j = 0; j <= m; j++)
+        G[i*(m+1) + j] = cos(w[kmax[i]]*j);
+    for ( i = l_kmax; i < nsys; i++)
+      for ( j = 0; j <= m; j++)
+        G[i*(m+1) + j] = -cos(w[kmin[i - l_kmax]]*j);
+    for ( i = 0; i < nsys*(m+1); i += m+1 )
+      G[i] /= r;
+    mattr(GT, G, nsys, m + 1);
+    rmmult(GG, G, GT, nsys, m + 1, nsys, cancel_handler, cancel_state);
+
+
+    rmmult(rhs, G, c, nsys, m + 1, 1, cancel_handler, cancel_state);
+    for ( i = 0; i < nsys; i++)
+      if ( i < l_kmax )
+        rhs[i] -= u[kmax[i]];
+      else
+        rhs[i] += l[kmin[i - l_kmax]];
+
+    for ( i = 0; i < nsys; ++i)
+      mu[i] = rhs[i];
+
+    solvps(GG, mu, nsys, cancel_handler, cancel_state);
+    log_printf("\nXXX KT-system with %d equations resolved.\n", nsys);
+    for ( i = 0, neg = 0; i < nsys; i++)
+      if ( mu[i] < 0 ) neg++;
+    log_printf("\nTotal number of negative multipliers = %3d\n\n", neg);
+    while (neg) {
+
+
+      umin = BIG_NUMBER;
+      for ( i = 0, j = 0; i < nsys; i++) {
+        if ( mu[i] >= 0 ) continue;
+        j++;
+        if ( mu[i] < umin ) {
+          imin = i;
+          umin = mu[i];
+        }
+      }
+
+      if ( umin > 0 )
+        break;
+      log_printf(" neg = %3d of %3d imin = %3d mu-min = %12.4e\n", j, nsys, imin, umin);
+
+      for ( i = imin; i < nsys - 1; i++) {
+        rhs[i] = rhs[i + 1];
+        for ( j = 0; j <= m; j++)
+          G[i*(m + 1) + j] =  G[(i + 1)*(m + 1) + j];
+      }
+
+      if ( imin < l_kmax ) {
+        for ( i = imin; i < l_kmax - 1; i++)
+          kmax[i] = kmax[i + 1];
+        l_kmax--;
+      } else {
+        for ( i = imin; i < nsys - 1; i++)
+          kmin[i - l_kmax] = kmin[i - l_kmax + 1];
+        l_kmin--;
+      }
+
+      --nsys;
+
+      mattr(GT, G, nsys, m + 1);
+      rmmult(GG, G, GT, nsys, m + 1, nsys, cancel_handler, cancel_state);
+      for ( i = 0; i < nsys; ++i)
+        mu[i] = rhs[i];
+      solvps(GG, mu, nsys, cancel_handler, cancel_state);
+    }
+
+    MallocArray<double> da(m + 1);
+    MallocArray<double> zd(nsys);
+
+    rmmult(da, GT, mu, m + 1, nsys, 1, cancel_handler, cancel_state);
+    for ( i = 0; i <= m; i++)
+      a[i] = c[i] - da[i];
+    rmmult(zd, G, a, nsys, m + 1, 1, cancel_handler, cancel_state);
+
+    MallocArray<double> zz(l_okmax + l_okmin);
+    Ggen(G, m, w, okmax, l_okmax, okmin, l_okmin);
+    rmmult(zz, G, a, l_okmax + l_okmin, m + 1, 1, cancel_handler, cancel_state);
+    uvo = lvo = eps;
+    k1 = k2 = -1;
+    ak1 = ak2 = -1;
+    if (l_okmax + l_okmin > 0)
+      log_printf("\nErrors on the previous set of freqs\n\n");
+    for (i = 0; i < l_okmax; i++) {
+      j = okmax[i];
+      cerr = zz[i] - u[j];
+      log_printf(" i %2d j %3d u %12.4e Ga %12.4e cerr %12.4e\n",
+             i, j, u[j], zz[i], cerr);
+      if ( cerr > uvo ) {
+        uvo = cerr;
+        k1 = i;
+        ak1 = j;
+      }
+    }
+    cerr = 0;
+    for (i = 0; i < l_okmin; i++) {
+      j = okmin[i];
+      cerr = l[j] + zz[i + l_okmax];
+      log_printf(" i %2d j %3d l %12.4e Ga %12.4e cerr %12.4e\n",
+             i, j, l[j], zz[i + l_okmax], cerr);
+      if ( cerr > lvo ) {
+        lvo = cerr;
+        k2 = i, ak2 = j;
+      }
+    }
+    if ( l_okmax + l_okmin > 0 ) {
+      log_printf("\n uvo = %12.4e k1 = %4d (%3d) ", uvo, k1, ak1);
+      log_printf("  lvo = %12.4e k2 = %4d (%3d) ", lvo, k2, ak2);
+      log_printf(" maxerr = %12.4e\n", uvo > lvo ? uvo : lvo);
+    }
+
+    log_printf("\nConstrained L2 band filter coefficients.\n");
+    log_printf("=====================================\n");
+
+#ifdef CL2BP_LOGGING
+    log_printf("\nZero order term %8.3lf\n", a[0]);
+    for ( i = 1; i <= m; i++) {
+      log_printf("%4d %8.3lf", i, a[i]);
+      if (i - 8*(i/8) == 0)
+        log_printf("\n");
+    }
+    log_printf("\n");
+#endif
+
+    //calcA(xx, a, m, L);
+    rmmult(xx, Ax, a, L + 1, m + 1, 1, cancel_handler, cancel_state);
+
+    if ( iter >= mit ) {
+      log_printf("Maximum iterations reached\n");
+      converged = false;
+    }
+  }
+  for (i = 0; i < m; i++) {
+    h[i] = a[m - i]/2;
+    h[m + i + 1] = a[i + 1]/2;
+  }
+  h[m] = a[0]*r/2;
+
+  return converged;
+}
diff --new-file --recursive --unified=3 -w signal.orig/src/cl2bp_lib.h signal/src/cl2bp_lib.h
--- signal.orig/src/cl2bp_lib.h 1969-12-31 19:00:00.000000000 -0500
+++ signal/src/cl2bp_lib.h 2009-10-17 18:19:46.000000000 -0400
@@ -0,0 +1,108 @@
+/*
+
+Copyright (c) 2008-2009, Evgeni A. Nurminski
+Copyright (c) 2008-2009, Pete Gonzalez
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Lesser 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 Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+Authors:
+
+Evgeni A. Nurminski <nurmi@...>
+Institute for Automation and Control Problems
+Far East branch of RAS, Vladivostok, Russia
+
+Pete Gonzalez <pgonzalez@...>
+Bluel Technologies Corporation
+
+*/
+
+#ifndef CL2BP_H
+#define CL2BP_H
+
+#include <assert.h>
+
+//-----------------------------------------------------------------------------------------------------------
+// If you want to debug the cl2bp algorithm, define the CL2BP_LOGGING symbol and provide an
+// implementation of cl2bp_log().
+#ifdef CL2BP_LOGGING
+extern void cl2bp_log(const char *message);
+#endif
+
+//-----------------------------------------------------------------------------------------------------------
+// This is a simple helper class that performs bounds-checking for debug builds, but reduces to an unchecked
+// malloc() array for release builds.
+template <class T>
+class MallocArray {
+  int length;
+  T *ptr;
+public:
+  T& operator[](int index) {
+    assert(index >= 0 && index < length);
+    return ptr[index];
+  }
+  T operator[](int index) const {
+    assert(index >= 0 && index < length);
+    return ptr[index];
+  }
+
+  int get_length() const { return length; }
+
+  void resize(int length_) {
+    assert(length_ >= 0 && length_ <= 512*1024*1024);  // verify that the array size is reasonable
+    length = length_;
+    ptr = (T *)realloc(ptr, length * sizeof(T));
+    memset(ptr, 0, length * sizeof(T));
+  }
+
+  MallocArray(int length_=0) {
+    ptr = 0;
+    length = 0;
+    if (length_) resize(length_);
+  }
+  ~MallocArray() { free(ptr); }
+private:
+  MallocArray(const MallocArray& src) { }  // copy constructor is unimplemented and disallowed
+};
+
+//-----------------------------------------------------------------------------------------------------------
+// Constrained L2 BandPass FIR filter design
+//  h       2*m+1 filter coefficients (output)
+//  m       degree of cosine polynomial
+//  w1,w2   fist and second band edges
+//  up      upper bounds
+//  lo      lower bounds
+//  L       grid size
+//  eps     stopping condition ("SN")
+//  mit     maximum number of iterations
+//
+// cl2bp() returns true if the solver converged, or false if the maximum number of iterations was reached.
+// If provided, the cancel_handler function pointer will be called periodically inside long-running loops,
+// giving an opportunity to abort (by throwing a C++ exception).  The cancel_state parameter is a
+// user-defined pointer, e.g. a button object, boolean flag, or other means of detecting a cancel request.
+//
+// Example usage:
+//   MallocArray<double> coefficients;
+//   double up[3] = { 0.02, 1.02, 0.02 };
+//   double lo[3] = { -0.02, 0.98, -0.02 };
+//   cl2bp(coefficients, 30, 0.3*M_PI, 0.6*M_PI, up, lo, 1<<11, 1.e-5, 20);
+//
+// The algorithm is described in this paper:
+//   I. W. Selesnick, M. Lang, and C. S. Burrus.  A modified algorithm for constrained least square
+//   design of multiband FIR filters without specified transition bands.  IEEE Trans. on Signal
+//   Processing, 46(2):497-501, February 1998.
+bool cl2bp(MallocArray<double>& h, int m, double w1, double w2, double up[3], double lo[3], int L,
+  double eps, int mit, void (*cancel_handler)(void *)=0, void *cancel_state=0);
+
+#endif
diff --new-file --recursive --unified=3 -w signal.orig/src/Makefile signal/src/Makefile
--- signal.orig/src/Makefile 2009-10-17 18:14:46.000000000 -0400
+++ signal/src/Makefile 2009-10-17 18:32:15.000000000 -0400
@@ -1,6 +1,18 @@
-all: remez.oct medfilt1.oct sosfilt.oct upfirdn.oct
+all: cl2bp.oct remez.oct medfilt1.oct sosfilt.oct upfirdn.oct
+
+CL2BP_DEFINES =
+# CL2BP_DEFINES = -DCL2BP_LOGGING
+
+cl2bp.o: cl2bp.cc cl2bp_lib.h
+ mkoctfile $(CL2BP_DEFINES) -c $<
+
+cl2bp_lib.o: cl2bp_lib.cc cl2bp_lib.h
+ mkoctfile $(CL2BP_DEFINES) -c $<
+
+cl2bp.oct: cl2bp.o cl2bp_lib.o
+ mkoctfile -s $^
 
 %.oct: %.cc
  mkoctfile -s $<
 
-clean: ; -rm *.o core octave-core *.oct *~
+clean: ; -rm -f *.o core octave-core *.oct *~

------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Pete Gonzalez]

Oops, I forgot to update the INDEX file.  The cl2bp function belongs in
the "FIR filter design" section.

However, I noticed that the INDEX file includes various functions which
are not actually part of the Octave distribution: fircls, fircls1,
firrcos, intfilt, etc.  Is this normal?

Cheers,
-Pete


Pete Gonzalez wrote:

> Søren Hauberg wrote:
>  > It seems nobody objects to your suggested code, so I guess it is fine
>  > for inclusion. Could you send a patch that integrates the code with the
>  > package, i.e. it should also make needed changes to 'Makefile' and
>  > 'INDEX' ?
>
> Attached is a patch against the svn HEAD revision for this directory:
>
> https://octave.svn.sourceforge.net/svnroot/octave/trunk/octave-forge/main/signal 
>
>
> I used "diff --new-file" to represent the contributed files.
>
> Cheers,
> -Pete
>
>

------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Søren Hauberg]

lør, 17 10 2009 kl. 18:54 -0400, skrev Pete Gonzalez:
> Oops, I forgot to update the INDEX file.  The cl2bp function belongs in
> the "FIR filter design" section.

I commited your patch; thanks! I made some small changes to fix some
compiler warnings. Could you check that they are okay? In
'cl2bp_lib.cc:cl2bp' I initialised

  imin = BIG_NUMBER
  k1 = k2 = -1

is that okay? I also added extra parenthesis in line 80 (in
'local_max'). Are they okay?

I added 'cl2bp' to the INDEX file.

> However, I noticed that the INDEX file includes various functions which
> are not actually part of the Octave distribution: fircls, fircls1,
> firrcos, intfilt, etc.  Is this normal?

These functions then show up as "not implemented" in the Octave-Forge
website. It highlights what is missing.

Søren


------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Pete Gonzalez]

Søren Hauberg wrote:
> I commited your patch; thanks! I made some small changes to fix some
> compiler warnings. Could you check that they are okay? In
> 'cl2bp_lib.cc:cl2bp' I initialised
>
>   imin = BIG_NUMBER
>   k1 = k2 = -1
 > is that okay?

Looks good.

> I also added extra parenthesis in line 80 (in
> 'local_max'). Are they okay?

The "&&" operator has a higher precedence than "||", so your parentheses
don't alter the semantics.  Since many people haven't memorized this
table of like 15 different precedence levels, it definitely improves
readability.

Of course, there are a lot of other obvious cleanups that could be made.
  I hope to do that eventually, but actually the code has already come
pretty far.  For programmers, complexity is the enemy, but for
mathematicians I guess it's a central feature of the job.  ;-)

I also noticed that you commented out the L2error() function.  This
function is required by the CL2BP_LOGGING switch, so I would suggest an
#ifdef instead of a comment.  (But I question the premise of
-Wunused-function, that these functions somehow imply a problem.)

One other note:  It would have been easier for me to analyze your
changes if you had applied them as a separate SVN revision.  Just a
suggestion.

Cheers,
-Pete

------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Søren Hauberg]

søn, 18 10 2009 kl. 20:30 -0400, skrev Pete Gonzalez:
> I also noticed that you commented out the L2error() function.  This
> function is required by the CL2BP_LOGGING switch, so I would suggest an
> #ifdef instead of a comment.  (But I question the premise of
> -Wunused-function, that these functions somehow imply a problem.)

Done

> One other note:  It would have been easier for me to analyze your
> changes if you had applied them as a separate SVN revision.  Just a
> suggestion.

Yeah, I wanted to do that, but I forgot to add the new files to SVN
before editing them :-(

Søren


------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev