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 <[email protected]>
Institute for Automation and Control Problems
Far East branch of RAS, Vladivostok, Russia
Pete Gonzalez <[email protected]>
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 approach...@*\n\
Inputs:@*\n\
@var{m}: degree of cosine polynomial, i.e. the number of output coefficients
will be @var{m}*...@*\n\
@var{w1}, @var{w2}: bandpass filter cutoffs in the range 0 <= @var{w1} <
@var{w2} <= pi,\n\
where pi is the Nyquist freque...@*\n\
@var{up}: vector of 3 upper bounds for [stopband1, passband, stopban...@*\n\
@var{lo}: vector of 3 lower bounds for [stopband1, passband, stopban...@*\n\
@var{gridsize}: search grid size; larger values may improve accuracy,\n\
but greatly increase calculation ti...@*\n\
Output:@*\n\
A vector of @var{m}*2+1 FIR coefficients, or an empty value if the solver
failed to conver...@*\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 <[email protected]>
Institute for Automation and Control Problems
Far East branch of RAS, Vladivostok, Russia
Pete Gonzalez <[email protected]>
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 <[email protected]>
Institute for Automation and Control Problems
Far East branch of RAS, Vladivostok, Russia
Pete Gonzalez <[email protected]>
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
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
