I'm trying to implement a MEX function for a LMS filter (which I've already implemented in MATLAB (*.m)) in order to improve my C/C++ and maybe learn some MEX. 🙂
Unfortunately (for my part) I'm still quite new to C/C++, and experience some difficulties with the implementation (which I thought would be piece of cake considering the easy implementation in *.m).
My *.m code for the LMS filter is as following
function [w,y,e,W] = adaptfilt_lms(x,dn,mu,M)N = length(x);y = zeros(N,1); % initialize filter output vector
w = zeros(M,1); % initialize filter coefficient vector
e = zeros(N,1); % initialize error vector
W = zeros(M,N); % filter coefficient matrix for coeff. history
for n = 1:N if n <= M % assume zero-samples if no data available
k = n:-1:1; xT = [x(k); zeros(M-numel(k),1)]; else xT = x(n-M+1:1:n); % M samples of x in reverse order
end y(n) = w'*xT; % filter output
e(n) = dn(n)-y(n); % error
w = w+2*mu*e(n)'*xT; % update filter coefficients
W(:,n) = w; % store current filter coefficients in matrix
endMy MEX implementation is as following:
#include <math.h>#include "mex.h"void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ /* Macros for input data */ #define xData prhs[0] /* Input signal to filter */ #define dData prhs[1] /* Desired output signal */ #define muData prhs[2] /* Step-size */ #define MData prhs[3] /* Filter order */ /* Macros for output data */ #define wData plhs[0] /* Final filter coefficients */ #define yData plhs[1] /* Filter output signal vector */ #define eData plhs[2] /* Error vector */ #define WData plhs[3] /* Filter coefficients matrix */ mxArray *xTemp; double *xr, *xi, *dr, *di, *yr, *yi; double *wr, *wi, *er, *ei, *Wr, *Wi; double *xtr, *xti, mu; int M, Mx, Md, Nx, Nd, m, n, k, nT; /* Get pointers to input data */ xData = prhs[0]; /* Get filter input signal */ xr = mxGetPr(xData); /* Real part data */ xi = mxGetPi(xData); /* Imaginary part data */ dData = prhs[1]; /* Desired signal */ dr = mxGetPr(dData); /* Real part data */ di = mxGetPi(dData); /* Imaginary part data */ muData = prhs[2]; /* Get step-size constant */ mu = mxGetScalar(muData); MData = prhs[3]; /* Get filter order */ M = (int)mxGetScalar(MData); /* Get dimensions of input data */ Mx = mxGetM(xData); /* Number of rows in x */ Nx = mxGetN(xData); /* Number of columns in x */ Md = mxGetM(dData); /* Number of rows in d */ Nd = mxGetN(dData); /* Number of columsn in d */ /* Temporary vector - size M-by-1 */ xTemp = mxCreateNumericMatrix(M, 1, mxDOUBLE_CLASS, mxCOMPLEX); xtr = mxGetPr(xTemp); /* Real part data */ xti = mxGetPi(xTemp); /* Imaginary part data */ /* Create output vector - size Mx-by-1 */ yData = mxCreateNumericMatrix(Mx, 1, mxDOUBLE_CLASS, mxCOMPLEX); yr = mxGetPr(yData); /* Pointer to the real part data of y */ yi = mxGetPi(yData); /* Pointer to the imaginary part data of y */ /* Create error vector - size Mx-by-1 */ eData = mxCreateNumericMatrix(Mx, 1, mxDOUBLE_CLASS, mxCOMPLEX); er = mxGetPr(eData); /* Pointer to the real part data of e */ ei = mxGetPi(eData); /* Pointer to the imaginary part data of e */ /* Create coeff. vector - size M-by-1 */ wData = mxCreateNumericMatrix(M, 1, mxDOUBLE_CLASS, mxCOMPLEX); wr = mxGetPr(wData); /* Pointer to the real part data of w */ wi = mxGetPi(wData); /* Pointer to the imaginary part data of w */ /* Create coeff. matrix - size M-by-Mx */ WData = mxCreateNumericMatrix(M, Mx, mxDOUBLE_CLASS, mxCOMPLEX); Wr = mxGetPr(WData); /* Pointer to the real part data of W */ Wi = mxGetPi(WData); /* Pointer to the imaginary data of W */ for(n = 0; n < Mx; n++) { /* Compute xT_m = [x[n] x[n-1] .... x[n-(M-1)]] */ for(m = 0; m < M; m++) { /* Assume zero's outside available data and zero-fill*/ if(n < (M-1)) { nT = n; for(k = 0; k < M; k++) { xtr[k] = xr[nT]; xti[k] = xi[nT]; if(nT == 0) break; else nT--; } } else /* Data available for all lags */ { xtr[m] = xr[n-m]; xti[m] = xi[n-m]; } } /* Compute y[n] = conj(w)*xT */ for(m = 0; m < M; m++) { yr[n] = yr[n]+(wr[m]*xtr[m]+wi[m]*xti[m]); /* Real */ yi[n] = yi[n]+(wr[m]*xti[m]-wi[m]*xtr[m]); /* Imag */ } /* Compute filter output error */ er[n] = dr[n]-yr[n]; /* Real part of e[n] */ ei[n] = di[n]-yi[n]; /* Imaginary part of e[n] */ /* Update filter coefficients w and store copy in W */ for(m = 0; m < M; m++) { /* Update filter coefficients */ /* w = w + 2*mu*conj(e)*xT */ wr[m] = wr[m]+2*mu*(er[n]*xtr[m]+ei[n]*xti[m]); /* Real */ wi[m] = wi[m]+2*mu*(er[n]*xti[n]-ei[n]*xtr[m]); /* Imag */ /* Store current version of coefficients in W */ Wr[m+(int)M*n] = wr[m]; /* Real part */ Wi[m+(int)M*n] = wi[m]; /* Imaginary part */ } } return; }I think the error originate from the part where I update the filter coefficients, but to be honest I'm not quite sure. Might be that the implementation can be simplified greatly as well (all tips and hints are greatly appreciated).
best regards,
dm
Best Answer