speech-tools/sigpr/pda/pcb_smoother.cc
2015-09-19 10:52:26 +02:00

302 lines
9.1 KiB
C++

/*************************************************************************/
/* */
/* Centre for Speech Technology Research */
/* University of Edinburgh, UK */
/* Copyright (c) 1997 */
/* All Rights Reserved. */
/* */
/* Permission is hereby granted, free of charge, to use and distribute */
/* this software and its documentation without restriction, including */
/* without limitation the rights to use, copy, modify, merge, publish, */
/* distribute, sublicense, and/or sell copies of this work, and to */
/* permit persons to whom this work is furnished to do so, subject to */
/* the following conditions: */
/* 1. The code must retain the above copyright notice, this list of */
/* conditions and the following disclaimer. */
/* 2. Any modifications must be clearly marked as such. */
/* 3. Original authors' names are not deleted. */
/* 4. The authors' names are not used to endorse or promote products */
/* derived from this software without specific prior written */
/* permission. */
/* */
/* THE UNIVERSITY OF EDINBURGH AND THE CONTRIBUTORS TO THIS WORK */
/* DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING */
/* ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT */
/* SHALL THE UNIVERSITY OF EDINBURGH NOR THE CONTRIBUTORS BE LIABLE */
/* FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES */
/* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN */
/* AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, */
/* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF */
/* THIS SOFTWARE. */
/* */
/*************************************************************************/
/* Author : Paul Bagshaw */
/* Date : 1993 */
/*-----------------------------------------------------------------------*/
/* */
/* The above copyright was given by Paul Bagshaw, he retains */
/* his original rights */
/* */
/*************************************************************************/
#include <cmath>
#include <cstdio>
#include "array_smoother.h"
#include "EST_cutils.h"
#define MAX_LEN 127
#define MODULE "array_smoother"
float median (int *counter, float valin, float valbuf[], int lmed, int mmed);
float hanning (int *counter, float valin, float valhan[], float win_coeff[],
struct Ms_Op *ms);
void mk_window_coeffs (int length, float win_coeff[]);
struct Ms_Op *default_ms_op(struct Ms_Op *ms);
void array_smoother (float *p_array, int arraylen, struct Ms_Op *ms)
{
int i, j, mid1, mid2 = 0, filler, nloops;
int C1, C2 = 0, C3 = 0, C4 = 0, c1, c2, c3, c4;
int delay, delx = 0, dely = 0;
int in = 0, out = 0;
float input, output;
float *inarray;
float xdel[2 * MAX_LEN - 2], ydel[2 * MAX_LEN - 2];
float medbuf1[MAX_LEN], medbuf2[MAX_LEN];
float hanbuf1[MAX_LEN], hanbuf2[MAX_LEN], win_coeffs[MAX_LEN];
float medval1, medval2, hanval1, hanval2, zatn;
inarray = new float[arraylen];
for (i = 0; i < arraylen; ++i)
inarray[i] = p_array[i];
if (ms == NULL)
{
ms = new Ms_Op;
default_ms_op(ms);
}
mk_window_coeffs (ms->window_length, win_coeffs);
/* determine the size and delay of each stage concerned */
mid1 = ms->first_median / 2;
C1 = delay = ms->first_median - 1;
if (ms->apply_hanning)
{
C2 = ms->window_length - 1;
delay = ms->first_median + ms->window_length - 2;
}
if (ms->smooth_double) {
mid2 = ms->second_median / 2;
C3 = ms->second_median - 1;
if (!ms->apply_hanning) {
delx = ms->first_median;
dely = ms->second_median;
}
else {
C4 = ms->window_length - 1;
delx = ms->first_median + ms->window_length - 1;
dely = ms->second_median + ms->window_length - 1;
}
delay = delx + dely - 2;
}
/* prepare for smoothing */
c1 = C1;
c2 = C2;
c3 = C3;
c4 = C4;
if (!ms->extrapolate) {
/* pad with breakers at the beginning */
for (i = 0; i < delay / 2; i++)
p_array[out++] = ms->breaker;
filler = 0;
nloops = arraylen;
}
else {
/* extrapolate by initialising filter with dummy breakers */
filler = delay / 2;
nloops = arraylen + delay;
}
/* smooth track element by track element */
for (j = 0; j < nloops; j++)
{
if (j < filler || j >= nloops - filler)
input = ms->breaker;
else
input = inarray[in++];
/* store input value if double smoothing */
if (ms->smooth_double) {
for (i = delx - 1; i > 0; i--)
xdel[i] = xdel[i - 1];
xdel[0] = input;
}
/* first median smoothing */
medval1 = median (&c1, input, medbuf1, ms->first_median, mid1);
if (c1 == -1)
{
output = medval1;
/* first hanning window (optional) */
if (ms->apply_hanning)
{
hanval1 = hanning (&c2, medval1, hanbuf1, win_coeffs, ms);
if (c2 == -1)
output = hanval1;
else
continue;
}
/* procedures for double smoothing (optional) */
if (ms->smooth_double)
{
/* compute rough component z(n) */
if (output != ms->breaker && xdel[delx - 1]
!= ms->breaker)
zatn = xdel[delx - 1] - output;
else
zatn = ms->breaker;
/* store results of first smoothing */
for (i = dely - 1; i > 0; i--)
ydel[i] = ydel[i - 1];
ydel[0] = output;
/* second median smoothing */
medval2 = median (&c3, zatn, medbuf2,
ms->second_median, mid2);
if (c3 == -1)
{
output = medval2;
/* second hanning smoothing (optional) */
if (ms->apply_hanning) {
hanval2 = hanning (&c4, medval2, hanbuf2,
win_coeffs, ms);
if (c4 == -1)
output = hanval2;
else
continue;
}
if (output != ms->breaker && ydel[dely - 1]
!= ms->breaker)
output += ydel[dely - 1];
else
output = ms->breaker;
}
else
continue;
}
/* write filtered result */
p_array[out++] = output;
}
}
if (!ms->extrapolate) /* pad with breakers at the end */
for (i = 0; i < delay / 2; i++)
p_array[out++] = ms->breaker;
delete inarray;
}
float median (int *counter, float valin, float valbuf[], int lmed, int mmed)
{
int i, j;
float tmp, filmed[MAX_LEN];
for (i = lmed - 1; i > 0; i--)
valbuf[i] = valbuf[i - 1];
valbuf[0] = valin;
if (*counter > 0)
{
(*counter)--;
return (0.0);
}
else
{
*counter = -1;
for (i = 0; i < lmed; i++)
filmed[i] = valbuf[i];
for (j = lmed - 1; j > 0; j--)
for (i = 0; i < j; i++)
if (filmed[i] > filmed[i + 1])
{
tmp = filmed[i + 1];
filmed[i + 1] = filmed[i];
filmed[i] = tmp;
}
return (filmed[mmed]);
}
}
#define TWO_PI 6.28318530717958647698
void mk_window_coeffs (int length, float win_coeff[])
{
int i;
double x;
for (i = 0; i < length; i++) {
x = TWO_PI * (i + 1.0) / (length + 1.0);
win_coeff[i] = (1.0 - (float) cos (x)) / (length + 1.0);
}
}
float hanning (int *counter, float valin, float valhan[], float win_coeff[],
struct Ms_Op *par)
{
int i, j, k = 0;
float valout = 0.0, weight[MAX_LEN];
for (i = par->window_length - 1; i > 0; i--)
valhan[i] = valhan[i - 1];
valhan[0] = valin;
if (*counter > 0) {
(*counter)--;
return (0.0);
}
else {
*counter = -1;
for (i = 0; i < par->window_length; i++)
if (valhan[i] == par->breaker)
k++;
if (!k)
for (i = 0; i < par->window_length; i++)
valout += valhan[i] * win_coeff[i];
else if (k <= par->window_length / 2 && par->extrapolate) {
mk_window_coeffs (par->window_length - k, weight);
for (i = 0, j = 0; i < par->window_length; i++)
if (valhan[i] != par->breaker)
valout += valhan[i] * weight[j++];
}
else
valout = par->breaker;
return (valout);
}
}
void initialise_parameters (struct Ms_Op *p_par)
{
p_par->smooth_double = 0;
p_par->apply_hanning = 0;
p_par->extrapolate = 0;
p_par->window_length = DEFAULT_WLEN;
p_par->first_median = DEFAULT_MED_1;
p_par->second_median = DEFAULT_MED_2;
return;
}
struct Ms_Op *default_ms_op(struct Ms_Op *ms)
{
ms->smooth_double = FALSE;
ms->apply_hanning = TRUE;
ms->extrapolate = TRUE;
ms->first_median = 11;
ms->second_median = 1;
ms->window_length = 7;
ms->breaker = -1.0;
return (ms);
}