Update of /cvsroot/pure-data/externals/iem/iem_adaptfilt/src In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv13932/externals/iem/iem_adaptfilt/src
Added Files: .DS_Store iem_adaptfilt.c iemlib.h makefile.txt makefile_lin makefile_win makefile_win.txt sign_CLNLMS.c sign_CNLMS.c sigNLMS.c sigNLMSCC.c Log Message: no message
--- NEW FILE: makefile_lin --- current: all
.SUFFIXES: .pd_linux
INCLUDE = -I. -I/usr/local/src/pd-0.37-1/src
LDFLAGS = -export-dynamic -shared LIB = -ldl -lm -lpthread
#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
SYSTEM = $(shell uname -m)
# the sources
SRC = sigNLMS.c \ sigNLMSCC.c \ sign_CNLMS.c \ sign_CLNLMS.c \ iem_adaptfilt.c
TARGET = iem_adaptfilt.pd_linux
OBJ = $(SRC:.c=.o)
# # ------------------ targets ------------------------------------ #
clean: rm ../../lib/$(TARGET) rm *.o
all: $(OBJ) @echo :: $(OBJ) ld $(LDFLAGS) -o $(TARGET) *.o $(LIB) strip --strip-unneeded $(TARGET) rm *.o
$(OBJ) : %.o : %.c touch $*.c cc $(CFLAGS) $(INCLUDE) -c -o $*.o $*.c
--- NEW FILE: makefile_win.txt ---
all: iem_adaptfilt.dll
VIS_CPP_PATH = "C:\Programme\Microsoft Visual Studio\Vc98"
PD_INST_PATH = "C:\Programme\pd-0.37-1"
PD_WIN_INCLUDE_PATH = /I. /I$(PD_INST_PATH)\src /I$(VIS_CPP_PATH)\include
PD_WIN_C_FLAGS = /nologo /W3 /WX /DMSW /DNT /DPD /DWIN32 /DWINDOWS /Ox -DPA_LITTLE_ENDIAN
PD_WIN_L_FLAGS = /nologo
PD_WIN_LIB = /NODEFAULTLIB:libc /NODEFAULTLIB:oldnames /NODEFAULTLIB:kernel /NODEFAULTLIB:uuid \ $(VIS_CPP_PATH)\lib\libc.lib \ $(VIS_CPP_PATH)\lib\oldnames.lib \ $(VIS_CPP_PATH)\lib\kernel32.lib \ $(VIS_CPP_PATH)\lib\wsock32.lib \ $(VIS_CPP_PATH)\lib\winmm.lib \ $(PD_INST_PATH)\bin\pthreadVC.lib \ $(PD_INST_PATH)\bin\pd.lib
SRC = sigNLMS.c \ sigNLMSCC.c \ sign_CNLMS.c \ sign_CLNLMS.c \ iem_adaptfilt.c
OBJ = $(SRC:.c=.obj)
.c.obj: cl $(PD_WIN_C_FLAGS) $(PD_WIN_INCLUDE_PATH) /c $*.c
iem_adaptfilt.dll: $(OBJ) link $(PD_WIN_L_FLAGS) /dll /export:iem_adaptfilt_setup \ /out:iem_adaptfilt.dll $(OBJ) $(PD_WIN_LIB) copy iem_adaptfilt.dll ....\lib\iem_adaptfilt.dll copy iem_adaptfilt.dll ......\iem_adaptfilt.dll
clean: del *.obj
--- NEW FILE: sign_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 sign_CLNLMS_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_sign_CLNLMS_kern;
typedef struct sign_CLNLMS { t_object x_obj; t_sign_CLNLMS_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_sign_CLNLMS;
t_class *sign_CLNLMS_class;
static void sign_CLNLMS_tick(t_sign_CLNLMS *x) { outlet_bang(x->x_out_compressing_bang); }
static t_float *sign_CLNLMS_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 sign_CLNLMS_beta(t_sign_CLNLMS *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 sign_CLNLMS_gamma(t_sign_CLNLMS *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 sign_CLNLMS_kappa(t_sign_CLNLMS *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 sign_CLNLMS_leakage(t_sign_CLNLMS *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 sign_CLNLMS_update(t_sign_CLNLMS *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 *sign_CLNLMS_perform_zero(t_int *w) { t_sign_CLNLMS *x = (t_sign_CLNLMS *)(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 *sign_CLNLMS_perform(t_int *w) { t_sign_CLNLMS *x = (t_sign_CLNLMS *)(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 sign_CLNLMSperfzero;// 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);
sign_CLNLMSperfzero:
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 sign_CLNLMS_dsp(t_sign_CLNLMS *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 = sign_CLNLMS_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(sign_CLNLMS_perform_zero, 2, x, n); else dsp_add(sign_CLNLMS_perform, 2, x, n); }
/* setup/setdown things */
static void sign_CLNLMS_free(t_sign_CLNLMS *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_sign_CLNLMS_kern));
clock_free(x->x_clock); }
static void *sign_CLNLMS_new(t_symbol *s, t_int argc, t_atom *argv) { t_sign_CLNLMS *x = (t_sign_CLNLMS *)pd_new(sign_CLNLMS_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_sign_CLNLMS_kern *)getbytes(x->x_n_io*sizeof(t_sign_CLNLMS_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)sign_CLNLMS_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 sign_CLNLMS_setup(void) { sign_CLNLMS_class = class_new(gensym("n_CLNLMS~"), (t_newmethod)sign_CLNLMS_new, (t_method)sign_CLNLMS_free, sizeof(t_sign_CLNLMS), 0, A_GIMME, 0); CLASS_MAINSIGNALIN(sign_CLNLMS_class, t_sign_CLNLMS, x_msi); class_addmethod(sign_CLNLMS_class, (t_method)sign_CLNLMS_dsp, gensym("dsp"), 0); class_addmethod(sign_CLNLMS_class, (t_method)sign_CLNLMS_update, gensym("update"), A_FLOAT, 0); // method: downsampling factor of learning (multiple of 2^N) class_addmethod(sign_CLNLMS_class, (t_method)sign_CLNLMS_beta, gensym("beta"), A_FLOAT, 0); //method: normalized learning rate class_addmethod(sign_CLNLMS_class, (t_method)sign_CLNLMS_gamma, gensym("gamma"), A_FLOAT, 0); // method: dithering noise related to signal class_addmethod(sign_CLNLMS_class, (t_method)sign_CLNLMS_kappa, gensym("kappa"), A_FLOAT, 0); // method: threshold for compressing w_coeff class_addmethod(sign_CLNLMS_class, (t_method)sign_CLNLMS_leakage, gensym("leakage"), A_FLOAT, 0); // method: leakage factor [0 1] for w update class_sethelpsymbol(sign_CLNLMS_class, gensym("iemhelp2/n_CLNLMS~")); }
--- NEW FILE: sigNLMSCC.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 sigNLMSCC { 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_sigNLMSCC;
t_class *sigNLMSCC_class;
static void sigNLMSCC_tick(t_sigNLMSCC *x) { outlet_bang(x->x_out_clipping_bang); }
static t_float *sigNLMSCC_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 sigNLMSCC_beta(t_sigNLMSCC *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 sigNLMSCC_gamma(t_sigNLMSCC *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 sigNLMSCC_update(t_sigNLMSCC *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 *sigNLMSCC_perform_zero(t_int *w) { t_sigNLMSCC *x = (t_sigNLMSCC *)(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 *sigNLMSCC_perform(t_int *w) { t_sigNLMSCC *x = (t_sigNLMSCC *)(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 sigNLMSCCperfzero;// this is Musil/Miller style if(!wmin_filt_coeff) goto sigNLMSCCperfzero; if(!wmax_filt_coeff) goto sigNLMSCCperfzero;// 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);
sigNLMSCCperfzero:
while(n--) { *filt_out++ = 0.0f; *err_out++ = 0.0f; } return(w+3); }
static void sigNLMSCC_dsp(t_sigNLMSCC *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 = sigNLMSCC_check_array(x->x_w_array_sym_name, x->x_n_order); x->x_wmin_array_mem_beg = sigNLMSCC_check_array(x->x_wmin_array_sym_name, x->x_n_order); x->x_wmax_array_mem_beg = sigNLMSCC_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(sigNLMSCC_perform_zero, 2, x, n); else dsp_add(sigNLMSCC_perform, 2, x, n); }
/* setup/setdown things */
static void sigNLMSCC_free(t_sigNLMSCC *x) {
freebytes(x->x_in_hist, 2*x->x_n_order*sizeof(t_float));
clock_free(x->x_clock); }
static void *sigNLMSCC_new(t_symbol *s, t_int argc, t_atom *argv) { t_sigNLMSCC *x = (t_sigNLMSCC *)pd_new(sigNLMSCC_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)sigNLMSCC_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 sigNLMSCC_setup(void) { sigNLMSCC_class = class_new(gensym("NLMSCC~"), (t_newmethod)sigNLMSCC_new, (t_method)sigNLMSCC_free, sizeof(t_sigNLMSCC), 0, A_GIMME, 0); CLASS_MAINSIGNALIN(sigNLMSCC_class, t_sigNLMSCC, x_msi); class_addmethod(sigNLMSCC_class, (t_method)sigNLMSCC_dsp, gensym("dsp"), 0); class_addmethod(sigNLMSCC_class, (t_method)sigNLMSCC_update, gensym("update"), A_FLOAT, 0); // method: downsampling factor of learning (multiple of 2^N) class_addmethod(sigNLMSCC_class, (t_method)sigNLMSCC_beta, gensym("beta"), A_FLOAT, 0); //method: normalized learning rate class_addmethod(sigNLMSCC_class, (t_method)sigNLMSCC_gamma, gensym("gamma"), A_FLOAT, 0); // method: dithering noise related to signal class_sethelpsymbol(sigNLMSCC_class, gensym("iemhelp2/NLMSCC~")); }
--- NEW FILE: makefile_win ---
all: iem_adaptfilt.dll
VIS_CPP_PATH = "C:\Programme\Microsoft Visual Studio\Vc98"
PD_INST_PATH = "C:\Programme\pd-0.37-3"
PD_WIN_INCLUDE_PATH = /I. /I$(PD_INST_PATH)\src /I$(VIS_CPP_PATH)\include
PD_WIN_C_FLAGS = /nologo /W3 /WX /DMSW /DNT /DPD /DWIN32 /DWINDOWS /Ox -DPA_LITTLE_ENDIAN
PD_WIN_L_FLAGS = /nologo
PD_WIN_LIB = /NODEFAULTLIB:libc /NODEFAULTLIB:oldnames /NODEFAULTLIB:kernel /NODEFAULTLIB:uuid \ $(VIS_CPP_PATH)\lib\libc.lib \ $(VIS_CPP_PATH)\lib\oldnames.lib \ $(VIS_CPP_PATH)\lib\kernel32.lib \ $(VIS_CPP_PATH)\lib\wsock32.lib \ $(VIS_CPP_PATH)\lib\winmm.lib \ $(PD_INST_PATH)\bin\pthreadVC.lib \ $(PD_INST_PATH)\bin\pd.lib
SRC = sigNLMS.c \ sigNLMSCC.c \ sign_CNLMS.c \ sign_CLNLMS.c \ iem_adaptfilt.c
OBJ = $(SRC:.c=.obj)
.c.obj: cl $(PD_WIN_C_FLAGS) $(PD_WIN_INCLUDE_PATH) /c $*.c
iem_adaptfilt.dll: $(OBJ) link $(PD_WIN_L_FLAGS) /dll /export:iem_adaptfilt_setup \ /out:iem_adaptfilt.dll $(OBJ) $(PD_WIN_LIB)
clean: del *.obj
--- NEW FILE: sigNLMS.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 sigNLMS { 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_sigNLMS;
t_class *sigNLMS_class;
static t_float *sigNLMS_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 sigNLMS_beta(t_sigNLMS *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 sigNLMS_gamma(t_sigNLMS *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 sigNLMS_update(t_sigNLMS *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 *sigNLMS_perform_zero(t_int *w) { t_sigNLMS *x = (t_sigNLMS *)(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 *sigNLMS_perform(t_int *w) { t_sigNLMS *x = (t_sigNLMS *)(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 sigNLMSperfzero;// 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);
sigNLMSperfzero:
while(n--) { *filt_out++ = 0.0f; *err_out++ = 0.0f; } return(w+3); }
static void sigNLMS_dsp(t_sigNLMS *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 = sigNLMS_check_array(x->x_w_array_sym_name, x->x_n_order);
if(!x->x_w_array_mem_beg) dsp_add(sigNLMS_perform_zero, 2, x, n); else dsp_add(sigNLMS_perform, 2, x, n); }
/* setup/setdown things */
static void sigNLMS_free(t_sigNLMS *x) { freebytes(x->x_in_hist, 2*x->x_n_order*sizeof(t_float)); }
static void *sigNLMS_new(t_symbol *s, t_int argc, t_atom *argv) { t_sigNLMS *x = (t_sigNLMS *)pd_new(sigNLMS_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 sigNLMS_setup(void) { sigNLMS_class = class_new(gensym("NLMS~"), (t_newmethod)sigNLMS_new, (t_method)sigNLMS_free, sizeof(t_sigNLMS), 0, A_GIMME, 0); CLASS_MAINSIGNALIN(sigNLMS_class, t_sigNLMS, x_msi); class_addmethod(sigNLMS_class, (t_method)sigNLMS_dsp, gensym("dsp"), 0); class_addmethod(sigNLMS_class, (t_method)sigNLMS_update, gensym("update"), A_FLOAT, 0); // method: downsampling factor of learning (multiple of 2^N) class_addmethod(sigNLMS_class, (t_method)sigNLMS_beta, gensym("beta"), A_FLOAT, 0); //method: normalized learning rate class_addmethod(sigNLMS_class, (t_method)sigNLMS_gamma, gensym("gamma"), A_FLOAT, 0); // method: dithering noise related to signal class_sethelpsymbol(sigNLMS_class, gensym("iemhelp2/NLMS~")); }
--- NEW FILE: makefile.txt --- current: all
.SUFFIXES: .pd_linux
INCLUDE = -I. -I/usr/local/src/pd-0.37-1/src
LDFLAGS = -export-dynamic -shared LIB = -ldl -lm -lpthread
#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
SYSTEM = $(shell uname -m)
# the sources
SRC = sigNLMS.c \ sigNLMSCC.c \ sign_CNLMS.c \ iem_adaptfilt.c
TARGET = iem_adaptfilt.pd_linux
OBJ = $(SRC:.c=.o)
# # ------------------ targets ------------------------------------ #
clean: rm $(TARGET) rm *.o
all: $(OBJ) @echo :: $(OBJ) ld $(LDFLAGS) -o $(TARGET) *.o $(LIB) strip --strip-unneeded $(TARGET) rm *.o
$(OBJ) : %.o : %.c touch $*.c cc $(CFLAGS) $(INCLUDE) -c -o $*.o $*.c
--- NEW FILE: .DS_Store --- (This appears to be a binary file; contents omitted.)
--- NEW FILE: sign_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 sign_CNLMS_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_sign_CNLMS_kern;
typedef struct sign_CNLMS { t_object x_obj; t_sign_CNLMS_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_sign_CNLMS;
t_class *sign_CNLMS_class;
static void sign_CNLMS_tick(t_sign_CNLMS *x) { outlet_bang(x->x_out_compressing_bang); }
static t_float *sign_CNLMS_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 sign_CNLMS_beta(t_sign_CNLMS *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 sign_CNLMS_gamma(t_sign_CNLMS *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 sign_CNLMS_kappa(t_sign_CNLMS *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 sign_CNLMS_update(t_sign_CNLMS *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 *sign_CNLMS_perform_zero(t_int *w) { t_sign_CNLMS *x = (t_sign_CNLMS *)(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 *sign_CNLMS_perform(t_int *w) { t_sign_CNLMS *x = (t_sign_CNLMS *)(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 sign_CNLMSperfzero;// 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);
sign_CNLMSperfzero:
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 sign_CNLMS_dsp(t_sign_CNLMS *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 = sign_CNLMS_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(sign_CNLMS_perform_zero, 2, x, n); else dsp_add(sign_CNLMS_perform, 2, x, n); }
/* setup/setdown things */
static void sign_CNLMS_free(t_sign_CNLMS *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_sign_CNLMS_kern));
clock_free(x->x_clock); }
static void *sign_CNLMS_new(t_symbol *s, t_int argc, t_atom *argv) { t_sign_CNLMS *x = (t_sign_CNLMS *)pd_new(sign_CNLMS_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_sign_CNLMS_kern *)getbytes(x->x_n_io*sizeof(t_sign_CNLMS_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)sign_CNLMS_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 sign_CNLMS_setup(void) { sign_CNLMS_class = class_new(gensym("n_CNLMS~"), (t_newmethod)sign_CNLMS_new, (t_method)sign_CNLMS_free, sizeof(t_sign_CNLMS), 0, A_GIMME, 0); CLASS_MAINSIGNALIN(sign_CNLMS_class, t_sign_CNLMS, x_msi); class_addmethod(sign_CNLMS_class, (t_method)sign_CNLMS_dsp, gensym("dsp"), 0); class_addmethod(sign_CNLMS_class, (t_method)sign_CNLMS_update, gensym("update"), A_FLOAT, 0); // method: downsampling factor of learning (multiple of 2^N) class_addmethod(sign_CNLMS_class, (t_method)sign_CNLMS_beta, gensym("beta"), A_FLOAT, 0); //method: normalized learning rate class_addmethod(sign_CNLMS_class, (t_method)sign_CNLMS_gamma, gensym("gamma"), A_FLOAT, 0); // method: dithering noise related to signal class_addmethod(sign_CNLMS_class, (t_method)sign_CNLMS_kappa, gensym("kappa"), A_FLOAT, 0); // method: threshold for compressing w_coeff class_sethelpsymbol(sign_CNLMS_class, gensym("iemhelp2/n_CNLMS~")); }
--- NEW FILE: iem_adaptfilt.c --- /* For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution.
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"
static t_class *iem_adaptfilt_class;
static void *iem_adaptfilt_new(void) { t_object *x = (t_object *)pd_new(iem_adaptfilt_class);
return (x); }
void sigNLMS_setup(void); void sigNLMSCC_setup(void); void sign_CNLMS_setup(void); void sign_CLNLMS_setup(void);
/* ------------------------ setup routine ------------------------- */
void iem_adaptfilt_setup(void) { sigNLMS_setup(); sigNLMSCC_setup(); sign_CNLMS_setup(); sign_CLNLMS_setup(); post("----------------------------------------------"); post("iem_adaptfilt (R-1.02) library loaded!"); post("(c) Markus Noisternig, Thomas Musil"); post(" {noisternig, musil}_AT_iem.at"); post(" IEM Graz, Austria"); post("----------------------------------------------"); }
--- NEW FILE: iemlib.h --- /* For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution.
iemlib2 written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2004 */
#ifndef __IEMLIB_H__ #define __IEMLIB_H__
#define IS_A_POINTER(atom,index) ((atom+index)->a_type == A_POINTER) #define IS_A_FLOAT(atom,index) ((atom+index)->a_type == A_FLOAT) #define IS_A_SYMBOL(atom,index) ((atom+index)->a_type == A_SYMBOL) #define IS_A_DOLLAR(atom,index) ((atom+index)->a_type == A_DOLLAR) #define IS_A_DOLLSYM(atom,index) ((atom+index)->a_type == A_DOLLSYM) #define IS_A_SEMI(atom,index) ((atom+index)->a_type == A_SEMI) #define IS_A_COMMA(atom,index) ((atom+index)->a_type == A_COMMA)
#ifdef NT int sys_noloadbang; //t_symbol *iemgui_key_sym=0; #include <io.h> #else extern int sys_noloadbang; //extern t_symbol *iemgui_key_sym; #include <unistd.h> #endif
#define DEFDELVS 64 #define XTRASAMPS 4 #define SAMPBLK 4
#define UNITBIT32 1572864. /* 3*2^19; bit 32 has place value 1 */
/* machine-dependent definitions. These ifdefs really should have been by CPU type and not by operating system! */ #ifdef IRIX /* big-endian. Most significant byte is at low address in memory */ #define HIOFFSET 0 /* word offset to find MSB */ #define LOWOFFSET 1 /* word offset to find LSB */ #define int32 long /* a data type that has 32 bits */ #else #ifdef MSW /* little-endian; most significant byte is at highest address */ #define HIOFFSET 1 #define LOWOFFSET 0 #define int32 long #else #ifdef __FreeBSD__ #include <machine/endian.h> #if BYTE_ORDER == LITTLE_ENDIAN #define HIOFFSET 1 #define LOWOFFSET 0 #else #define HIOFFSET 0 /* word offset to find MSB */ #define LOWOFFSET 1 /* word offset to find LSB */ #endif /* BYTE_ORDER */ #include <sys/types.h> #define int32 int32_t #endif #ifdef __linux__
#include <endian.h>
#if !defined(__BYTE_ORDER) || !defined(__LITTLE_ENDIAN) #error No byte order defined #endif
#if __BYTE_ORDER == __LITTLE_ENDIAN #define HIOFFSET 1 #define LOWOFFSET 0 #else #define HIOFFSET 0 /* word offset to find MSB */ #define LOWOFFSET 1 /* word offset to find LSB */ #endif /* __BYTE_ORDER */
#include <sys/types.h> #define int32 int32_t
#else #ifdef __APPLE__ #define HIOFFSET 0 /* word offset to find MSB */ #define LOWOFFSET 1 /* word offset to find LSB */ #define int32 int /* a data type that has 32 bits */
#endif /* __APPLE__ */ #endif /* __linux__ */ #endif /* MSW */ #endif /* SGI */
union tabfudge { double tf_d; int32 tf_i[2]; };
#define IEM_DENORMAL(f) ((((*(unsigned int*)&(f))&0x60000000)==0) || \ (((*(unsigned int*)&(f))&0x60000000)==0x60000000)) /* more stringent test: anything not between 1e-19 and 1e19 in absolute val */
#endif