Writing MEX code for LMS filter: some difficulties with my implementation

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
end

My 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

[Kaustubha Govind]

What error do you see?

[dm]

When testing with a random input signal, i.e. x = rand(N,1)+1i*rand(N,1), and letting the desired signal be i.e. d = filter([0.5,0.3,0.1],1,x), the error that I see is that the output deviate quite a lot from the output when calling the MATLAB implementation.

When testing with a set of simulated input/output signals (16QAM, ~20.000 samples) of a RF power amplifier (data acquired from an envelope simulation of a RFPA in Agilent ADS), the MATLAB function yields quite good results (or at least a decreasing error for increasing number of samples processed), whereas the MEX function diverge after quite few samples (between 100 to 1000).

Another funny thing: running the following code yields a different output from the MEX function for each run, even when no parameters are changed. The MATLAB function yields the same result each time (thus being robust). So it might be some numerical inaccuracies here as well which trick me, but I can't tell.

clear all; clc;

load ads_pavg20; % load set of input/output data

g = meangain(x,y,50); % find mean gain of RFPA 16QAM

d = g*x; % desired signal

M = 16;

mu = 0.00008; % LMS step size

[w1 y1 e1 W1] = adaptfilt_lms(x,d,mu,M); % MATLAB function

[w2 y2 e2 W2] = adaptfilt_lms_mex(x,d,mu,M); % MEX function

Answer 1

[dm]

I found the error;

wi[m] = wi[m]+2*mu*(er[n]*xti[n]-ei[n]*xtr[m]); /*

should be

wi[m] = wi[m]+2*mu*(er[n]*xti[m]-ei[n]*xtr[m]); /*

Managed to use index "n" where I should have used "m" in one array.

Now just some small tweaking left I assume.

[Fahd Omrani]

for j=1:26

alphabets(j)=char(64+j);

end