Mercurial > hg > ltpda

comparison src/ltpda_smoother/ltpda_smoother.c @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 #include <math.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include <mex.h>
6
7 #include "version.h"
8 #include "ltpda_smoother.h"
9 #include "quick.c"
10
11 #define DEBUG 0
12 /*
13 * A smoother function that uses various methods to smooth data.
14 *
15 * M Hewitson 28-08-01
16 *
17 * $Id: ltpda_smoother.c,v 1.4 2008/09/25 08:16:16 hewitson Exp $
18 */
19
20
21 /*
22 function sy = ltpda_smoother(y, bw, ol, method);
23
24 */
25 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
26 {
27 int status;
28 int buff_len;
29
30 /* outputs */
31 double *outputPtr;
32 double *nxx;
33 int a;
34
35 /* inputs */
36 double *xx;
37 int bw;
38 double ol;
39 int nx;
40 char *method;
41
42 /* parse input functions */
43
44 if( (nrhs == 0) || (nlhs == 0) )
45 {
46 print_usage(VERSION);
47 }
48
49 if( (nrhs == 4) && (nlhs == 1) )/* let's go */
50 {
51 /*----------------- set inputs*/
52 xx = mxGetPr(prhs[0]);
53 bw = (int)floor(mxGetScalar(prhs[1]));
54 ol = mxGetScalar(prhs[2]);
55 nx = mxGetNumberOfElements(prhs[0]);
56
57 /* Read the method*/
58 buff_len = (mxGetM(prhs[3]) * mxGetN(prhs[3])) + 1;
59 method = mxCalloc(buff_len, sizeof(char));
60 status = mxGetString(prhs[3], method, buff_len);
61 if(status != 0)
62 mexWarnMsgTxt("Not enough space. String is truncated.");
63
64 /* create output vector*/
65 nxx = (double*)mxCalloc(nx, sizeof(double));
66
67 /* nfest*/
68 smooth(nxx, xx, nx, bw, ol, method);
69
70 /* output noise-floor vector*/
71 plhs[0] = mxCreateDoubleMatrix(1,nx,mxREAL);
72 outputPtr = mxGetPr(plhs[0]);
73 for(a=0; a<nx; a++)
74 outputPtr[a] = nxx[a];
75
76 mxFree(nxx);
77 }
78 else
79 {
80 print_usage(VERSION);
81 }
82 }
83
84 void print_usage(char *version)
85 {
86 mexPrintf("ltpda_smoother version %s\n", version);
87 mexPrintf(" usage: function sy = ltpda_smoother(y, bw, ol, method); \n");
88 mexErrMsgTxt("### incorrect usage");
89 }
90 /*
91 * Smoother
92 */
93 int smooth(double *nxx, double *xx, int nx, int bw, double ol, char *method)
94 {
95 int k, j, idx;
96 int hbw;
97 double *seg;
98
99 seg = (double*)mxCalloc(bw+1, sizeof(double));
100
101 /* go through each element */
102 for(k=0; k<nx; k++)
103 {
104 /* get segment*/
105 for(j=0; j<=bw; j++)
106 {
107 idx = k+j-bw/2;
108 if (idx<0)
109 idx = 0;
110 if (idx>=nx)
111 idx = nx-1;
112 seg[j] = xx[idx];
113 }
114 /* sort segment*/
115 quickSort(seg, bw-1);
116 /* stop index*/
117 idx = (int)floor(ol*bw);
118 /* Which smoothing method? */
119 if (strcmp(method, "median")==0)
120 {
121 /* make median estimate of selected samples */
122 nxx[k] = median(seg, idx);
123 }
124 else if (strcmp(method, "mean")==0)
125 {
126 /* make mean estimate of selected samples */
127 nxx[k] = mean(seg, idx);
128 }
129 else if (strcmp(method, "min")==0)
130 {
131 /* make min estimate of selected samples */
132 nxx[k] = smin(seg, idx);
133 }
134 else if (strcmp(method, "max")==0)
135 {
136 /* make max estimate of selected samples */
137 nxx[k] = smax(seg, idx);
138 }
139 }
140
141 /* Clean up */
142 mxFree(seg);
143
144 return 0;
145 }
146
147
148 double smin(double *x, int nx)
149 {
150 double m;
151 int j;
152
153 m = x[0];
154
155 for (j=1; j<nx; j++)
156 if (x[j] < m)
157 m = x[j];
158
159 return m;
160 }
161
162 double smax(double *x, int nx)
163 {
164 double m;
165 int j;
166
167 m = x[0];
168
169 for (j=1; j<nx; j++)
170 if (x[j] > m)
171 m = x[j];
172
173 return m;
174 }
175
176 double mean(double *x, int nx)
177 {
178 int j;
179 double sx = 0.0;
180
181 for (j=0; j<nx; j++)
182 sx += x[j];
183
184 return sx/nx;
185 }
186
187 double median(double *x, int nx)
188 {
189 double m;
190
191 if (nx%2 == 0) /* even*/
192 m = (x[nx/2] + x[nx/2 -1])/2.0;
193 else
194 m = x[(nx-1)/2];
195
196 return m;
197 }