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 <[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
//===========================================================================================================
#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 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;
}
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
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
