Matlab MEX Files

1
Matlab MEX Files

• Prof. David V. Anderson
• Electrical and Computer Engineering
• Georgia Institute of Technology

2
What are Mex Files?

• MEX Matlab EXecutable
• Dynamically Linked Libraries
• Used like a .m function
• Written in C (or Fortran)

3
When/Why Do I Use MEX Fcns?

• To speed up slow functions
• (Use profile command to find trouble spots)
• To gradually convert from Matlab to C
• Convert single functions at a time and test
  incrementally
• To use existing C code with Matlab

4
Example LMS filtering for system identification
e(n)
x(n)
d(n)
S

-
W
wi(n) wi(n-1) 2µx(n-i)e(n)
5
LMS Matlab Code

• create random input sequence
• x randn(1,1000)
• "unknown" filter
• h 1 -1 .5 -.5 .2 .1 .1
• filtered ("desired") signal
• d filter(h,1,x)
• initial filter weight guess
• N 10
• w zeros(1,N)
• adaptation factor
• mu 0.008
• now implement the adaptive filter
• for n Nlength(x),
• produce filtered output sample

6
Profile

• gtgt profile on
• gtgt taco
• gtgt profile report
• 100 of the total time in this function was spent
  on the following lines
•  
• 19   produce filtered output sample
• 0.030 37 20  dhat(n)  w  x(n-1n-(N-1))'  
  
• 21   update the filter coefficients
   
• 22  e(n)  d(n) - dhat(n)
• 0.050 63 23  w  w2mux(n-1n-(N-1))e(n)
   
• 24  end

7
Steps to Convert

• Rewrite Matlab code with sub-functions to
  optimize
• function dhat,e,w lms(x,d,mu,w)
• if (nargin lt 4)
• initial filter weight guess
• N 10
• w zeros(1,N)
• else
• N length(w)
• end
• now implement the adaptive filter
• for n Nlength(x),
• produce filtered output sample
• dhat(n) w x(n-1n-(N-1))'
• update the filter coefficients
• e(n) d(n) - dhat(n)
• w w2mux(n-1n-(N-1))e(n)
• end

8
Components of MEX-files

• Note, much of the remaining material for this
  talk is taken from the MathWorks help document
  1615, page 1
• Components of MEX-files
• There are a number of issues that need to be
  understood in order to write your own C
  MEX-files.
• The source code for a MEX-file consists of
  Computational routine - which contains the actual
  code that you want implemented in the MEX-file.
  It is called by the gateway routine (for short
  MEX-files, the computational routine code can be
  included in the gateway).
• Gateway routine - which interfaces the
  computational routine with MATLAB It calls the
  computational routine.

9
Components of MEX-files cont.

• All C MEX-file source code must contain the
  statement include "mex.h"
  so that the entry point and the interface
  routines are declared properly.
• The entry point to the gateway routine must be
  named mexFunction. This function acts similarly
  to the function main() in a standalone C program,
  in that it is where the function MEX-file starts
  executing when MATLAB calls it. void
  mexFunction(int nlhs, mxArray plhs,
  int nrhs, const mxArray prhs) NOTE Notice
  the const keyword here, the data associated with
  prhs can be viewed, however, we promise not to
  change the data that prhs points to.

10
Simple MEX-file

• include "mex.h"
• void mexFunction(int nlhs,
• mxArray plhs,
• int nrhs,
• const mxArray prhs)
• mexPrintf("Hello world\n")

11
Matlab Arrays

• The mxArray declaration corresponds to the
  internal data structure that MATLAB uses to
  represent arrays. The mxArray structure contains
  among other things
• The MATLAB variables name
• Its dimensions
• Its type
• Whether the variable is real or complex.
• If the variable contains complex numbers as
  elements, the MATLAB array includes vectors
  containing the real and imaginary parts.
• The data in the matrix is stored in these arrays
  columnwise

12
Array Example

• gtgt A 1.0000 1.0000i 2.0000
  3.0000
• 4.0000 5.0000 1.0000i
  6.0000
• would be stored in the matrix (or C structure) as
  mxArray.name "A" mxArray.m 2 mxArray.n 3
  mxArray.pr0 1.0 mxArray.pi0
  1.0mxArray.pr1 4.0 mxArray.pi1
  0.0mxArray.pr2 2.0 mxArray.pi2
  0.0mxArray.pr3 5.0 mxArray.pi3
  1.0mxArray.pr4 3.0 mxArray.pi4
  0.0mxArray.pr5 6.0 mxArray.pi5 0.0
• Note that if a matrix is purely real (all
  imaginary elements are zero), then the pi pointer
  equals NULL.

13
MEX-file Inputs and Outputs

• void mexFunction(int nlhs, mxArray plhs,
  int nrhs, const mxArray prhs)

14
MEX-file Shell

• Now write a function that serves as the interface
  to Matlab.
• Does not need to be changed much when you write
  new MEX-files
