Hi all!
I am trying to write a Matlab mex-file, and I find the c-coding part a bit challenging. I can get the code below to compile, but it doesn't do what is expected:
Code:
#include <math.h>
#include "mex.h"
/*Subroutine to calculate the Fuzzy Sample Entropy*/
void FSE(double *cop, int M, double r, double c, mwSize N, double *FuzzSampEnt) {
mwIndex m, i, ii, j, jj, k, l;
double expo, d, *seg_i, *seg_j, count, sumCim, *Cim, CIM[2];
expo=(log(log(pow(2, c)))/log(r));
for (m=M;m<M+2;m++) {
seg_i=mxCalloc(m, sizeof(double));
seg_j=mxCalloc(m, sizeof(double));
Cim=mxCalloc(N-m, sizeof(double));
for (i=0; i<N-m; i++) {
count=0;
for (ii=0; ii<m; ii++){
seg_i[ii]=cop[i+ii];//TROUBLE
}
for (j=0; j<N-m; j++) {
for (jj=0; jj<m; jj++) {
seg_j[jj]=cop[j+jj];//TROUBLE
}
d=0;
/*Finding the maximum absolut value of seg_i-seg_j*/
for (k=0; k<m; k++) {
if (d < fabs(seg_i[k]-seg_j[k])) {
d=fabs(seg_i[k]-seg_j[k]);
}
}
/*if d=0, self match*/
if (d>0) {
count+=exp(-(pow(d, expo/c)));
}
}
Cim[i]=count/(N-m-1);
}
sumCim=0;
for(l=0;l<N-m;l++){
sumCim+=Cim[l];
}
if (m==M){
CIM[0]=sumCim/(N-m);
}
else if (m==M+1){
CIM[1]=sumCim/(N-m);
}
}
*FuzzSampEnt=seg_j[k];//-log(CIM[1]/ CIM[0]);
mxFree(seg_i);
mxFree(seg_j);
mxFree(Cim);
}
/*the gateway function*/
void mexFunction(int nlhs, mxArray *plhs[ ],
int nrhs, const mxArray *prhs[ ]) {
int M;
double *cop, *FuzzSampEnt;
double c, r;
mwSize N;
if(nrhs!=4) {
mexErrMsgTxt("Four inputs required (cop,M,r,c).");
}
/*Assign pointers*/
cop=mxGetPr(prhs[0]);
N=mxGetM(prhs[0]);
M=(int)(mxGetScalar(prhs[1]));
r=mxGetScalar(prhs[2]);
c=mxGetScalar(prhs[3]);
plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
FuzzSampEnt = mxGetPr(plhs[0]);
/*Call the subroutine */
FSE(cop, M, r, c, N, FuzzSampEnt);
return;
}
The trouble seems to be here: seg_i[ii]=cop[i+ii]; cop[i+ii] seems to be what it should be, but I can't seem to get it's value to seg_i[ii]. Also I'm not sure if it is necessary to use pointers here at all. I know that seg_i and seg_j are either M or M+1 in length and Cim is either N-M or N-M-1 in length.
The idea would be to make the code as fast as possible, as fast as possible.. So how would you deal with the seg_i, seg_j and Cim? Can you see what I'm doing wrong in here?
Thank You!
-Aino