Update of /cvsroot/pure-data/externals/iem/iem_adaptfilt/src In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv7807
Modified Files: iem_adaptfilt.c makefile_lin makefile.txt makefile_win makefile_win.txt Added Files: n_CLNLMS~.c n_CNLMS~.c NLMS~.c NLMSCC~.c Log Message: changed files to new naming convention
Index: makefile_lin =================================================================== RCS file: /cvsroot/pure-data/externals/iem/iem_adaptfilt/src/makefile_lin,v retrieving revision 1.1 retrieving revision 1.2 diff -C2 -d -r1.1 -r1.2 *** makefile_lin 2 Aug 2006 14:02:28 -0000 1.1 --- makefile_lin 11 Jan 2007 17:40:10 -0000 1.2 *************** *** 3,7 **** .SUFFIXES: .pd_linux
! INCLUDE = -I. -I/usr/local/src/pd-0.37-1/src
LDFLAGS = -export-dynamic -shared --- 3,7 ---- .SUFFIXES: .pd_linux
! INCLUDE = -I. -I/usr/local/src/pd/src
LDFLAGS = -export-dynamic -shared *************** *** 10,14 **** #select either the DBG and OPT compiler flags below:
! CFLAGS = -DPD -DUNIX -W -Werror -Wno-unused \ -Wno-parentheses -Wno-switch -O6 -funroll-loops -fomit-frame-pointer \ -DDL_OPEN --- 10,14 ---- #select either the DBG and OPT compiler flags below:
! CFLAGS = -DPD -DUNIX -W -Wno-unused \ -Wno-parentheses -Wno-switch -O6 -funroll-loops -fomit-frame-pointer \ -DDL_OPEN *************** *** 18,25 **** # the sources
! SRC = sigNLMS.c \ ! sigNLMSCC.c \ ! sign_CNLMS.c \ ! sign_CLNLMS.c \ iem_adaptfilt.c
--- 18,25 ---- # the sources
! SRC = NLMS~.c \ ! NLMSCC~.c \ ! n_CNLMS~.c \ ! n_CLNLMS~.c \ iem_adaptfilt.c
Index: makefile_win.txt =================================================================== RCS file: /cvsroot/pure-data/externals/iem/iem_adaptfilt/src/makefile_win.txt,v retrieving revision 1.1 retrieving revision 1.2 diff -C2 -d -r1.1 -r1.2 *** makefile_win.txt 2 Aug 2006 14:02:28 -0000 1.1 --- makefile_win.txt 11 Jan 2007 17:40:10 -0000 1.2 *************** *** 22,29 ****
! SRC = sigNLMS.c \ ! sigNLMSCC.c \ ! sign_CNLMS.c \ ! sign_CLNLMS.c \ iem_adaptfilt.c
--- 22,29 ----
! SRC = NLMS~.c \ ! NLMSCC~.c \ ! n_CNLMS~.c \ ! n_CLNLMS~.c \ iem_adaptfilt.c
--- NEW FILE: n_CNLMS~.c --- /* For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution.
n_CNLMS multichannel-constrained (non leaky) normalized LMS algorithm lib iem_adaptfilt written by Markus Noisternig & Thomas Musil noisternig_AT_iem.at; musil_AT_iem.at (c) Institute of Electronic Music and Acoustics, Graz Austria 2005 */
#ifdef NT #pragma warning( disable : 4244 ) #pragma warning( disable : 4305 ) #endif
#include "m_pd.h" #include "iemlib.h" #include <math.h> #include <stdio.h> #include <string.h>
/* ----------------------- n_CNLMS~ ------------------------------ */ /* -- multi-channel Constraint Normalized Least Mean Square (linear adaptive FIR-filter) -- */
/* -- first input: reference signal -- */ /* -- second input: desired signal -- */ /* -- -- */
/* for further information on adaptive filter design we refer to */ /* [1] Haykin, "Adaptive Filter Theory", 4th ed, Prentice Hall */ /* [2] Benesty, "Adaptive Signal Processing", Springer */
typedef struct n_CNLMS_tilde_kern { t_symbol *x_w_array_sym_name; t_float *x_w_array_mem_beg; t_float *x_in_ptr_beg;// memory: sig-in vector t_float *x_out_ptr_beg;// memory: sig-out vector t_float *x_in_hist;// start point double buffer for sig-in history } t_n_CNLMS_tilde_kern;
typedef struct n_CNLMS_tilde { t_object x_obj; t_n_CNLMS_tilde_kern *x_my_kern; t_float *x_des_in_ptr_beg;// memory: desired-in vector t_float *x_err_out_ptr_beg;// memory: error-out vector t_int x_n_io;// number of in-channels and filtered out-channels t_int x_rw_index;// read-write-index t_int x_n_order;// filter order t_int x_update;// rounded by 2^n, yields downsampling of update rate t_float x_beta;// learn rate [0 .. 2] t_float x_gamma;// normalization t_float x_kappa;// constraint: threshold of energy (clipping) t_outlet *x_out_compressing_bang; t_clock *x_clock; t_float x_msi; } t_n_CNLMS_tilde;
t_class *n_CNLMS_tilde_class;
static void n_CNLMS_tilde_tick(t_n_CNLMS_tilde *x) { outlet_bang(x->x_out_compressing_bang); }
static t_float *n_CNLMS_tilde_check_array(t_symbol *array_sym_name, t_int length) { t_int n_points; t_garray *a; t_float *vec;
if(!(a = (t_garray *)pd_findbyclass(array_sym_name, garray_class))) { error("%s: no such array for n_CNLMS~", array_sym_name->s_name); return((t_float *)0); } else if(!garray_getfloatarray(a, &n_points, &vec)) { error("%s: bad template for n_CNLMS~", array_sym_name->s_name); return((t_float *)0); } else if(n_points < length) { error("%s: bad array-size for n_CNLMS~: %d", array_sym_name->s_name, n_points); return((t_float *)0); } else { return(vec); } }
static void n_CNLMS_tilde_beta(t_n_CNLMS_tilde *x, t_floatarg f) // learn rate { if(f < 0.0f) f = 0.0f; if(f > 2.0f) f = 2.0f;
x->x_beta = f; }
static void n_CNLMS_tilde_gamma(t_n_CNLMS_tilde *x, t_floatarg f) // regularization (dither) { if(f < 0.0f) f = 0.0f; if(f > 1.0f) f = 1.0f;
x->x_gamma = f; }
static void n_CNLMS_tilde_kappa(t_n_CNLMS_tilde *x, t_floatarg f) // threshold for w_coeff { if(f < 0.0001f) f = 0.0001f; if(f > 10000.0f) f = 10000.0f;
x->x_kappa = f; }
static void n_CNLMS_tilde_update(t_n_CNLMS_tilde *x, t_floatarg f) // downsampling of learn rate { t_int i=1, u = (t_int)f;
if(u < 0) u = 0; else { while(i <= u) // convert u for 2^N i *= 2; // round downward i /= 2; u = i; } x->x_update = u; }
/* ============== DSP ======================= */
static t_int *n_CNLMS_tilde_perform_zero(t_int *w) { t_n_CNLMS_tilde *x = (t_n_CNLMS_tilde *)(w[1]); t_int n = (t_int)(w[2]);
t_int n_io = x->x_n_io; t_float *out; t_int i, j;
out = x->x_err_out_ptr_beg; for(i=0; i<n; i++) *out++ = 0.0f; for(j=0; j<n_io; j++) { out = x->x_my_kern[j].x_out_ptr_beg; for(i=0; i<n; i++) *out++ = 0.0f; } return (w+3); }
static t_int *n_CNLMS_tilde_perform(t_int *w) { t_n_CNLMS_tilde *x = (t_n_CNLMS_tilde *)(w[1]); t_int n = (t_int)(w[2]); t_int n_order = x->x_n_order; /* filter-order */ t_int rw_index2, rw_index = x->x_rw_index; t_int n_io = x->x_n_io; t_float *in;// first sig in t_float din;// second sig in t_float *filt_out;// first sig out t_float *err_out, err_sum;// second sig out t_float *read_in_hist; t_float *w_filt_coeff; t_float my, my_err, sum; t_float beta = x->x_beta; t_float hgamma, gamma = x->x_gamma; t_float hkappa, kappa = x->x_kappa; t_int i, j, k, update_counter; t_int update = x->x_update; t_int ord8=n_order&0xfffffff8; t_int ord_residual=n_order&0x7; t_int compressed = 0;
for(k=0; k<n_io; k++) { if(!x->x_my_kern[k].x_w_array_mem_beg) goto n_CNLMS_tildeperfzero;// this is Musil/Miller style }
hgamma = gamma * gamma * (float)n_order; //hkappa = kappa * kappa * (float)n_order; hkappa = kappa;// kappa regards to energy value, else use line above
for(i=0, update_counter=0; i<n; i++)// history and (block-)convolution { rw_index2 = rw_index + n_order;
for(k=0; k<n_io; k++)// times n_io { x->x_my_kern[k].x_in_hist[rw_index] = x->x_my_kern[k].x_in_ptr_beg[i]; // save inputs into variabel & history x->x_my_kern[k].x_in_hist[rw_index+n_order] = x->x_my_kern[k].x_in_ptr_beg[i]; } din = x->x_des_in_ptr_beg[i];
// begin convolution err_sum = din; for(k=0; k<n_io; k++)// times n_io { sum = 0.0f; w_filt_coeff = x->x_my_kern[k].x_w_array_mem_beg; // Musil's special convolution buffer struct read_in_hist = &x->x_my_kern[k].x_in_hist[rw_index2]; for(j=0; j<ord8; j+=8) // loop unroll 8 taps { sum += w_filt_coeff[0] * read_in_hist[0]; sum += w_filt_coeff[1] * read_in_hist[-1]; sum += w_filt_coeff[2] * read_in_hist[-2]; sum += w_filt_coeff[3] * read_in_hist[-3]; sum += w_filt_coeff[4] * read_in_hist[-4]; sum += w_filt_coeff[5] * read_in_hist[-5]; sum += w_filt_coeff[6] * read_in_hist[-6]; sum += w_filt_coeff[7] * read_in_hist[-7]; w_filt_coeff += 8; read_in_hist -= 8; } for(j=0; j<ord_residual; j++) // for filter order < 2^N sum += w_filt_coeff[j] * read_in_hist[-j];
x->x_my_kern[k].x_out_ptr_beg[i] = sum; err_sum -= sum; } x->x_err_out_ptr_beg[i] = err_sum; // end convolution
if(update) // downsampling of learn rate { update_counter++; if(update_counter >= update) { update_counter = 0;
for(k=0; k<n_io; k++)// times n_io { sum = 0.0f;// calculate energy for last n-order samples in filter read_in_hist = &x->x_my_kern[k].x_in_hist[rw_index2]; for(j=0; j<ord8; j+=8) // unrolling quadrature calc { sum += read_in_hist[0] * read_in_hist[0]; sum += read_in_hist[-1] * read_in_hist[-1]; sum += read_in_hist[-2] * read_in_hist[-2]; sum += read_in_hist[-3] * read_in_hist[-3]; sum += read_in_hist[-4] * read_in_hist[-4]; sum += read_in_hist[-5] * read_in_hist[-5]; sum += read_in_hist[-6] * read_in_hist[-6]; sum += read_in_hist[-7] * read_in_hist[-7]; read_in_hist -= 8; } for(j=0; j<ord_residual; j++) // residual sum += read_in_hist[-j] * read_in_hist[-j]; // [-j] only valid for Musil's double buffer structure sum += hgamma; // convert gamma corresponding to filter order my = beta / sum;// calculate mue
my_err = my * err_sum; w_filt_coeff = x->x_my_kern[k].x_w_array_mem_beg; read_in_hist = &x->x_my_kern[k].x_in_hist[rw_index2]; sum = 0.0f; for(j=0; j<ord8; j+=8) // unrolling quadrature calc { w_filt_coeff[0] += read_in_hist[0] * my_err; sum += w_filt_coeff[0] * w_filt_coeff[0]; w_filt_coeff[1] += read_in_hist[-1] * my_err; sum += w_filt_coeff[1] * w_filt_coeff[1]; w_filt_coeff[2] += read_in_hist[-2] * my_err; sum += w_filt_coeff[2] * w_filt_coeff[2]; w_filt_coeff[3] += read_in_hist[-3] * my_err; sum += w_filt_coeff[3] * w_filt_coeff[3]; w_filt_coeff[4] += read_in_hist[-4] * my_err; sum += w_filt_coeff[4] * w_filt_coeff[4]; w_filt_coeff[5] += read_in_hist[-5] * my_err; sum += w_filt_coeff[5] * w_filt_coeff[5]; w_filt_coeff[6] += read_in_hist[-6] * my_err; sum += w_filt_coeff[6] * w_filt_coeff[6]; w_filt_coeff[7] += read_in_hist[-7] * my_err; sum += w_filt_coeff[7] * w_filt_coeff[7]; w_filt_coeff += 8; read_in_hist -= 8; } for(j=0; j<ord_residual; j++) // residual { w_filt_coeff[j] += read_in_hist[-j] * my_err; sum += w_filt_coeff[j] * w_filt_coeff[j]; } if(sum > hkappa) { compressed = 1; my = sqrt(hkappa/sum); w_filt_coeff = x->x_my_kern[k].x_w_array_mem_beg; for(j=0; j<ord8; j+=8) // unrolling quadrature calc { w_filt_coeff[0] *= my; w_filt_coeff[1] *= my; w_filt_coeff[2] *= my; w_filt_coeff[3] *= my; w_filt_coeff[4] *= my; w_filt_coeff[5] *= my; w_filt_coeff[6] *= my; w_filt_coeff[7] *= my; w_filt_coeff += 8; } for(j=0; j<ord_residual; j++) // residual w_filt_coeff[j] *= my; } } } } rw_index++; if(rw_index >= n_order) rw_index -= n_order; }
x->x_rw_index = rw_index; // back to start
if(compressed) clock_delay(x->x_clock, 0);
return(w+3);
n_CNLMS_tildeperfzero:
err_out = x->x_err_out_ptr_beg; for(i=0; i<n; i++) *err_out++ = 0.0f; for(j=0; j<n_io; j++) { filt_out = x->x_my_kern[j].x_out_ptr_beg; for(i=0; i<n; i++) *filt_out++ = 0.0f; }
return(w+3); }
static void n_CNLMS_tilde_dsp(t_n_CNLMS_tilde *x, t_signal **sp) { t_int i, n = sp[0]->s_n; t_int ok_w = 1; t_int m = x->x_n_io;
for(i=0; i<m; i++) x->x_my_kern[i].x_in_ptr_beg = sp[i]->s_vec; x->x_des_in_ptr_beg = sp[m]->s_vec; for(i=0; i<m; i++) x->x_my_kern[i].x_out_ptr_beg = sp[i+m+1]->s_vec; x->x_err_out_ptr_beg = sp[2*m+1]->s_vec;
for(i=0; i<m; i++) { x->x_my_kern[i].x_w_array_mem_beg = n_CNLMS_tilde_check_array(x->x_my_kern[i].x_w_array_sym_name, x->x_n_order); if(!x->x_my_kern[i].x_w_array_mem_beg) ok_w = 0; }
if(!ok_w) dsp_add(n_CNLMS_tilde_perform_zero, 2, x, n); else dsp_add(n_CNLMS_tilde_perform, 2, x, n); }
/* setup/setdown things */
static void n_CNLMS_tilde_free(t_n_CNLMS_tilde *x) { t_int i, n_io=x->x_n_io, n_order=x->x_n_order;
for(i=0; i<n_io; i++) freebytes(x->x_my_kern[i].x_in_hist, 2*x->x_n_order*sizeof(t_float)); freebytes(x->x_my_kern, n_io*sizeof(t_n_CNLMS_tilde_kern));
clock_free(x->x_clock); }
static void *n_CNLMS_tilde_new(t_symbol *s, t_int argc, t_atom *argv) { t_n_CNLMS_tilde *x = (t_n_CNLMS_tilde *)pd_new(n_CNLMS_tilde_class); char buffer[400]; t_int i, n_order=39, n_io=1; t_symbol *w_name; t_float beta=0.1f; t_float gamma=0.00001f; t_float kappa = 1.0f;
if((argc >= 6) && IS_A_FLOAT(argv,0) && //IS_A_FLOAT/SYMBOL from iemlib.h IS_A_FLOAT(argv,1) && IS_A_FLOAT(argv,2) && IS_A_FLOAT(argv,3) && IS_A_FLOAT(argv,4) && IS_A_SYMBOL(argv,5)) { n_io = (t_int)atom_getintarg(0, argc, argv); n_order = (t_int)atom_getintarg(1, argc, argv); beta = (t_float)atom_getfloatarg(2, argc, argv); gamma = (t_float)atom_getfloatarg(3, argc, argv); kappa = (t_float)atom_getfloatarg(4, argc, argv); w_name = (t_symbol *)atom_getsymbolarg(5, argc, argv);
if(beta < 0.0f) beta = 0.0f; if(beta > 2.0f) beta = 2.0f;
if(gamma < 0.0f) gamma = 0.0f; if(gamma > 1.0f) gamma = 1.0f;
if(kappa < 0.0001f) kappa = 0.0001f; if(kappa > 10000.0f) kappa = 10000.0f;
if(n_order < 2) n_order = 2; if(n_order > 11111) n_order = 11111;
if(n_io < 1) n_io = 1; if(n_io > 60) n_io = 60;
for(i=0; i<n_io; i++) inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_signal, &s_signal); for(i=0; i<=n_io; i++) outlet_new(&x->x_obj, &s_signal);
x->x_out_compressing_bang = outlet_new(&x->x_obj, &s_bang);
x->x_msi = 0; x->x_n_io = n_io; x->x_n_order = n_order; x->x_update = 0; x->x_beta = beta; x->x_gamma = gamma; x->x_kappa = kappa; x->x_my_kern = (t_n_CNLMS_tilde_kern *)getbytes(x->x_n_io*sizeof(t_n_CNLMS_tilde_kern)); for(i=0; i<n_io; i++) { sprintf(buffer, "%d_%s", i+1, w_name->s_name); x->x_my_kern[i].x_w_array_sym_name = gensym(buffer); x->x_my_kern[i].x_w_array_mem_beg = (t_float *)0; x->x_my_kern[i].x_in_hist = (t_float *)getbytes(2*x->x_n_order*sizeof(t_float)); } x->x_clock = clock_new(x, (t_method)n_CNLMS_tilde_tick);
return(x); } else { post("n_CNLMSC~-ERROR: need 5 float- + 1 symbol-arguments:"); post(" number_of_filters + order_of_filters + learnrate_beta + security_value_gamma + threshold_kappa + array_name_taps"); return(0); } }
void n_CNLMS_tilde_setup(void) { n_CNLMS_tilde_class = class_new(gensym("n_CNLMS~"), (t_newmethod)n_CNLMS_tilde_new, (t_method)n_CNLMS_tilde_free, sizeof(t_n_CNLMS_tilde), 0, A_GIMME, 0); CLASS_MAINSIGNALIN(n_CNLMS_tilde_class, t_n_CNLMS_tilde, x_msi); class_addmethod(n_CNLMS_tilde_class, (t_method)n_CNLMS_tilde_dsp, gensym("dsp"), 0); class_addmethod(n_CNLMS_tilde_class, (t_method)n_CNLMS_tilde_update, gensym("update"), A_FLOAT, 0); // method: downsampling factor of learning (multiple of 2^N) class_addmethod(n_CNLMS_tilde_class, (t_method)n_CNLMS_tilde_beta, gensym("beta"), A_FLOAT, 0); //method: normalized learning rate class_addmethod(n_CNLMS_tilde_class, (t_method)n_CNLMS_tilde_gamma, gensym("gamma"), A_FLOAT, 0); // method: dithering noise related to signal class_addmethod(n_CNLMS_tilde_class, (t_method)n_CNLMS_tilde_kappa, gensym("kappa"), A_FLOAT, 0); // method: threshold for compressing w_coeff
//class_sethelpsymbol(n_CNLMS_tilde_class, gensym("iemhelp2/n_CNLMS~")); }
--- NEW FILE: n_CLNLMS~.c --- /* For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution.
n_CLNLMS multichannel-constrained leaky normalized LMS algorithm lib iem_adaptfilt written by Markus Noisternig & Thomas Musil noisternig_AT_iem.at; musil_AT_iem.at (c) Institute of Electronic Music and Acoustics, Graz Austria 2005 */
#ifdef NT #pragma warning( disable : 4244 ) #pragma warning( disable : 4305 ) #endif
#include "m_pd.h" #include "iemlib.h" #include <math.h> #include <stdio.h> #include <string.h>
/* ----------------------- n_CLNLMS~ ------------------------------ */ /* -- multiple Constraint LEAKY Normalized Least Mean Square (linear adaptive FIR-filter) -- */
//* -- first input: reference signal -- */ /* -- second input: desired signal -- */ /* -- -- */
/* for further information on adaptive filter design we refer to */ /* [1] Haykin, "Adaptive Filter Theory", 4th ed, Prentice Hall */ /* [2] Benesty, "Adaptive Signal Processing", Springer */
typedef struct n_CLNLMS_tilde_kern { t_symbol *x_w_array_sym_name; t_float *x_w_array_mem_beg; t_float *x_in_ptr_beg;// memory: sig-in vector t_float *x_out_ptr_beg;// memory: sig-out vector t_float *x_in_hist;// start point double buffer for sig-in history } t_n_CLNLMS_tilde_kern;
typedef struct n_CLNLMS_tilde { t_object x_obj; t_n_CLNLMS_tilde_kern *x_my_kern; t_float *x_des_in_ptr_beg;// memory: desired-in vector t_float *x_err_out_ptr_beg;// memory: error-out vector t_int x_n_io;// number of in-channels and filtered out-channels t_int x_rw_index;// read-write-index t_int x_n_order;// filter order t_int x_update;// rounded by 2^n, yields downsampling of learn-rate t_float x_beta;// learn rate [0 .. 2] t_float x_gamma;// normalization t_float x_kappa;// constreint: treshold of energy (clipping) t_float x_leakage;// leakage-Faktor for NLMS t_outlet *x_out_compressing_bang; t_clock *x_clock; t_float x_msi; } t_n_CLNLMS_tilde;
t_class *n_CLNLMS_tilde_class;
static void n_CLNLMS_tilde_tick(t_n_CLNLMS_tilde *x) { outlet_bang(x->x_out_compressing_bang); }
static t_float *n_CLNLMS_tilde_check_array(t_symbol *array_sym_name, t_int length) { t_int n_points; t_garray *a; t_float *vec;
if(!(a = (t_garray *)pd_findbyclass(array_sym_name, garray_class))) { error("%s: no such array for n_CLNLMS~", array_sym_name->s_name); return((t_float *)0); } else if(!garray_getfloatarray(a, &n_points, &vec)) { error("%s: bad template for n_CLNLMS~", array_sym_name->s_name); return((t_float *)0); } else if(n_points < length) { error("%s: bad array-size for n_CLNLMS~: %d", array_sym_name->s_name, n_points); return((t_float *)0); } else { return(vec); } }
static void n_CLNLMS_tilde_beta(t_n_CLNLMS_tilde *x, t_floatarg f) // learn rate { if(f < 0.0f) f = 0.0f; if(f > 2.0f) f = 2.0f;
x->x_beta = f; }
static void n_CLNLMS_tilde_gamma(t_n_CLNLMS_tilde *x, t_floatarg f) // regularization (dither) { if(f < 0.0f) f = 0.0f; if(f > 1.0f) f = 1.0f;
x->x_gamma = f; }
static void n_CLNLMS_tilde_kappa(t_n_CLNLMS_tilde *x, t_floatarg f) // threshold for w_coeff { if(f < 0.0001f) f = 0.0001f; if(f > 10000.0f) f = 10000.0f;
x->x_kappa = f; }
static void n_CLNLMS_tilde_leakage(t_n_CLNLMS_tilde *x, t_floatarg f) // leakage of NLMS { if(f < 0.0001f) f = 0.0001f; if(f > 1.0f) f = 1.0f;
x->x_leakage = f; }
static void n_CLNLMS_tilde_update(t_n_CLNLMS_tilde *x, t_floatarg f) // downsample learn rate { t_int i=1, u = (t_int)f;
if(u < 0) u = 0; else { while(i <= u) // convert u for 2^N i *= 2; // round down i /= 2; u = i; } x->x_update = u; }
/* ============== DSP ======================= */
static t_int *n_CLNLMS_tilde_perform_zero(t_int *w) { t_n_CLNLMS_tilde *x = (t_n_CLNLMS_tilde *)(w[1]); t_int n = (t_int)(w[2]);
t_int n_io = x->x_n_io; t_float *out; t_int i, j;
out = x->x_err_out_ptr_beg; for(i=0; i<n; i++) *out++ = 0.0f; for(j=0; j<n_io; j++) { out = x->x_my_kern[j].x_out_ptr_beg; for(i=0; i<n; i++) *out++ = 0.0f; } return (w+3); }
static t_int *n_CLNLMS_tilde_perform(t_int *w) { t_n_CLNLMS_tilde *x = (t_n_CLNLMS_tilde *)(w[1]); t_int n = (t_int)(w[2]); t_int n_order = x->x_n_order; /* number of filter-order */ t_int rw_index2, rw_index = x->x_rw_index; t_int n_io = x->x_n_io; t_float *in;// first sig in t_float din;// second sig in t_float *filt_out;// first sig out t_float *err_out, err_sum;// second sig out t_float *read_in_hist; t_float *w_filt_coeff; t_float my, my_err, sum; t_float beta = x->x_beta; t_float hgamma, gamma = x->x_gamma; t_float hkappa, kappa = x->x_kappa; t_float hleakage, leakage = x->x_leakage; t_int i, j, k, update_counter; t_int update = x->x_update; t_int ord8=n_order&0xfffffff8; t_int ord_residual=n_order&0x7; t_int compressed = 0;
for(k=0; k<n_io; k++) { if(!x->x_my_kern[k].x_w_array_mem_beg) goto n_CLNLMS_tildeperfzero;// this is Musil/Miller style }
hgamma = gamma * gamma * (float)n_order; //hkappa = kappa * kappa * (float)n_order; hkappa = kappa; // kappa regards to energy value, else use line above
for(i=0, update_counter=0; i<n; i++)// history and (block-)convolution { rw_index2 = rw_index + n_order;
for(k=0; k<n_io; k++)// times n_io { x->x_my_kern[k].x_in_hist[rw_index] = x->x_my_kern[k].x_in_ptr_beg[i]; // save inputs into variabel & history x->x_my_kern[k].x_in_hist[rw_index+n_order] = x->x_my_kern[k].x_in_ptr_beg[i]; } din = x->x_des_in_ptr_beg[i];
// begin convolution err_sum = din; for(k=0; k<n_io; k++)// times n_io { sum = 0.0f; w_filt_coeff = x->x_my_kern[k].x_w_array_mem_beg; // Musil's special convolution buffer struct read_in_hist = &x->x_my_kern[k].x_in_hist[rw_index2]; for(j=0; j<ord8; j+=8) // loop unroll 8 taps { sum += w_filt_coeff[0] * read_in_hist[0]; sum += w_filt_coeff[1] * read_in_hist[-1]; sum += w_filt_coeff[2] * read_in_hist[-2]; sum += w_filt_coeff[3] * read_in_hist[-3]; sum += w_filt_coeff[4] * read_in_hist[-4]; sum += w_filt_coeff[5] * read_in_hist[-5]; sum += w_filt_coeff[6] * read_in_hist[-6]; sum += w_filt_coeff[7] * read_in_hist[-7]; w_filt_coeff += 8; read_in_hist -= 8; } for(j=0; j<ord_residual; j++) // for filter order < 2^N sum += w_filt_coeff[j] * read_in_hist[-j];
x->x_my_kern[k].x_out_ptr_beg[i] = sum; err_sum -= sum; } x->x_err_out_ptr_beg[i] = err_sum; // end convolution
if(update) // downsampling of learn rate { update_counter++; if(update_counter >= update) { update_counter = 0;
for(k=0; k<n_io; k++)// times n_io { sum = 0.0f;// calculate energy for last n-order samples in filter read_in_hist = &x->x_my_kern[k].x_in_hist[rw_index2]; for(j=0; j<ord8; j+=8) // unrolling quadrature calc { sum += read_in_hist[0] * read_in_hist[0]; sum += read_in_hist[-1] * read_in_hist[-1]; sum += read_in_hist[-2] * read_in_hist[-2]; sum += read_in_hist[-3] * read_in_hist[-3]; sum += read_in_hist[-4] * read_in_hist[-4]; sum += read_in_hist[-5] * read_in_hist[-5]; sum += read_in_hist[-6] * read_in_hist[-6]; sum += read_in_hist[-7] * read_in_hist[-7]; read_in_hist -= 8; } for(j=0; j<ord_residual; j++) // residual sum += read_in_hist[-j] * read_in_hist[-j]; // [-j] only valid for Musil's double buffer structure sum += hgamma; // convert gamma corresponding to filter order my = beta / sum;// calculate mue
my_err = my * err_sum; w_filt_coeff = x->x_my_kern[k].x_w_array_mem_beg; read_in_hist = &x->x_my_kern[k].x_in_hist[rw_index2]; sum = 0.0f; for(j=0; j<ord8; j+=8) // unrolling quadrature calc { w_filt_coeff[0] = leakage * w_filt_coeff[0] + read_in_hist[0] * my_err; sum += w_filt_coeff[0] * w_filt_coeff[0]; w_filt_coeff[1] = leakage * w_filt_coeff[1] + read_in_hist[-1] * my_err; sum += w_filt_coeff[1] * w_filt_coeff[1]; w_filt_coeff[2] = leakage * w_filt_coeff[2] + read_in_hist[-2] * my_err; sum += w_filt_coeff[2] * w_filt_coeff[2]; w_filt_coeff[3] = leakage * w_filt_coeff[3] + read_in_hist[-3] * my_err; sum += w_filt_coeff[3] * w_filt_coeff[3]; w_filt_coeff[4] = leakage * w_filt_coeff[4] + read_in_hist[-4] * my_err; sum += w_filt_coeff[4] * w_filt_coeff[4]; w_filt_coeff[5] = leakage * w_filt_coeff[5] + read_in_hist[-5] * my_err; sum += w_filt_coeff[5] * w_filt_coeff[5]; w_filt_coeff[6] = leakage * w_filt_coeff[6] + read_in_hist[-6] * my_err; sum += w_filt_coeff[6] * w_filt_coeff[6]; w_filt_coeff[7] = leakage * w_filt_coeff[7] + read_in_hist[-7] * my_err; sum += w_filt_coeff[7] * w_filt_coeff[7]; w_filt_coeff += 8; read_in_hist -= 8; } for(j=0; j<ord_residual; j++) // residual { w_filt_coeff[j] = leakage * w_filt_coeff[j] + read_in_hist[-j] * my_err; sum += w_filt_coeff[j] * w_filt_coeff[j]; } if(sum > hkappa) { compressed = 1; my = sqrt(hkappa/sum); w_filt_coeff = x->x_my_kern[k].x_w_array_mem_beg; for(j=0; j<ord8; j+=8) // unrolling quadrature calc { w_filt_coeff[0] *= my; w_filt_coeff[1] *= my; w_filt_coeff[2] *= my; w_filt_coeff[3] *= my; w_filt_coeff[4] *= my; w_filt_coeff[5] *= my; w_filt_coeff[6] *= my; w_filt_coeff[7] *= my; w_filt_coeff += 8; } for(j=0; j<ord_residual; j++) // residual w_filt_coeff[j] *= my; } } } } rw_index++; if(rw_index >= n_order) rw_index -= n_order; }
x->x_rw_index = rw_index; // wieder in die garage stellen
if(compressed) clock_delay(x->x_clock, 0);
return(w+3);
n_CLNLMS_tildeperfzero:
err_out = x->x_err_out_ptr_beg; for(i=0; i<n; i++) *err_out++ = 0.0f; for(j=0; j<n_io; j++) { filt_out = x->x_my_kern[j].x_out_ptr_beg; for(i=0; i<n; i++) *filt_out++ = 0.0f; }
return(w+3); }
static void n_CLNLMS_tilde_dsp(t_n_CLNLMS_tilde *x, t_signal **sp) { t_int i, n = sp[0]->s_n; t_int ok_w = 1; t_int m = x->x_n_io;
for(i=0; i<m; i++) x->x_my_kern[i].x_in_ptr_beg = sp[i]->s_vec; x->x_des_in_ptr_beg = sp[m]->s_vec; for(i=0; i<m; i++) x->x_my_kern[i].x_out_ptr_beg = sp[i+m+1]->s_vec; x->x_err_out_ptr_beg = sp[2*m+1]->s_vec;
for(i=0; i<m; i++) { x->x_my_kern[i].x_w_array_mem_beg = n_CLNLMS_tilde_check_array(x->x_my_kern[i].x_w_array_sym_name, x->x_n_order); if(!x->x_my_kern[i].x_w_array_mem_beg) ok_w = 0; }
if(!ok_w) dsp_add(n_CLNLMS_tilde_perform_zero, 2, x, n); else dsp_add(n_CLNLMS_tilde_perform, 2, x, n); }
/* setup/setdown things */
static void n_CLNLMS_tilde_free(t_n_CLNLMS_tilde *x) { t_int i, n_io=x->x_n_io, n_order=x->x_n_order;
for(i=0; i<n_io; i++) freebytes(x->x_my_kern[i].x_in_hist, 2*x->x_n_order*sizeof(t_float)); freebytes(x->x_my_kern, n_io*sizeof(t_n_CLNLMS_tilde_kern));
clock_free(x->x_clock); }
static void *n_CLNLMS_tilde_new(t_symbol *s, t_int argc, t_atom *argv) { t_n_CLNLMS_tilde *x = (t_n_CLNLMS_tilde *)pd_new(n_CLNLMS_tilde_class); char buffer[400]; t_int i, n_order=39, n_io=1; t_symbol *w_name; t_float beta=0.1f; t_float gamma=0.00001f; t_float kappa = 1.0f; t_float leakage = 0.99f;
if((argc >= 7) && IS_A_FLOAT(argv,0) && //IS_A_FLOAT/SYMBOL from iemlib.h IS_A_FLOAT(argv,1) && IS_A_FLOAT(argv,2) && IS_A_FLOAT(argv,3) && IS_A_FLOAT(argv,4) && IS_A_FLOAT(argv,5) && IS_A_SYMBOL(argv,6)) { n_io = (t_int)atom_getintarg(0, argc, argv); n_order = (t_int)atom_getintarg(1, argc, argv); beta = (t_float)atom_getfloatarg(2, argc, argv); gamma = (t_float)atom_getfloatarg(3, argc, argv); kappa = (t_float)atom_getfloatarg(4, argc, argv); leakage = (t_float)atom_getfloatarg(5, argc, argv); w_name = (t_symbol *)atom_getsymbolarg(6, argc, argv);
if(beta < 0.0f) beta = 0.0f; if(beta > 2.0f) beta = 2.0f;
if(gamma < 0.0f) gamma = 0.0f; if(gamma > 1.0f) gamma = 1.0f;
if(kappa < 0.0001f) kappa = 0.0001f; if(kappa > 10000.0f) kappa = 10000.0f;
if(leakage < 0.0001f) leakage = 0.0001f; if(leakage > 1.0f) leakage = 1.0f;
if(n_order < 2) n_order = 2; if(n_order > 11111) n_order = 11111;
if(n_io < 1) n_io = 1; if(n_io > 60) n_io = 60;
for(i=0; i<n_io; i++) inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_signal, &s_signal); for(i=0; i<=n_io; i++) outlet_new(&x->x_obj, &s_signal);
x->x_out_compressing_bang = outlet_new(&x->x_obj, &s_bang);
x->x_msi = 0; x->x_n_io = n_io; x->x_n_order = n_order; x->x_update = 0; x->x_beta = beta; x->x_gamma = gamma; x->x_kappa = kappa; x->x_leakage = leakage; x->x_my_kern = (t_n_CLNLMS_tilde_kern *)getbytes(x->x_n_io*sizeof(t_n_CLNLMS_tilde_kern)); for(i=0; i<n_io; i++) { sprintf(buffer, "%d_%s", i+1, w_name->s_name); x->x_my_kern[i].x_w_array_sym_name = gensym(buffer); x->x_my_kern[i].x_w_array_mem_beg = (t_float *)0; x->x_my_kern[i].x_in_hist = (t_float *)getbytes(2*x->x_n_order*sizeof(t_float)); } x->x_clock = clock_new(x, (t_method)n_CLNLMS_tilde_tick);
return(x); } else { post("n_CLNLMSC~-ERROR: need 6 float- + 1 symbol-arguments:"); post(" number_of_filters + order_of_filters + learnrate_beta + security_value_gamma + threshold_kappa + leakage_factor_lambda + array_name_taps"); return(0); } }
void n_CLNLMS_tilde_setup(void) { n_CLNLMS_tilde_class = class_new(gensym("n_CLNLMS~"), (t_newmethod)n_CLNLMS_tilde_new, (t_method)n_CLNLMS_tilde_free, sizeof(t_n_CLNLMS_tilde), 0, A_GIMME, 0); CLASS_MAINSIGNALIN(n_CLNLMS_tilde_class, t_n_CLNLMS_tilde, x_msi); class_addmethod(n_CLNLMS_tilde_class, (t_method)n_CLNLMS_tilde_dsp, gensym("dsp"), 0); class_addmethod(n_CLNLMS_tilde_class, (t_method)n_CLNLMS_tilde_update, gensym("update"), A_FLOAT, 0); // method: downsampling factor of learning (multiple of 2^N) class_addmethod(n_CLNLMS_tilde_class, (t_method)n_CLNLMS_tilde_beta, gensym("beta"), A_FLOAT, 0); //method: normalized learning rate class_addmethod(n_CLNLMS_tilde_class, (t_method)n_CLNLMS_tilde_gamma, gensym("gamma"), A_FLOAT, 0); // method: dithering noise related to signal class_addmethod(n_CLNLMS_tilde_class, (t_method)n_CLNLMS_tilde_kappa, gensym("kappa"), A_FLOAT, 0); // method: threshold for compressing w_coeff class_addmethod(n_CLNLMS_tilde_class, (t_method)n_CLNLMS_tilde_leakage, gensym("leakage"), A_FLOAT, 0); // method: leakage factor [0 1] for w update
//class_sethelpsymbol(n_CLNLMS_tilde_class, gensym("iemhelp2/n_CLNLMS~")); }
Index: makefile_win =================================================================== RCS file: /cvsroot/pure-data/externals/iem/iem_adaptfilt/src/makefile_win,v retrieving revision 1.1 retrieving revision 1.2 diff -C2 -d -r1.1 -r1.2 *** makefile_win 2 Aug 2006 14:02:28 -0000 1.1 --- makefile_win 11 Jan 2007 17:40:10 -0000 1.2 *************** *** 22,29 ****
! SRC = sigNLMS.c \ ! sigNLMSCC.c \ ! sign_CNLMS.c \ ! sign_CLNLMS.c \ iem_adaptfilt.c
--- 22,29 ----
! SRC = NLMS~.c \ ! NLMSCC~.c \ ! n_CNLMS~.c \ ! n_CLNLMS~.c \ iem_adaptfilt.c
--- NEW FILE: NLMS~.c --- /* For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution.
NLMS normalized least mean square (LMS) algorithm lib iem_adaptfilt written by Markus Noisternig & Thomas Musil noisternig_AT_iem.at; musil_AT_iem.at (c) Institute of Electronic Music and Acoustics, Graz Austria 2005 */
#ifdef NT #pragma warning( disable : 4244 ) #pragma warning( disable : 4305 ) #endif
#include "m_pd.h" #include "iemlib.h" #include <math.h> #include <stdio.h> #include <string.h>
/* ----------------------- NLMS~ ------------------------------ */ /* -- Normalized Least Mean Square (linear adaptive FIR-filter) -- */ /* -- first input: reference signal -- */ /* -- second input: desired signal -- */ /* -- -- */
/* for further information on adaptive filter design we refer to */ /* [1] Haykin, "Adaptive Filter Theory", 4th ed, Prentice Hall */ /* [2] Benesty, "Adaptive Signal Processing", Springer */
typedef struct NLMS_tilde { t_object x_obj; t_symbol *x_w_array_sym_name; t_float *x_w_array_mem_beg; t_float *x_io_ptr_beg[4];// memory: 2 sig-in and 2 sig-out vectors t_float *x_in_hist;// start point double buffer for sig-in history t_int x_rw_index;// read-write-index t_int x_n_order;// order of filter t_int x_update;// 2^n rounded value, downsampling of update speed t_float x_beta;// learn rate [0 .. 2] t_float x_gamma;// regularization t_float x_msi; } t_NLMS_tilde;
t_class *NLMS_tilde_class;
static t_float *NLMS_tilde_check_array(t_symbol *array_sym_name, t_int length) { t_int n_points; t_garray *a; t_float *vec;
if(!(a = (t_garray *)pd_findbyclass(array_sym_name, garray_class))) { error("%s: no such array for NLMS~", array_sym_name->s_name); return((t_float *)0); } else if(!garray_getfloatarray(a, &n_points, &vec)) { error("%s: bad template for NLMS~", array_sym_name->s_name); return((t_float *)0); } else if(n_points < length) { error("%s: bad array-size for NLMS~: %d", array_sym_name->s_name, n_points); return((t_float *)0); } else { return(vec); } }
static void NLMS_tilde_beta(t_NLMS_tilde *x, t_floatarg f) // learn rate { if(f < 0.0f) f = 0.0f; if(f > 2.0f) f = 2.0f;
x->x_beta = f; }
static void NLMS_tilde_gamma(t_NLMS_tilde *x, t_floatarg f) // regularization factor (dither) { if(f < 0.0f) f = 0.0f; if(f > 1.0f) f = 1.0f;
x->x_gamma = f; }
static void NLMS_tilde_update(t_NLMS_tilde *x, t_floatarg f) // downsample learn-rate { t_int i=1, u = (t_int)f;
if(u < 0) u = 0; else { while(i <= u) // convert u for 2^N i *= 2; // round downwards i /= 2; u = i; } x->x_update = u; }
/* ============== DSP ======================= */
static t_int *NLMS_tilde_perform_zero(t_int *w) { t_NLMS_tilde *x = (t_NLMS_tilde *)(w[1]); t_int n = (t_int)(w[2]);
t_float **io = x->x_io_ptr_beg; t_float *out; t_int i, j;
for(j=0; j<2; j++)/* output-vector-row */ { out = io[j+2]; for(i=0; i<n; i++) { *out++ = 0.0f; } } return (w+3); }
static t_int *NLMS_tilde_perform(t_int *w) { t_NLMS_tilde *x = (t_NLMS_tilde *)(w[1]); t_int n = (t_int)(w[2]); t_int n_order = x->x_n_order; /* number of filter-order */ t_int rw_index = x->x_rw_index; t_float *in = x->x_io_ptr_beg[0];// first sig in t_float *desired_in = x->x_io_ptr_beg[1], din;// second sig in t_float *filt_out = x->x_io_ptr_beg[2];// first sig out t_float *err_out = x->x_io_ptr_beg[3], eout;// second sig out t_float *write_in_hist1 = x->x_in_hist; t_float *write_in_hist2 = write_in_hist1+n_order; t_float *read_in_hist = write_in_hist2; t_float *w_filt_coeff = x->x_w_array_mem_beg; t_float my, my_err, sum; t_float beta = x->x_beta; t_float gamma = x->x_gamma; t_int i, j, update_counter; t_int update = x->x_update; t_int ord8=n_order&0xfffffff8; t_int ord_residual=n_order&0x7;
if(!w_filt_coeff) goto NLMS_tildeperfzero;// this is quick&dirty Musil/Miller style
for(i=0, update_counter=0; i<n; i++)// store history and convolve { write_in_hist1[rw_index] = in[i]; // save inputs to variable & history write_in_hist2[rw_index] = in[i]; din = desired_in[i];
// begin convolution sum = 0.0f; w_filt_coeff = x->x_w_array_mem_beg; // Musil's special convolution buffer struct read_in_hist = &write_in_hist2[rw_index]; for(j=0; j<ord8; j+=8) // loop unroll 8 taps { sum += w_filt_coeff[0] * read_in_hist[0]; sum += w_filt_coeff[1] * read_in_hist[-1]; sum += w_filt_coeff[2] * read_in_hist[-2]; sum += w_filt_coeff[3] * read_in_hist[-3]; sum += w_filt_coeff[4] * read_in_hist[-4]; sum += w_filt_coeff[5] * read_in_hist[-5]; sum += w_filt_coeff[6] * read_in_hist[-6]; sum += w_filt_coeff[7] * read_in_hist[-7]; w_filt_coeff += 8; read_in_hist -= 8; } for(j=0; j<ord_residual; j++) // for filter order < 2^N sum += w_filt_coeff[j] * read_in_hist[-j];
filt_out[i] = sum; eout = din - filt_out[i]; // buffer-struct for further use err_out[i] = eout;
if(update) // downsampling for learn rate { update_counter++; if(update_counter >= update) { update_counter = 0;
sum = 0.0f;// calculate energy for last n-order samples in filter read_in_hist = &write_in_hist2[rw_index]; for(j=0; j<ord8; j+=8) // unrolling quadrature calc { sum += read_in_hist[0] * read_in_hist[0]; sum += read_in_hist[-1] * read_in_hist[-1]; sum += read_in_hist[-2] * read_in_hist[-2]; sum += read_in_hist[-3] * read_in_hist[-3]; sum += read_in_hist[-4] * read_in_hist[-4]; sum += read_in_hist[-5] * read_in_hist[-5]; sum += read_in_hist[-6] * read_in_hist[-6]; sum += read_in_hist[-7] * read_in_hist[-7]; read_in_hist -= 8; } for(j=0; j<ord_residual; j++) // residual sum += read_in_hist[-j] * read_in_hist[-j]; // [-j] only valid for Musil's double buffer structure sum += gamma * gamma * (float)n_order; // convert gamma corresponding to filter order my = beta / sum;// calculate mue
my_err = my * eout; w_filt_coeff = x->x_w_array_mem_beg; // coefficient constraints read_in_hist = &write_in_hist2[rw_index]; for(j=0; j<n_order; j++) // without unroll w_filt_coeff[j] += read_in_hist[-j] * my_err; } } rw_index++; if(rw_index >= n_order) rw_index -= n_order; }
x->x_rw_index = rw_index; // back to start
return(w+3);
NLMS_tildeperfzero:
while(n--) { *filt_out++ = 0.0f; *err_out++ = 0.0f; } return(w+3); }
static void NLMS_tilde_dsp(t_NLMS_tilde *x, t_signal **sp) { t_int i, n = sp[0]->s_n;
for(i=0; i<4; i++) // store io_vec x->x_io_ptr_beg[i] = sp[i]->s_vec;
x->x_w_array_mem_beg = NLMS_tilde_check_array(x->x_w_array_sym_name, x->x_n_order);
if(!x->x_w_array_mem_beg) dsp_add(NLMS_tilde_perform_zero, 2, x, n); else dsp_add(NLMS_tilde_perform, 2, x, n); }
/* setup/setdown things */
static void NLMS_tilde_free(t_NLMS_tilde *x) { freebytes(x->x_in_hist, 2*x->x_n_order*sizeof(t_float)); }
static void *NLMS_tilde_new(t_symbol *s, t_int argc, t_atom *argv) { t_NLMS_tilde *x = (t_NLMS_tilde *)pd_new(NLMS_tilde_class); t_int i, n_order=39; t_symbol *w_name; t_float beta=0.1f; t_float gamma=0.00001f;
if((argc >= 4) && IS_A_FLOAT(argv,0) && //IS_A_FLOAT/SYMBOL from iemlib.h IS_A_FLOAT(argv,1) && IS_A_FLOAT(argv,2) && IS_A_SYMBOL(argv,3)) { n_order = (t_int)atom_getintarg(0, argc, argv); beta = (t_float)atom_getfloatarg(1, argc, argv); gamma = (t_float)atom_getfloatarg(2, argc, argv); w_name = (t_symbol *)atom_getsymbolarg(3, argc, argv);
if(beta < 0.0f) beta = 0.0f; if(beta > 2.0f) beta = 2.0f;
if(gamma < 0.0f) gamma = 0.0f; if(gamma > 1.0f) gamma = 1.0f;
if(n_order < 2) n_order = 2; if(n_order > 11111) n_order = 11111;
inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_signal, &s_signal); outlet_new(&x->x_obj, &s_signal); outlet_new(&x->x_obj, &s_signal);
x->x_msi = 0; x->x_n_order = n_order; x->x_update = 0; x->x_beta = beta; x->x_gamma = gamma; // 2 times in and one time desired_in memory allocation (history) x->x_in_hist = (t_float *)getbytes(2*x->x_n_order*sizeof(t_float));
// table-symbols will be linked to their memory in future (dsp_routine) x->x_w_array_sym_name = gensym(w_name->s_name); x->x_w_array_mem_beg = (t_float *)0;
return(x); } else { post("NLMS~-ERROR: need 3 float- + 1 symbol-arguments:"); post(" order_of_filter + learnrate_beta + security_value + array_name_taps"); return(0); } }
void NLMS_tilde_setup(void) { NLMS_tilde_class = class_new(gensym("NLMS~"), (t_newmethod)NLMS_tilde_new, (t_method)NLMS_tilde_free, sizeof(t_NLMS_tilde), 0, A_GIMME, 0); CLASS_MAINSIGNALIN(NLMS_tilde_class, t_NLMS_tilde, x_msi); class_addmethod(NLMS_tilde_class, (t_method)NLMS_tilde_dsp, gensym("dsp"), 0); class_addmethod(NLMS_tilde_class, (t_method)NLMS_tilde_update, gensym("update"), A_FLOAT, 0); // method: downsampling factor of learning (multiple of 2^N) class_addmethod(NLMS_tilde_class, (t_method)NLMS_tilde_beta, gensym("beta"), A_FLOAT, 0); //method: normalized learning rate class_addmethod(NLMS_tilde_class, (t_method)NLMS_tilde_gamma, gensym("gamma"), A_FLOAT, 0); // method: dithering noise related to signal
//class_sethelpsymbol(NLMS_tilde_class, gensym("iemhelp2/NLMS~")); }
--- NEW FILE: NLMSCC~.c --- /* For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution.
NLMSCC normalized LMS algorithm with coefficient constraints lib iem_adaptfilt written by Markus Noisternig & Thomas Musil noisternig_AT_iem.at; musil_AT_iem.at (c) Institute of Electronic Music and Acoustics, Graz Austria 2005 */
#ifdef NT #pragma warning( disable : 4244 ) #pragma warning( disable : 4305 ) #endif
#include "m_pd.h" #include "iemlib.h" #include <math.h> #include <stdio.h> #include <string.h>
/* ----------------------- NLMSCC~ ------------------------------ */ /* -- Normalized Least Mean Square (linear adaptive FIR-filter) -- */ /* -- with Coefficient Constraint /* -- first input: reference signal -- */ /* -- second input: desired signal -- */ /* -- -- */ /* for further information on adaptive filter design we refer to */ /* [1] Haykin, "Adaptive Filter Theory", 4th ed, Prentice Hall */ /* [2] Benesty, "Adaptive Signal Processing", Springer */ /* */
typedef struct NLMSCC_tilde { t_object x_obj; t_symbol *x_w_array_sym_name; t_float *x_w_array_mem_beg; t_symbol *x_wmin_array_sym_name; t_float *x_wmin_array_mem_beg; t_symbol *x_wmax_array_sym_name; t_float *x_wmax_array_mem_beg; t_float *x_io_ptr_beg[4];// memory: 2 sig-in and 2 sig-out vectors t_float *x_in_hist;// start point double buffer for sig-in history t_int x_rw_index;// read-write-index t_int x_n_order;// order of filter t_int x_update;// 2^n rounded value, downsampling of update speed t_float x_beta;// learn rate [0 .. 2] t_float x_gamma;// regularization t_outlet *x_out_clipping_bang; t_clock *x_clock; t_float x_msi; } t_NLMSCC_tilde;
t_class *NLMSCC_tilde_class;
static void NLMSCC_tilde_tick(t_NLMSCC_tilde *x) { outlet_bang(x->x_out_clipping_bang); }
static t_float *NLMSCC_tilde_check_array(t_symbol *array_sym_name, t_int length) { t_int n_points; t_garray *a; t_float *vec;
if(!(a = (t_garray *)pd_findbyclass(array_sym_name, garray_class))) { error("%s: no such array for NLMSCC~", array_sym_name->s_name); return((t_float *)0); } else if(!garray_getfloatarray(a, &n_points, &vec)) { error("%s: bad template for NLMSCC~", array_sym_name->s_name); return((t_float *)0); } else if(n_points < length) { error("%s: bad array-size for NLMSCC~: %d", array_sym_name->s_name, n_points); return((t_float *)0); } else { return(vec); } }
static void NLMSCC_tilde_beta(t_NLMSCC_tilde *x, t_floatarg f) // learn rate { if(f < 0.0f) f = 0.0f; if(f > 2.0f) f = 2.0f;
x->x_beta = f; }
static void NLMSCC_tilde_gamma(t_NLMSCC_tilde *x, t_floatarg f) // regularization factor (dither) { if(f < 0.0f) f = 0.0f; if(f > 1.0f) f = 1.0f;
x->x_gamma = f; }
static void NLMSCC_tilde_update(t_NLMSCC_tilde *x, t_floatarg f) // downsample of learn-rate { t_int i=1, u = (t_int)f;
if(u < 0) u = 0; else { while(i <= u) // convert u for 2^N i *= 2; // round downwards i /= 2; u = i; } x->x_update = u; }
/* ============== DSP ======================= */
static t_int *NLMSCC_tilde_perform_zero(t_int *w) { t_NLMSCC_tilde *x = (t_NLMSCC_tilde *)(w[1]); t_int n = (t_int)(w[2]);
t_float **io = x->x_io_ptr_beg; t_float *out; t_int i, j;
for(j=0; j<2; j++)/* output-vector-row */ { out = io[j+2]; for(i=0; i<n; i++) { *out++ = 0.0f; } } return (w+3); }
static t_int *NLMSCC_tilde_perform(t_int *w) { t_NLMSCC_tilde *x = (t_NLMSCC_tilde *)(w[1]); t_int n = (t_int)(w[2]); t_int n_order = x->x_n_order; /* filter-order */ t_int rw_index = x->x_rw_index; t_float *in = x->x_io_ptr_beg[0];// first sig in t_float *desired_in = x->x_io_ptr_beg[1], din;// second sig in t_float *filt_out = x->x_io_ptr_beg[2];// first sig out t_float *err_out = x->x_io_ptr_beg[3], eout;// second sig out t_float *write_in_hist1 = x->x_in_hist; t_float *write_in_hist2 = write_in_hist1+n_order; t_float *read_in_hist = write_in_hist2; t_float *w_filt_coeff = x->x_w_array_mem_beg; t_float *wmin_filt_coeff = x->x_wmin_array_mem_beg; t_float *wmax_filt_coeff = x->x_wmax_array_mem_beg; t_float my, my_err, sum; t_float beta = x->x_beta; t_float gamma = x->x_gamma; t_int i, j, update_counter; t_int update = x->x_update; t_int ord8=n_order&0xfffffff8; t_int ord_residual=n_order&0x7; t_int clipped = 0;
if(!w_filt_coeff) goto NLMSCC_tildeperfzero;// this is Musil/Miller style if(!wmin_filt_coeff) goto NLMSCC_tildeperfzero; if(!wmax_filt_coeff) goto NLMSCC_tildeperfzero;// if not constrained, perform zero
for(i=0, update_counter=0; i<n; i++)// store in history and convolve { write_in_hist1[rw_index] = in[i]; // save inputs into variabel & history write_in_hist2[rw_index] = in[i]; din = desired_in[i];
// begin convolution sum = 0.0f; w_filt_coeff = x->x_w_array_mem_beg; // Musil's special convolution buffer struct read_in_hist = &write_in_hist2[rw_index]; for(j=0; j<ord8; j+=8) // loop unroll 8 taps { sum += w_filt_coeff[0] * read_in_hist[0]; sum += w_filt_coeff[1] * read_in_hist[-1]; sum += w_filt_coeff[2] * read_in_hist[-2]; sum += w_filt_coeff[3] * read_in_hist[-3]; sum += w_filt_coeff[4] * read_in_hist[-4]; sum += w_filt_coeff[5] * read_in_hist[-5]; sum += w_filt_coeff[6] * read_in_hist[-6]; sum += w_filt_coeff[7] * read_in_hist[-7]; w_filt_coeff += 8; read_in_hist -= 8; } for(j=0; j<ord_residual; j++) // for filter order < 2^N sum += w_filt_coeff[j] * read_in_hist[-j];
filt_out[i] = sum; eout = din - filt_out[i]; // buffer-struct for further use err_out[i] = eout;
if(update) // downsampling for learn rate { update_counter++; if(update_counter >= update) { update_counter = 0;
sum = 0.0f;// calculate energy for last n-order samples in filter read_in_hist = &write_in_hist2[rw_index]; for(j=0; j<ord8; j+=8) // unrolling quadrature calc { sum += read_in_hist[0] * read_in_hist[0]; sum += read_in_hist[-1] * read_in_hist[-1]; sum += read_in_hist[-2] * read_in_hist[-2]; sum += read_in_hist[-3] * read_in_hist[-3]; sum += read_in_hist[-4] * read_in_hist[-4]; sum += read_in_hist[-5] * read_in_hist[-5]; sum += read_in_hist[-6] * read_in_hist[-6]; sum += read_in_hist[-7] * read_in_hist[-7]; read_in_hist -= 8; } for(j=0; j<ord_residual; j++) // residual sum += read_in_hist[-j] * read_in_hist[-j]; // [-j] only valid for Musil's double buffer structure sum += gamma * gamma * (float)n_order; // convert gamma corresponding to filter order my = beta / sum;// calculate mue
my_err = my * eout; w_filt_coeff = x->x_w_array_mem_beg; // coefficient constraints wmin_filt_coeff = x->x_wmin_array_mem_beg; wmax_filt_coeff = x->x_wmax_array_mem_beg; read_in_hist = &write_in_hist2[rw_index]; for(j=0; j<n_order; j++) // without unroll { w_filt_coeff[j] += read_in_hist[-j] * my_err; if(w_filt_coeff[j] > wmax_filt_coeff[j]) { w_filt_coeff[j] = wmax_filt_coeff[j]; clipped = 1; } else if(w_filt_coeff[j] < wmin_filt_coeff[j]) { w_filt_coeff[j] = wmin_filt_coeff[j]; clipped = 1; } } } } rw_index++; if(rw_index >= n_order) rw_index -= n_order; }
x->x_rw_index = rw_index; // back to start
if(clipped) clock_delay(x->x_clock, 0); return(w+3);
NLMSCC_tildeperfzero:
while(n--) { *filt_out++ = 0.0f; *err_out++ = 0.0f; } return(w+3); }
static void NLMSCC_tilde_dsp(t_NLMSCC_tilde *x, t_signal **sp) { t_int i, n = sp[0]->s_n;
for(i=0; i<4; i++) // store io_vec x->x_io_ptr_beg[i] = sp[i]->s_vec;
x->x_w_array_mem_beg = NLMSCC_tilde_check_array(x->x_w_array_sym_name, x->x_n_order); x->x_wmin_array_mem_beg = NLMSCC_tilde_check_array(x->x_wmin_array_sym_name, x->x_n_order); x->x_wmax_array_mem_beg = NLMSCC_tilde_check_array(x->x_wmax_array_sym_name, x->x_n_order);
if(!(x->x_w_array_mem_beg && x->x_wmin_array_mem_beg && x->x_wmax_array_mem_beg)) dsp_add(NLMSCC_tilde_perform_zero, 2, x, n); else dsp_add(NLMSCC_tilde_perform, 2, x, n); }
/* setup/setdown things */
static void NLMSCC_tilde_free(t_NLMSCC_tilde *x) {
freebytes(x->x_in_hist, 2*x->x_n_order*sizeof(t_float));
clock_free(x->x_clock); }
static void *NLMSCC_tilde_new(t_symbol *s, t_int argc, t_atom *argv) { t_NLMSCC_tilde *x = (t_NLMSCC_tilde *)pd_new(NLMSCC_tilde_class); t_int i, n_order=39; t_symbol *w_name; t_symbol *wmin_name; t_symbol *wmax_name; t_float beta=0.1f; t_float gamma=0.00001f;
if((argc >= 6) && IS_A_FLOAT(argv,0) && //IS_A_FLOAT/SYMBOL from iemlib.h IS_A_FLOAT(argv,1) && IS_A_FLOAT(argv,2) && IS_A_SYMBOL(argv,3) && IS_A_SYMBOL(argv,4) && IS_A_SYMBOL(argv,5)) { n_order = (t_int)atom_getintarg(0, argc, argv); beta = (t_float)atom_getfloatarg(1, argc, argv); gamma = (t_float)atom_getfloatarg(2, argc, argv); w_name = (t_symbol *)atom_getsymbolarg(3, argc, argv); wmin_name = (t_symbol *)atom_getsymbolarg(4, argc, argv); wmax_name = (t_symbol *)atom_getsymbolarg(5, argc, argv);
if(beta < 0.0f) beta = 0.0f; if(beta > 2.0f) beta = 2.0f;
if(gamma < 0.0f) gamma = 0.0f; if(gamma > 1.0f) gamma = 1.0f;
if(n_order < 2) n_order = 2; if(n_order > 11111) n_order = 11111;
inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_signal, &s_signal); outlet_new(&x->x_obj, &s_signal); outlet_new(&x->x_obj, &s_signal); x->x_out_clipping_bang = outlet_new(&x->x_obj, &s_bang);
x->x_msi = 0; x->x_n_order = n_order; x->x_update = 0; x->x_beta = beta; x->x_gamma = gamma; // 2 times in and one time desired_in memory allocation (history) x->x_in_hist = (t_float *)getbytes(2*x->x_n_order*sizeof(t_float));
// table-symbols will be linked to their memory in future (dsp_routine) x->x_w_array_sym_name = gensym(w_name->s_name); x->x_w_array_mem_beg = (t_float *)0; x->x_wmin_array_sym_name = gensym(wmin_name->s_name); x->x_wmin_array_mem_beg = (t_float *)0; x->x_wmax_array_sym_name = gensym(wmax_name->s_name); x->x_wmax_array_mem_beg = (t_float *)0;
x->x_clock = clock_new(x, (t_method)NLMSCC_tilde_tick);
return(x); } else { post("NLMSCC~-ERROR: need 3 float- + 3 symbol-arguments:"); post(" order_of_filter + learnrate_beta + security_value + array_name_taps + array_name_tap_min + array_name_tap_max"); return(0); } }
void NLMSCC_tilde_setup(void) { NLMSCC_tilde_class = class_new(gensym("NLMSCC~"), (t_newmethod)NLMSCC_tilde_new, (t_method)NLMSCC_tilde_free, sizeof(t_NLMSCC_tilde), 0, A_GIMME, 0); CLASS_MAINSIGNALIN(NLMSCC_tilde_class, t_NLMSCC_tilde, x_msi); class_addmethod(NLMSCC_tilde_class, (t_method)NLMSCC_tilde_dsp, gensym("dsp"), 0); class_addmethod(NLMSCC_tilde_class, (t_method)NLMSCC_tilde_update, gensym("update"), A_FLOAT, 0); // method: downsampling factor of learning (multiple of 2^N) class_addmethod(NLMSCC_tilde_class, (t_method)NLMSCC_tilde_beta, gensym("beta"), A_FLOAT, 0); //method: normalized learning rate class_addmethod(NLMSCC_tilde_class, (t_method)NLMSCC_tilde_gamma, gensym("gamma"), A_FLOAT, 0); // method: dithering noise related to signal
//class_sethelpsymbol(NLMSCC_tilde_class, gensym("iemhelp2/NLMSCC~")); }
Index: makefile.txt =================================================================== RCS file: /cvsroot/pure-data/externals/iem/iem_adaptfilt/src/makefile.txt,v retrieving revision 1.1 retrieving revision 1.2 diff -C2 -d -r1.1 -r1.2 *** makefile.txt 2 Aug 2006 14:02:28 -0000 1.1 --- makefile.txt 11 Jan 2007 17:40:10 -0000 1.2 *************** *** 18,24 **** # the sources
! SRC = sigNLMS.c \ ! sigNLMSCC.c \ ! sign_CNLMS.c \ iem_adaptfilt.c
--- 18,24 ---- # the sources
! SRC = NLMS~.c \ ! NLMSCC~.c \ ! n_CNLMS~.c \ iem_adaptfilt.c
Index: iem_adaptfilt.c =================================================================== RCS file: /cvsroot/pure-data/externals/iem/iem_adaptfilt/src/iem_adaptfilt.c,v retrieving revision 1.1 retrieving revision 1.2 diff -C2 -d -r1.1 -r1.2 *** iem_adaptfilt.c 2 Aug 2006 14:02:28 -0000 1.1 --- iem_adaptfilt.c 11 Jan 2007 17:40:10 -0000 1.2 *************** *** 24,31 **** }
! void sigNLMS_setup(void); ! void sigNLMSCC_setup(void); ! void sign_CNLMS_setup(void); ! void sign_CLNLMS_setup(void);
/* ------------------------ setup routine ------------------------- */ --- 24,31 ---- }
! void NLMS_tilde_setup(void); ! void NLMSCC_tilde_setup(void); ! void n_CNLMS_tilde_setup(void); ! void n_CLNLMS_tilde_setup(void);
/* ------------------------ setup routine ------------------------- */ *************** *** 33,40 **** void iem_adaptfilt_setup(void) { ! sigNLMS_setup(); ! sigNLMSCC_setup(); ! sign_CNLMS_setup(); ! sign_CLNLMS_setup(); post("----------------------------------------------"); --- 33,40 ---- void iem_adaptfilt_setup(void) { ! NLMS_tilde_setup(); ! NLMSCC_tilde_setup(); ! n_CNLMS_tilde_setup(); ! n_CLNLMS_tilde_setup(); post("----------------------------------------------");