|
Packit |
67cb25 |
#include <stdio.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_filter.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_rng.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_randist.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
main(void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = 1000; /* length of time series */
|
|
Packit |
67cb25 |
const size_t K = 25; /* window size */
|
|
Packit |
67cb25 |
const double t = 4.0; /* number of scale factors for outlier detection */
|
|
Packit |
67cb25 |
gsl_vector *x = gsl_vector_alloc(N); /* input vector */
|
|
Packit |
67cb25 |
gsl_vector *y = gsl_vector_alloc(N); /* output (filtered) vector */
|
|
Packit |
67cb25 |
gsl_vector *xmedian = gsl_vector_alloc(N); /* window medians */
|
|
Packit |
67cb25 |
gsl_vector *xsigma = gsl_vector_alloc(N); /* window scale estimates */
|
|
Packit |
67cb25 |
gsl_vector_int *ioutlier = gsl_vector_int_alloc(N); /* outlier detected? */
|
|
Packit |
67cb25 |
gsl_filter_impulse_workspace * w = gsl_filter_impulse_alloc(K);
|
|
Packit |
67cb25 |
gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
|
|
Packit |
67cb25 |
size_t noutlier;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* generate input signal */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = 10.0 * sin(2.0 * M_PI * i / (double) N);
|
|
Packit |
67cb25 |
double ei = gsl_ran_gaussian(r, 2.0);
|
|
Packit |
67cb25 |
double u = gsl_rng_uniform(r);
|
|
Packit |
67cb25 |
double outlier = (u < 0.01) ? 15.0*GSL_SIGN(ei) : 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set(x, i, xi + ei + outlier);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* apply impulse detection filter */
|
|
Packit |
67cb25 |
gsl_filter_impulse(GSL_FILTER_END_TRUNCATE, GSL_FILTER_SCALE_QN, t, x, y,
|
|
Packit |
67cb25 |
xmedian, xsigma, &noutlier, ioutlier, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* print results */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i);
|
|
Packit |
67cb25 |
double yi = gsl_vector_get(y, i);
|
|
Packit |
67cb25 |
double xmedi = gsl_vector_get(xmedian, i);
|
|
Packit |
67cb25 |
double xsigmai = gsl_vector_get(xsigma, i);
|
|
Packit |
67cb25 |
int outlier = gsl_vector_int_get(ioutlier, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
printf("%zu %f %f %f %f %d\n",
|
|
Packit |
67cb25 |
i,
|
|
Packit |
67cb25 |
xi,
|
|
Packit |
67cb25 |
yi,
|
|
Packit |
67cb25 |
xmedi + t * xsigmai,
|
|
Packit |
67cb25 |
xmedi - t * xsigmai,
|
|
Packit |
67cb25 |
outlier);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free(x);
|
|
Packit |
67cb25 |
gsl_vector_free(y);
|
|
Packit |
67cb25 |
gsl_vector_free(xmedian);
|
|
Packit |
67cb25 |
gsl_vector_free(xsigma);
|
|
Packit |
67cb25 |
gsl_vector_int_free(ioutlier);
|
|
Packit |
67cb25 |
gsl_filter_impulse_free(w);
|
|
Packit |
67cb25 |
gsl_rng_free(r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return 0;
|
|
Packit |
67cb25 |
}
|