Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/KPMstats/parzenC.c @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 /* C mex version of parzen.m
wolffd@0 2 [B,B2] = parzen(feat, mu, Sigma, Nproto);
wolffd@0 3 */
wolffd@0 4 #include "mex.h"
wolffd@0 5 #include <stdio.h>
wolffd@0 6 #include <math.h>
wolffd@0 7
wolffd@0 8 #define PI 3.141592654
wolffd@0 9
wolffd@0 10 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0 11 int D, M, Q, T, d, m, q, t;
wolffd@0 12 double *data, *mu, *SigmaPtr, *N, Sigma;
wolffd@0 13 double *B, *dist, *B2, tmp;
wolffd@0 14 const int* dim_mu;
wolffd@0 15 double const1, const2, sum_m, sum_d, diff;
wolffd@0 16 int Dt, DMq, Dm, MQt, Mq;
wolffd@0 17 int dims_B2[3];
wolffd@0 18
wolffd@0 19 int ndim_mu, i, save_B2;
wolffd@0 20
wolffd@0 21 data = mxGetPr(prhs[0]);
wolffd@0 22 mu = mxGetPr(prhs[1]);
wolffd@0 23 SigmaPtr = mxGetPr(prhs[2]);
wolffd@0 24 Sigma = *SigmaPtr;
wolffd@0 25 N = mxGetPr(prhs[3]);
wolffd@0 26
wolffd@0 27 D = mxGetM(prhs[0]);
wolffd@0 28 T = mxGetN(prhs[0]);
wolffd@0 29
wolffd@0 30 ndim_mu = mxGetNumberOfDimensions(prhs[1]);
wolffd@0 31 dim_mu = mxGetDimensions(prhs[1]);
wolffd@0 32 D = dim_mu[0];
wolffd@0 33 M = dim_mu[1];
wolffd@0 34 /* printf("parzenC: nlhs=%d, D=%d, M=%d, T=%d\n", nlhs, D, M, T); */
wolffd@0 35
wolffd@0 36 /* If mu is mu(d,m,o,p), then [d M Q] = size(mu) in matlab sets Q=o*p,
wolffd@0 37 i.e.. the size of all conditioning variabeles */
wolffd@0 38 Q = 1;
wolffd@0 39 for (i = 2; i < ndim_mu; i++) {
wolffd@0 40 /* printf("dim_mu[%d]=%d\n", i, dim_mu[i]); */
wolffd@0 41 Q = Q*dim_mu[i];
wolffd@0 42 }
wolffd@0 43
wolffd@0 44 /* M = max(N) */
wolffd@0 45 M = -1000000;
wolffd@0 46 for (i=0; i < Q; i++) {
wolffd@0 47 /* printf("N[%d]=%d\n", i, (int) N[i]); */
wolffd@0 48 if (N[i] > M) {
wolffd@0 49 M = (int) N[i];
wolffd@0 50 }
wolffd@0 51 }
wolffd@0 52
wolffd@0 53 /* printf("parzenC: nlhs=%d, D=%d, Q=%d, M=%d, T=%d\n", nlhs, D, Q, M, T); */
wolffd@0 54
wolffd@0 55 plhs[0] = mxCreateDoubleMatrix(Q,T, mxREAL);
wolffd@0 56 B = mxGetPr(plhs[0]);
wolffd@0 57
wolffd@0 58 if (nlhs >= 2)
wolffd@0 59 save_B2 = 1;
wolffd@0 60 else
wolffd@0 61 save_B2 = 0;
wolffd@0 62
wolffd@0 63 if (save_B2) {
wolffd@0 64 /* printf("parzenC saving B2\n"); */
wolffd@0 65 /*plhs[1] = mxCreateDoubleMatrix(M*Q*T,1, mxREAL);*/
wolffd@0 66 dims_B2[0] = M;
wolffd@0 67 dims_B2[1] = Q;
wolffd@0 68 dims_B2[2] = T;
wolffd@0 69 plhs[1] = mxCreateNumericArray(3, dims_B2, mxDOUBLE_CLASS, mxREAL);
wolffd@0 70 B2 = mxGetPr(plhs[1]);
wolffd@0 71 } else {
wolffd@0 72 /* printf("parzenC not saving B2\n"); */
wolffd@0 73 }
wolffd@0 74 /*
wolffd@0 75 plhs[2] = mxCreateDoubleMatrix(M*Q*T,1, mxREAL);
wolffd@0 76 dist = mxGetPr(plhs[2]);
wolffd@0 77 */
wolffd@0 78 const1 = pow(2*PI*Sigma, -D/2.0);
wolffd@0 79 const2 = -(1/(2*Sigma));
wolffd@0 80
wolffd@0 81 for (t=0; t < T; t++) {
wolffd@0 82 /* printf("t=%d!\n",t); */
wolffd@0 83 Dt = D*t;
wolffd@0 84 MQt = M*Q*t;
wolffd@0 85 for (q=0; q < Q; q++) {
wolffd@0 86 sum_m = 0;
wolffd@0 87 DMq = D*M*q;
wolffd@0 88 Mq = M*q;
wolffd@0 89
wolffd@0 90 for (m=0; m < (int)N[q]; m++) {
wolffd@0 91 sum_d = 0;
wolffd@0 92 Dm = D*m;
wolffd@0 93 for (d=0; d < D; d++) {
wolffd@0 94 /* diff = data(d,t) - mu(d,m,q) */
wolffd@0 95 /*diff = data[d + D*t] - mu[d + D*m + D*M*q]; */
wolffd@0 96 diff = data[d + Dt] - mu[d + Dm + DMq];
wolffd@0 97 sum_d = sum_d + diff*diff;
wolffd@0 98 }
wolffd@0 99 /* dist[m,q,t] = dist[m + M*q + M*Q*t] = dist[m + Mq + MQt] = sum_d */
wolffd@0 100 tmp = const1 * exp(const2*sum_d);
wolffd@0 101 sum_m = sum_m + tmp;
wolffd@0 102 if (save_B2)
wolffd@0 103 B2[m + Mq + MQt] = tmp;
wolffd@0 104 }
wolffd@0 105
wolffd@0 106 if (N[q]>0) {
wolffd@0 107 B[q + Q*t] = (1.0/N[q]) * sum_m;
wolffd@0 108 } else {
wolffd@0 109 B[q + Q*t] = 0.0;
wolffd@0 110 }
wolffd@0 111 }
wolffd@0 112 }
wolffd@0 113 }
wolffd@0 114
wolffd@0 115
wolffd@0 116