• Calls your processing routine

15
Memory Usage Notes

• It is important that you NOT use malloc(),
  calloc(), etc.
• Use mxCalloc() or mxMalloc
• Also, you must free up variables when you are
  done
• Use mxDestroyArray
• Do not use this on the input arguments
• It is possible to allocate static memory
• You should also register an exit function if you
  do

16
Special MEX Features

• You can access data in the context of the calling
  function
• You can call Matlab functions from MEX functions
• C-MEX functions can access any type of data that
  Matlab uses

17
lms_mex.c

• include "mex.h"
• void dolms(double x, double d, double dhat,
  double e, double w, double mu, int N, int
  xlen)
• / Macros /
• define max(a,b) ((a) gt (b) ? (a) (b))
• define DEFAULT_FILTER_SIZE 10
• static mxArray persistent_wArray NULL
• static int initialized 0
• void cleanup(void)
• if (initialized)
• mexPrintf("MEX-file is terminating,
  destroying array\n")
• mxDestroyArray(persistent_wArray)

18
lms_mex.c cont.

• /
• Dereference input arguments
• /
• / Get x, xlen /
• x mxGetPr(prhs0)
• m mxGetM(prhs0)
• n mxGetN(prhs0)
• if (m ! 1 n ! 1)
• mexErrMsgTxt("x must be a column or row
  vector.\n\n Usage\n dhat,e,w
  lms(x,d,mu,w)")
• xlen max(m,n)
• / Get d, dlen /
• d mxGetPr(prhs1)
• m mxGetM(prhs1)
• n mxGetN(prhs1)
• if (m ! 1 n ! 1)
• mexErrMsgTxt("d must be a column or row
  vector.\n\n Usage\n dhat,e,w
  lms(x,d,mu,w)")
• dlen max(m,n)

19
lms_mex.c cont.

• if (nrhs 4)
• m mxGetM(prhs3)
• n mxGetN(prhs3)
• if (m ! 1 n ! 1)
• mexErrMsgTxt("w must be a column or row
  vector.\n\n Usage\n dhat,e,w
  lms(x,d,mu,w)")
• N max(m,n)
• / if we already have a persistent array, we
  need to deallocate it unless
• it just happens to be the same size as the
  input weights /
• if (initialized)
• if (N ! mxGetM(persistent_wArray))
• mxDestroyArray(persistent_wArray)
• initialized 0
• mexPrintf("Resizing persistent array\n")
• /
• Initialize Persistant Array

20
lms_mex.c cont

• w mxGetPr(persistent_wArray)
• / If initial filter taps were passed then copy
  them to the persistent array /
• if (nrhs 4)
• w0 mxGetPr(prhs3)
• for(j0 jltN j)
• wn w0n
• / Create matrices for the return argument. /
• dhatArray mxCreateDoubleMatrix(dlen, 1,
  mxREAL)
• dhat mxGetPr(dhatArray)
• eArray mxCreateDoubleMatrix(dlen, 1, mxREAL)
• e mxGetPr(eArray)
• / mexPrintf("Getting ready to call dolms\n")
  /
• dolms(x, d, dhat, e, w, 2.0mu, N, xlen)

21
lms_mex.c cont

• if (nlhs gt 1)
• plhs1 eArray
• else
• mxDestroyArray(eArray)
• if (nlhs gt 2)
• plhs2 mxCreateDoubleMatrix(N, 1, mxREAL)
• wout mxGetPr(plhs2)
• for(j0jltNj)
• woutj wj

22
dolms.c

• void dolms(double x, double d, double dhat,
  double e, double w, double mu, int N, int xlen)
• //
• // lms(x,d,dhat,e,w,xlen,N,mu)
• // double x / pointer to
  input data buffer /
• // double d / pointer to
  desired signal buffer /
• // double w / weight vector
  /
• // double dhat / estimate of d
  produced by filter /
• // double e / error /
• // int xlen / length of
  input data buffer /
• // int N / filter length
  /
• // double mu / mu 2 x's
  the conventional mu /
• / Note the input data is in a buffer that is
  actually LN long /
• / so there will be enough data to fill the
  filter taps as well. /
• register double s / summer used in
  filter /
• register int j,i / loop counters /

23
dolms.c cont.

• /
• Matlab equivalent code
• for n Nlength(x),
• produce filtered output sample
• dhat(n) w x(n-1n-(N-1))'
• update the filter coefficients
• e(n) d(n) - dhat(n)
• w w2mux(n-1n-(N-1))e(n)
• end
• /
• for(iN ilt xlen i)
• for(wn w, x1xi, j0,s0.0 jltN j)
• s (wn) (x1--)
• ei di - s
• dhati s
• s mu ei
• for(j0 jltN j)

24
The Finished Product!!!

• Notice that the execution time went from 8.62
  seconds to 0.03 secondsan improvement of a
  factor of 287

25
(No Transcript)