[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Possible contribution for the "signal" package
|
From: |
Pete Gonzalez |
|
Subject: |
Possible contribution for the "signal" package |
|
Date: |
Mon, 28 Sep 2009 08:23:25 -0400 |
|
User-agent: |
Thunderbird 2.0.0.23 (Windows/20090812) |
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.
Here's a brief explanation of why this code is interesting:
Books often create the impression that Finite Impulse Response (FIR)
filters are inferior to Infinite Impulse Response (IIR), since IIR
filters are more advanced. (More advanced is better, right?) But
actually FIR filters have a number of advantages:
1.6.1 What are the advantages of FIR Filters compared to IIR filters?
http://dspguru.com/info/faqs/fir/basics.htm
There are two traditional approaches for designing FIR filters. The
simple approach is called the windowing or "frequency sampling" method.
It obtains the coefficients via an inverse DFT, possibly applying a
window function to improve behavior between the DFT points. This is
fir1.m/fir2.m in Octave. It's very stable, but produces suboptimal
results because there is no formal error criterion.
The other traditional approach is the Parks-McClellan solver published
in 1972, based on the Remez exchange algorithm. This is remez.cc in
Octave. It produces significantly better results (optimal in the
Chebyshev sense), but the algorithm frequently fails to converge. It
also requires careful manual tuning of transition bands, which was
unacceptable for our application, since we needed to generate FIR
filters at runtime (i.e. without human oversight).
Despite these shortcomings, these two approaches are the canon of
classic texts like Lyons' "Understanding DSP", Rorabaugh's "DSP Primer",
etc. But in the 1990's a number of innovative new algorithms were
developed for FIR design. Sidney Burrus summarizes the results in his
online book draft:
"Constrained Approximation and Mixed Criteria"
http://cnx.org/content/m16923/latest/
If you don't want to babysit your algorithm or carefully tweak its
results, the most interesting innovation is the Constrained Least
Squares (CLS) method. It has three key advantages: First, although
it's not proven, in practice CLS always converges for reasonable error
constraints. Second, you don't need to specify transition bands -- it
implicitly finds the smallest transition bands that satisfy the
specified error, which makes it very "user friendly". Lastly, the error
metric provides a theoretical balance between two extremes (see Figure 1
in the above link), which is helpful for general applications.
Unfortunately Octave's "signal" package does not currently include a CLS
choice for FIR design. We have created an optimized C++ implementation
which we're offering to contribute as an Octave extension under the LGPL
license. As far as I know, this is the first implementation of its
kind. It's the bandpass version, which will also solve highpass and
lowpass designs (simply by setting a cutoff to 0 or pi). It would not
be difficult to generalize the code for the multiband piecewise
scenario, but I suspect that relatively few applications would actually
benefit from that.
Anyway, I think this is a valuable addition to an engineer's toolbox,
and I hope that you find it useful. I didn't see an official guideline
for Octave submissions, so if you need any minor changes to comply with
your conventions, let me know.
Please CC any replies to the mailing list, since unknown e-mail
addresses often end up in my spam folder. ;-)
Thanks,
Pete Gonzalez
/*
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 <address@hidden>
Institute for Automation and Control Problems
Far East branch of RAS, Vladivostok, Russia
Pete Gonzalez <address@hidden>
Bluel Technologies Corporation
*/
#include <octave/oct.h>
//===========================================================================================================
#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>
#define BIG_NUMBER 100000
#if 0 // enable/disable debug logging
//-----------------------------------------------------------------------------------------------------------
static void debug_print(const char *format, ...) {
va_list argptr;
va_start(argptr, format);
static 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
#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 count;
T *ptr;
public:
T& operator[](int index) {
assert(index >= 0 && index < count);
return ptr[index];
}
const T& operator[](int index) const {
assert(index >= 0 && index < count);
return ptr[index];
}
int get_count() const { return count; }
void resize(int count_) {
assert(count_ >= 0 && count_ <= 10*1024*1024); // verify that the array
size is reasonable
count = count_;
ptr = (T *)realloc(ptr, count * sizeof(T));
memset(ptr, 0, count * sizeof(T));
}
MallocArray(int count_) {
ptr = 0;
resize(count_);
}
~MallocArray() { free(ptr); }
private:
MallocArray(const MallocArray& src) { } // copy constructor is unimplemented
and disallowed
};
static void cl2bp(double *h, int m, double w1, double w2, double up[3], double
lo[3], int L, double eps, int mit);
//-----------------------------------------------------------------------------------------------------------
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) {
double t;
int j, k;
int a_p;
int a_q;
int a_r;
int a_s;
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) {
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) {
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) {
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(MallocArray<double>& 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;
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];
}
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
MallocArray<double> Ax(0);
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 <= mit; 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;
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);
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);
}
MallocArray<double> da(m + 1);
MallocArray<double> zd(nsys);
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);
G.resize((l_okmax + l_okmin)*(m + 1));
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);
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");
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");
//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;
}
//===========================================================================================================
/* == Octave interface starts here ====================================== */
DEFUN_DLD (cl2bp, args, ,
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} address@hidden =} cl2bp (@var{m}, @var{w1},
@var{w2}, @var{up}, @var{lo} [, @var{gridsize}])\n\
Constrained L2 BandPass FIR filter address@hidden
Inputs:@*\n\
@var{m}: degree of cosine polynomial, i.e. the number of output coefficients
will be @address@hidden
@var{w1}, @var{w2}: bandpass filter cutoffs in the range 0 <= @var{w1} <
@var{w2} <= pi, where pi is the Nyquist address@hidden
@var{up}: vector of 3 upper bounds for [stopband1, passband, address@hidden
@var{lo}: vector of 3 lower bounds for [stopband1, passband, address@hidden
@var{gridsize}: search grid size; larger values may improve accuracy, but
greatly increase calculation time.\n\
The default/recommended grid size is address@hidden
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")
{
octave_value_list retval;
int i;
int nargin = args.length();
if (nargin < 5 || nargin > 6) {
print_usage ();
return retval;
}
int m = NINT(args(0).double_value());
if (m < 1 || m > 2000) {
error("cl2bp: The m count must be between 1 and 2000");
return retval;
}
double w1 = args(1).double_value();
double w2 = args(2).double_value();
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());
ColumnVector lo_vector(args(4).vector_value());
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 = NINT(args(5).double_value());
if (L > 100000) {
error("cl2bp: The gridsize L must be less than 100000");
return retval;
}
// MallocArray<double> Z(2*L-1-2*m);
if (2*L-1-2*m <= 0) {
error("cl2bp: The gridsize is too small for the given number of
coefficients");
return retval;
}
MallocArray<double> h(m*2+1);
cl2bp(h, m, w1, w2, up, lo, L, 1.e-5, 20);
ColumnVector h_vector(h.get_count());
for (int i=0; i<h.get_count(); ++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);
*/
- Possible contribution for the "signal" package,
Pete Gonzalez <=