cov_nd.c

toolbox of BVQX, This is the access between BV and matlab. It will help you to analysis data from BV

/*

calculation of covariance and correlation for pairs of vectors in N-D matrices

% Version:  v0.5c
% Build:    6120415
% Date:     Dec-04 2006, 3:15 PM CET
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

*/

#include "mex.h"
#include "matrix.h"
#include "math.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	int docr, c, ic, nd, nl, dr[64], rd, rs, ts;
	const int *d1, *d2;
	double dnl, dnx, *p1, *p2, *t1, *t2;
	long double m1, m2, x1, x2, cv;
	double mp;

	if (nrhs != 2 || nlhs < 1 || nlhs > 2)
		mexErrMsgTxt("Bad number of input/output arguments.");

	if (nlhs > 1)
		docr = 1;
	    else	
		docr = 0;

	nd = mxGetNumberOfDimensions(prhs[0]);
        if (nd != mxGetNumberOfDimensions(prhs[1]) || nd > 64)
		mexErrMsgTxt("Matrices must have the same size.");

	d1 = mxGetDimensions(prhs[0]);
	d2 = mxGetDimensions(prhs[1]);
	for(c=0;c<nd;++c)
		if (d1[c] != d2[c] || d1[c] == 0)
			mexErrMsgTxt("Matrices must have the same size and not be empty.");

	for(c=0; c<2; ++c)
		if (!mxIsNumeric(prhs[c]) ||  mxIsComplex(prhs[c]) ||
		     mxIsSparse(prhs[c])  || !mxIsDouble(prhs[c]))
			mexErrMsgTxt("Matrices must be numeric, real, full and double.");

	rd = nd-1;
	if (d1[rd] == 1)
		mexErrMsgTxt("The last matrix dimensions must not be singleton.");

	p1 = (double *) mxGetPr(prhs[0]);
	p2 = (double *) mxGetPr(prhs[1]);

	rs = 1;
	for(c=0;c<rd;++c) {
		dr[c] = d1[c];
		rs *= dr[c];
	}
	nl = d1[rd];
	ts = rs * nl;
	dnl = (double) nl;
	dnx = dnl - 1.0;

	plhs[0] = mxCreateNumericArray(rd, dr, mxDOUBLE_CLASS, mxREAL);
	t1 = (double *) mxGetPr(plhs[0]);

	if (!docr) {

		for(c=0;c<rs;++c) {
			m1 = 0.0;
			m2 = 0.0;
			mp = 0.0;
			for(ic=c;ic<ts;ic+=rs) {
				m1 += p1[ic];
				m2 += p2[ic];
				mp += p1[ic]*p2[ic];
			}
			m1 /= dnl;
			m2 /= dnl;
			mp /= dnl;
			t1[c] = (double) ((mp - m1*m2) * (dnl/dnx));
		}

	} else {

		plhs[1] = mxCreateNumericArray(rd, dr, mxDOUBLE_CLASS, mxREAL);
		t2 = (double *) mxGetPr(plhs[1]);
		for(c=0;c<rs;++c) {
			m1 = 0.0;
			m2 = 0.0;
			x1 = 0.0;
			x2 = 0.0;
			mp = 0.0;
			for(ic=c;ic<ts;ic+=rs) {
				m1 += p1[ic];
				m2 += p2[ic];
				x1 += p1[ic]*p1[ic];
				x2 += p2[ic]*p2[ic];
				mp += p1[ic]*p2[ic];
			}
			m1 /= dnl;
			m2 /= dnl;
			x1 /= dnl;
			x2 /= dnl;
			mp /= dnl;
			cv = (mp - m1*m2) * (dnl/dnx);
			t1[c] = (double) cv;
			t2[c] = (double) (cv / sqrt((x1-m1*m1)*(x2-m2*m2)) / (dnl/dnx));
		}
	}
}