1 /*************************************************************************/
2 /* */
3 /* Centre for Speech Technology Research */
4 /* University of Edinburgh, UK */
5 /* Copyright (c) 1997 */
6 /* All Rights Reserved. */
7 /* */
8 /* Permission is hereby granted, free of charge, to use and distribute */
9 /* this software and its documentation without restriction, including */
10 /* without limitation the rights to use, copy, modify, merge, publish, */
11 /* distribute, sublicense, and/or sell copies of this work, and to */
12 /* permit persons to whom this work is furnished to do so, subject to */
13 /* the following conditions: */
14 /* 1. The code must retain the above copyright notice, this list of */
15 /* conditions and the following disclaimer. */
16 /* 2. Any modifications must be clearly marked as such. */
17 /* 3. Original authors' names are not deleted. */
18 /* 4. The authors' names are not used to endorse or promote products */
19 /* derived from this software without specific prior written */
20 /* permission. */
21 /* */
22 /* THE UNIVERSITY OF EDINBURGH AND THE CONTRIBUTORS TO THIS WORK */
23 /* DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING */
24 /* ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT */
25 /* SHALL THE UNIVERSITY OF EDINBURGH NOR THE CONTRIBUTORS BE LIABLE */
26 /* FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES */
27 /* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN */
28 /* AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, */
29 /* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF */
30 /* THIS SOFTWARE. */
31 /* */
32 /*************************************************************************/
33 /* Author : Paul Bagshaw */
34 /* Date : 1993 */
35 /*-----------------------------------------------------------------------*/
36 /* */
37 /* The above copyright was given by Paul Bagshaw, he retains */
38 /* his original rights */
39 /* */
40 /*************************************************************************/
41 #include <cmath>
42 #include <cstdio>
43 #include "array_smoother.h"
44 #include "EST_cutils.h"
45
46 #define MAX_LEN 127
47
48 #define MODULE "array_smoother"
49
50 float median (int *counter, float valin, float valbuf[], int lmed, int mmed);
51 float hanning (int *counter, float valin, float valhan[], float win_coeff[],
52 struct Ms_Op *ms);
53 void mk_window_coeffs (int length, float win_coeff[]);
54
55 struct Ms_Op *default_ms_op(struct Ms_Op *ms);
56
array_smoother(float * p_array,int arraylen,struct Ms_Op * ms)57 void array_smoother (float *p_array, int arraylen, struct Ms_Op *ms)
58 {
59 int i, j, mid1, mid2 = 0, filler, nloops;
60 int C1, C2 = 0, C3 = 0, C4 = 0, c1, c2, c3, c4;
61 int delay, delx = 0, dely = 0;
62 int in = 0, out = 0;
63 float input, output;
64 float *inarray;
65 float xdel[2 * MAX_LEN - 2], ydel[2 * MAX_LEN - 2];
66 float medbuf1[MAX_LEN], medbuf2[MAX_LEN];
67 float hanbuf1[MAX_LEN], hanbuf2[MAX_LEN], win_coeffs[MAX_LEN];
68 float medval1, medval2, hanval1, hanval2, zatn;
69
70 inarray = new float[arraylen];
71 for (i = 0; i < arraylen; ++i)
72 inarray[i] = p_array[i];
73
74 if (ms == NULL)
75 {
76 ms = new Ms_Op;
77 default_ms_op(ms);
78 }
79
80 mk_window_coeffs (ms->window_length, win_coeffs);
81 /* determine the size and delay of each stage concerned */
82 mid1 = ms->first_median / 2;
83 C1 = delay = ms->first_median - 1;
84 if (ms->apply_hanning)
85 {
86 C2 = ms->window_length - 1;
87 delay = ms->first_median + ms->window_length - 2;
88 }
89 if (ms->smooth_double) {
90 mid2 = ms->second_median / 2;
91 C3 = ms->second_median - 1;
92 if (!ms->apply_hanning) {
93 delx = ms->first_median;
94 dely = ms->second_median;
95 }
96 else {
97 C4 = ms->window_length - 1;
98 delx = ms->first_median + ms->window_length - 1;
99 dely = ms->second_median + ms->window_length - 1;
100 }
101 delay = delx + dely - 2;
102 }
103 /* prepare for smoothing */
104 c1 = C1;
105 c2 = C2;
106 c3 = C3;
107 c4 = C4;
108 if (!ms->extrapolate) {
109 /* pad with breakers at the beginning */
110 for (i = 0; i < delay / 2; i++)
111 p_array[out++] = ms->breaker;
112 filler = 0;
113 nloops = arraylen;
114 }
115 else {
116 /* extrapolate by initialising filter with dummy breakers */
117 filler = delay / 2;
118 nloops = arraylen + delay;
119 }
120 /* smooth track element by track element */
121 for (j = 0; j < nloops; j++)
122 {
123 if (j < filler || j >= nloops - filler)
124 input = ms->breaker;
125 else
126 input = inarray[in++];
127
128 /* store input value if double smoothing */
129 if (ms->smooth_double) {
130 for (i = delx - 1; i > 0; i--)
131 xdel[i] = xdel[i - 1];
132 xdel[0] = input;
133 }
134 /* first median smoothing */
135
136 medval1 = median (&c1, input, medbuf1, ms->first_median, mid1);
137
138 if (c1 == -1)
139 {
140 output = medval1;
141 /* first hanning window (optional) */
142 if (ms->apply_hanning)
143 {
144 hanval1 = hanning (&c2, medval1, hanbuf1, win_coeffs, ms);
145 if (c2 == -1)
146 output = hanval1;
147 else
148 continue;
149 }
150 /* procedures for double smoothing (optional) */
151 if (ms->smooth_double)
152 {
153 /* compute rough component z(n) */
154 if (output != ms->breaker && xdel[delx - 1]
155 != ms->breaker)
156 zatn = xdel[delx - 1] - output;
157 else
158 zatn = ms->breaker;
159 /* store results of first smoothing */
160 for (i = dely - 1; i > 0; i--)
161 ydel[i] = ydel[i - 1];
162 ydel[0] = output;
163 /* second median smoothing */
164 medval2 = median (&c3, zatn, medbuf2,
165 ms->second_median, mid2);
166 if (c3 == -1)
167 {
168 output = medval2;
169 /* second hanning smoothing (optional) */
170 if (ms->apply_hanning) {
171 hanval2 = hanning (&c4, medval2, hanbuf2,
172 win_coeffs, ms);
173 if (c4 == -1)
174 output = hanval2;
175 else
176 continue;
177 }
178 if (output != ms->breaker && ydel[dely - 1]
179 != ms->breaker)
180 output += ydel[dely - 1];
181 else
182 output = ms->breaker;
183 }
184 else
185 continue;
186 }
187 /* write filtered result */
188 p_array[out++] = output;
189 }
190 }
191 if (!ms->extrapolate) /* pad with breakers at the end */
192 for (i = 0; i < delay / 2; i++)
193 p_array[out++] = ms->breaker;
194
195 delete inarray;
196 }
197
median(int * counter,float valin,float valbuf[],int lmed,int mmed)198 float median (int *counter, float valin, float valbuf[], int lmed, int mmed)
199 {
200 int i, j;
201 float tmp, filmed[MAX_LEN];
202
203 for (i = lmed - 1; i > 0; i--)
204 valbuf[i] = valbuf[i - 1];
205 valbuf[0] = valin;
206
207 if (*counter > 0)
208 {
209 (*counter)--;
210 return (0.0);
211 }
212 else
213 {
214 *counter = -1;
215
216 for (i = 0; i < lmed; i++)
217 filmed[i] = valbuf[i];
218
219 for (j = lmed - 1; j > 0; j--)
220 for (i = 0; i < j; i++)
221 if (filmed[i] > filmed[i + 1])
222 {
223 tmp = filmed[i + 1];
224 filmed[i + 1] = filmed[i];
225 filmed[i] = tmp;
226 }
227 return (filmed[mmed]);
228 }
229
230 }
231
232 #define TWO_PI 6.28318530717958647698
233
mk_window_coeffs(int length,float win_coeff[])234 void mk_window_coeffs (int length, float win_coeff[])
235 {
236 int i;
237 double x;
238
239 for (i = 0; i < length; i++) {
240 x = TWO_PI * (i + 1.0) / (length + 1.0);
241 win_coeff[i] = (1.0 - (float) cos (x)) / (length + 1.0);
242 }
243
244 }
245
hanning(int * counter,float valin,float valhan[],float win_coeff[],struct Ms_Op * par)246 float hanning (int *counter, float valin, float valhan[], float win_coeff[],
247 struct Ms_Op *par)
248 {
249 int i, j, k = 0;
250 float valout = 0.0, weight[MAX_LEN];
251
252 for (i = par->window_length - 1; i > 0; i--)
253 valhan[i] = valhan[i - 1];
254 valhan[0] = valin;
255 if (*counter > 0) {
256 (*counter)--;
257 return (0.0);
258 }
259 else {
260 *counter = -1;
261 for (i = 0; i < par->window_length; i++)
262 if (valhan[i] == par->breaker)
263 k++;
264 if (!k)
265 for (i = 0; i < par->window_length; i++)
266 valout += valhan[i] * win_coeff[i];
267 else if (k <= par->window_length / 2 && par->extrapolate) {
268 mk_window_coeffs (par->window_length - k, weight);
269 for (i = 0, j = 0; i < par->window_length; i++)
270 if (valhan[i] != par->breaker)
271 valout += valhan[i] * weight[j++];
272 }
273 else
274 valout = par->breaker;
275 return (valout);
276 }
277
278 }
279
initialise_parameters(struct Ms_Op * p_par)280 void initialise_parameters (struct Ms_Op *p_par)
281 {
282 p_par->smooth_double = 0;
283 p_par->apply_hanning = 0;
284 p_par->extrapolate = 0;
285 p_par->window_length = DEFAULT_WLEN;
286 p_par->first_median = DEFAULT_MED_1;
287 p_par->second_median = DEFAULT_MED_2;
288 return;
289 }
290
default_ms_op(struct Ms_Op * ms)291 struct Ms_Op *default_ms_op(struct Ms_Op *ms)
292 {
293 ms->smooth_double = FALSE;
294 ms->apply_hanning = TRUE;
295 ms->extrapolate = TRUE;
296 ms->first_median = 11;
297 ms->second_median = 1;
298 ms->window_length = 7;
299 ms->breaker = -1.0;
300 return (ms);
301 }
302