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