libhat-0.3/ 0000755 0316553 0005316 00000000000 13333033026 012146 5 ustar micarla def-cg libhat-0.3/src/ 0000755 0316553 0005316 00000000000 13333033025 012734 5 ustar micarla def-cg libhat-0.3/src/fa.h 0000666 0316553 0005316 00000000240 13153300017 013471 0 ustar micarla def-cg double get_principal_line_fast(double *x, int count, double q_min, double q_max);
double get_principal_line(double *x, int count, double q_min, double q_max);
libhat-0.3/src/fa.c 0000644 0316553 0005316 00000010401 13256664333 013502 0 ustar micarla def-cg /*
This file is part of hat.
hat is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
hat is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with hat. If not, see .
*/
#include
#include
#include
#include
#include
#include "fa_utils.h"
double
get_amplitude(double *buf, int count, double q)
{
int j;
double amp;
double complex z, z_aux;
double complex twiddle, twiddle_0;
twiddle_0 = cexp(2. * M_PI * I * (double complex)(q));
z = 0.;
twiddle = 1.;
for (j = 0; j < count; j++) {
z_aux = buf[j];
if (isnan(z_aux))
continue;
z_aux *= 1. - (double complex)cos(2. * M_PI * (double)j / (double)(count - 1));
z_aux *= twiddle;
z += z_aux;
twiddle *= twiddle_0;
}
amp = cabs(z);
return amp;
}
double
get_principal_line_guess(double *x, int count, double q_min, double q_max)
{
int i;
int max_idx;
double max, max_p1, max_m1;
double amp, amp_previous;
double q;
max_idx = -1.;
max = 0.;
max_p1 = NAN;
max_m1 = NAN;
amp_previous = NAN;
q_min *= count;
q_max *= count;
if (q_min < 1)
q_min = 1;
if (q_max > count / 2 - 1)
q_max = count / 2 - 1;
for (i = q_min; i <= q_max; i++) {
q = (double)i / (double)count;
amp = get_amplitude(x, count, q);
if (amp > max) {
max = amp;
max_idx = i;
max_m1 = amp_previous; /* Point before the maximum */
}
amp_previous = amp;
if (i == max_idx + 1) /* Point after the maximum */
max_p1 = amp;
}
if (max_idx == q_min || max_idx == q_max || max_idx == -1)
return NAN;
q = max_idx / (double)count;
return q;
}
struct Data {
double *buf; /* Data to be fitted */
int turns; /* Number of turns */
};
double
test_amplitude(double q, void *param)
{
int count;
double *buf;
double amp;
buf = ((struct Data *)param)->buf;
count = ((struct Data *)param)->turns;
amp = get_amplitude(buf, count, q);
return -1. * amp;
}
double
get_principal_line_laskar(double *data, int turns, double guess)
{
int status;
int iter;
gsl_min_fminimizer *minimizer;
double lower_bound, upper_bound;
gsl_function f;
struct Data fit_data;
fit_data.buf = data;
fit_data.turns = turns;
lower_bound = guess - 1. / (double)(turns);
upper_bound = guess + 1. / (double)(turns);
if (test_amplitude(guess, &fit_data) >= test_amplitude(lower_bound, &fit_data) ||
test_amplitude(guess, &fit_data) >= test_amplitude(upper_bound, &fit_data)) {
return NAN;
}
f.function = &test_amplitude;
f.params = &fit_data;
minimizer = gsl_min_fminimizer_alloc(gsl_min_fminimizer_brent);
gsl_min_fminimizer_set(minimizer, &f, guess,
lower_bound, upper_bound);
iter = 0;
do {
iter++;
status = gsl_min_fminimizer_iterate(minimizer);
if (status) {
guess = NAN;
break;
}
guess = gsl_min_fminimizer_x_minimum(minimizer);
lower_bound = gsl_min_fminimizer_x_lower(minimizer);
upper_bound = gsl_min_fminimizer_x_upper(minimizer);
status = gsl_min_test_interval(lower_bound, upper_bound, 1e-6, 1e-6);
} while (status == GSL_CONTINUE && iter < 10);
gsl_min_fminimizer_free(minimizer);
return guess;
}
double
get_principal_line(double *x, int count, double q_min, double q_max)
{
double q;
q = get_principal_line_guess(x, count, q_min, q_max);
if (isnan(q))
return NAN;
q = get_principal_line_laskar(x, count, q);
return q;
}
libhat-0.3/src/damping.h 0000666 0316553 0005316 00000000110 13150575037 014533 0 ustar micarla def-cg double fit_damping(double *buf, int turns, double freq, double phase);
libhat-0.3/src/fa_utils.h 0000666 0316553 0005316 00000001003 13261427653 014727 0 ustar micarla def-cg double complex ortho_projection(double *buf, int count, double q);
double complex ortho_projection_damped(double *buf, int count,
double q, double damping);
void init_hanning(int count);
void window_data(double *data, double *windowed_data, int turns);
void undamp(double *data, int turns, double k);
void purge_line(double *data, int count, double q,
double complex z, double damping);
void subtract_and_undamp(double *data, int turns, double *guess);
libhat-0.3/src/damping.c 0000666 0316553 0005316 00000003350 13333032002 014516 0 ustar micarla def-cg /*
This file is part of hat.
hat is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
hat is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with hat. If not, see .
*/
#include
#include "amoeba.h"
struct Data {
double *y; /* Data to be fitted */
int turns; /* Number of turns */
double omega;
double phi;
};
double
f(double *x, void *data)
{
double *y;
double aux;
double omega, phi, amp, tau;
double residual;
int turns, t;
turns = ((struct Data *)data)->turns;
y = ((struct Data *)data)->y;
omega = ((struct Data *)data)->omega;
phi = ((struct Data *)data)->phi;
amp = x[0];
tau = x[1];
residual = 0;
for (t = 0; t < turns; t++) {
aux = amp * cos(omega * (double)t + phi) * exp((double)t * tau) - y[t];
residual += aux * aux;
}
return residual;
}
double
fit_damping(double *buf, int turns, double freq, double phase)
{
struct Data data;
double simplex[3][2] = {{0., 0.}, {1., 0.}, {0., -1./100.}};
int ilo;
data.turns = turns;
data.y = buf;
data.omega = freq * 2. * M_PI;
data.phi = phase * 2. * M_PI;
ilo = amoeba(simplex, 1e-6, f, (void *)&data);
return -1. / simplex[ilo][1];
}
libhat-0.3/src/pyhat.c 0000644 0316553 0005316 00000014233 13261427507 014244 0 ustar micarla def-cg /*
This file is part of hat.
hat is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
hat is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with hat. If not, see .
*/
#include
#include
#include
#include
#include "fa_utils.h"
#include "damping.h"
#include "fa.h"
static PyObject*
find_dominant_line_naff_wrapper(PyObject *dummy, PyObject *args)
{
PyObject *arg1 = NULL, *arg2 = NULL;
PyObject *data = NULL, *freq_bound = NULL;
double q;
npy_intp *count;
double *data_buf, *freq_bound_buf;
if (!PyArg_ParseTuple(args, "OO", &arg1, &arg2))
return NULL;
data = PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_IN_ARRAY);
freq_bound = PyArray_FROM_OTF(arg2, NPY_DOUBLE, NPY_IN_ARRAY);
if (data == NULL || freq_bound == NULL)
return NULL;
count = PyArray_DIMS(data); //number of dimensions
data_buf = (double *)PyArray_DATA(data);
freq_bound_buf = (double *)PyArray_DATA(freq_bound);
q = get_principal_line(data_buf, (int)count[0],
freq_bound_buf[0], freq_bound_buf[1]);
Py_DECREF(data);
Py_DECREF(freq_bound);
return Py_BuildValue("d", q);
}
static PyObject*
eval_line_naff_wrapper(PyObject *dummy, PyObject *args)
{
PyObject *arg1 = NULL, *arg2 = NULL;
PyObject *data = NULL;
PyObject *z_obj;
double q;
double complex z;
double *data_buf;
npy_intp *count;
if (!PyArg_ParseTuple(args, "OO", &arg1, &arg2))
return NULL;
q = PyFloat_AsDouble(arg2);
data = PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_IN_ARRAY);
if (data == NULL)
return NULL;
count = PyArray_DIMS(data); //number of dimensions
data_buf = (double *)PyArray_DATA(data);
z = ortho_projection(data_buf, (int)count[0], q);
z_obj = PyComplex_FromDoubles(creal(z), cimag(z));
Py_DECREF(data);
return z_obj;
}
static PyObject*
eval_damping_wrapper(PyObject *dummy, PyObject *args)
{
PyObject *arg1 = NULL, *arg2 = NULL, *arg3 = NULL;
PyObject *data = NULL;
double q, phase, damping;
double *data_buf;
npy_intp *count;
if (!PyArg_ParseTuple(args, "OOO", &arg1, &arg2, &arg3))
return NULL;
q = PyFloat_AsDouble(arg2);
phase = PyFloat_AsDouble(arg3);
data = PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_IN_ARRAY);
if (data == NULL)
return NULL;
count = PyArray_DIMS(data); //number of dimensions
data_buf = (double *)PyArray_DATA(data);
damping = fit_damping(data_buf, (int)count[0], q, phase);
Py_DECREF(data);
return Py_BuildValue("d", damping);
}
static PyObject*
eval_damped_line_naff_wrapper(PyObject *dummy, PyObject *args)
{
PyObject *arg1 = NULL, *arg2 = NULL, *arg3 = NULL;
PyObject *data = NULL;
double q, damping;
// double re, im;
double complex z;
PyObject *z_obj;
double *data_buf;
npy_intp *count;
if (!PyArg_ParseTuple(args, "OOO", &arg1, &arg2, &arg3))
return NULL;
q = PyFloat_AsDouble(arg2);
damping = PyFloat_AsDouble(arg3);
data = PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_IN_ARRAY);
if (data == NULL)
return NULL;
count = PyArray_DIMS(data); //number of dimensions
data_buf = (double *)PyArray_DATA(data);
z = ortho_projection_damped(data_buf, (int)count[0], q, damping);
z_obj = PyComplex_FromDoubles(creal(z), cimag(z));
Py_DECREF(data);
return z_obj;
//Py_DECREF(data);
//return Py_BuildValue("dd", re, im);
}
static PyObject*
purge_line_wrapper(PyObject *dummy, PyObject *args)
{
PyObject *arg1 = NULL, *arg2 = NULL;
PyObject *arg3 = NULL, *arg4 = NULL;
PyObject *data = NULL;
Py_complex z_obj;
double q, damping;
double complex z;
double *data_buf;
npy_intp *count;
if (!PyArg_ParseTuple(args, "OOOO", &arg1, &arg2, &arg3, &arg4))
return NULL;
q = PyFloat_AsDouble(arg2);
z_obj = PyComplex_AsCComplex(arg3);
damping = PyFloat_AsDouble(arg4);
data = PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_IN_ARRAY);
if (data == NULL)
return NULL;
count = PyArray_DIMS(data); //number of dimensions
data_buf = (double *)PyArray_DATA(data);
z = z_obj.real + I * z_obj.imag;
purge_line(data_buf, (int)count[0], q, z, damping);
Py_DECREF(data);
Py_RETURN_NONE;
}
static struct PyMethodDef methods[] = {
{"find_dominant_line_naff",
find_dominant_line_naff_wrapper, METH_VARARGS,
"Return the dominant spectral line in the"
"given bound, using a Laskar like method"},
{"eval_line_naff",
eval_line_naff_wrapper, METH_VARARGS,
"Evaluate amplitude and phase of a spectral line"},
{"eval_damped_line_naff",
eval_damped_line_naff_wrapper, METH_VARARGS,
"Evaluate amplitude and phase of a spectral line in"
"presence of damping"},
{"eval_damping",
eval_damping_wrapper, METH_VARARGS,
"Evaluate damping"},
{"purge_line",
purge_line_wrapper, METH_VARARGS,
"Remove a spectral line from the signal"},
{NULL, NULL, 0, NULL}
};
PyMODINIT_FUNC
inithat(void)
{
char *doc = "HAT is my personal attempt at an Harmonic Analysis Toolkit\n"
"Copyright Michele Carla' 2018";
(void)Py_InitModule3("hat", methods, doc);
import_array();
}
libhat-0.3/src/fa_utils.c 0000666 0316553 0005316 00000004042 13316705410 014717 0 ustar micarla def-cg /*
This file is part of hat.
hat is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
hat is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with hat. If not, see .
*/
#include
#include
double complex
ortho_projection(double *buf, int count, double q)
{
int i;
double x;
double complex z;
q *= -2. * M_PI;
z = 0;
for (i = 0; i < count; i++) {
x = buf[i];
x *= 1. - cos(2. * M_PI * (double)i / (double)count);
z += x * cexp(I * (double complex)q * (double complex)i);
}
z *= 2. / (double)count;
return z;
}
double complex
ortho_projection_damped(double *buf, int count, double q, double damping)
{
int i;
double x;
double complex z;
damping = -1. * damping;
q *= -2. * M_PI;
z = 0;
for (i = 0; i < count; i++) {
x = buf[i];
x *= exp((double)i / damping);
x *= 1. - cos(2. * M_PI * (double)i / (double)count);
z += x * cexp(I * (double complex)q * (double complex)i);
}
z *= 2. / (double)count;
return z;
}
void
undamp(double *data, int turns, double k)
{
int t;
for (t = 0; t < turns; t++) {
data[t] *= exp(-k * t);
}
}
void
purge_line(double *data, int count, double q,
double complex z, double damping)
{
int i;
double aux;
q *= 2. * M_PI;
for (i = 0; i < count; i++) {
aux = creal(z * cexp(I * q * (double)i));
aux *= exp((double)i / damping);
data[i] -= aux;
}
}
libhat-0.3/src/pyhat.o 0000644 0316553 0005316 00000022400 13333033025 014237 0 ustar micarla def-cg ELF > @ @ UHH H= HEHE H} u
HEH5 H HEH} u=HEH HPHEHHEH HuHEH@H@0HUHи HEHPH H9uHEH H HEH HPHEHHEH HuHEH@H@0HUHHEH HPHEHHEH HuHEH@H@0HUHH Hu
H H = t=H H ЉH H Ѻ H5 HǸ H H H Ѓw@H H H ЉH H Ѻ H5 HǸ pH H H ЉE} u%H H H5 HǸ 0}t%H H H5 HǸ UHSHXH}HuHE HE HE HE HMHUHEH5 HǸ u
H H( H HH Hh H HHEA A HHEH H( H HH Hh H HHEA A HHEH} tH} u
HEH@ HEHEH@HEHEH@HEHEHHE HEH HEH EHEH HPHEHHEH HuHEH@H@0HUHHEH HPHEHHEH HuHEH@H@0HUHEH= HX[UHSHxH}HuHE HE HE HMHUHEH5 HǸ u
HEH EH H( H HH Hh H HHEA A HHEH} u
HEH@ HEHEH@HEHEH EHEH EMHEHEHEHEME HEHEH HPHEHHEH HuHEH@H@0HUHHEHx[UHSHhH}HuHE HE HE HE H]HMHUHEIH5 HǸ u
HEH EHEH EH H( H HH Hh H HHEA A HHEH} u
HEH@ HEHEH@HEHEH MEHEH EHEH HPHEHHEH HuHEH@H@0HUHEH= Hh[UHSH H}HuHE HE HE HE H]HMHUHEIH5 HǸ u
, HEH EHEH EH H( H HH Hh H HHEA A HHEH} u
HEH@ HEHEH@HEHEH MEHEH pxHpHEHxHEME HEHEH HPHEHHEH HuHEH@H@0HUHHEHĈ [UHATSH HxHpHE HE HE HE HE HuH]HMHUHpIIH5 HǸ u
HEH EHEH `hH`HEHhHEHEH EH H( H HH Hh H HHEA A HHEH} u
H HEH@ HEHEH@HEHEIĻ HEH¸ HHHHHH fW PXHPHXLHHHHHHHHHHHXXMEHEH ]MUEHEH HEH HPHEHHEH HuHEH@H@0HUHH H HPH HH Hİ [A\UHHH HEHEA HH5 H= y H H H5 H numpy.core.multiarray _ARRAY_API module compiled against ABI version %x but this version of numpy is %x module compiled against API version %x but this version of numpy is %x FATAL: module compiled as unknown endian FATAL: module compiled as little endian, but detected different endianness at runtime OO d OOO OOOO find_dominant_line_naff Return the dominant spectral line in thegiven bound, using a Laskar like method eval_line_naff Evaluate amplitude and phase of a spectral line eval_damped_line_naff Evaluate amplitude and phase of a spectral line inpresence of damping eval_damping Evaluate damping purge_line Remove a spectral line from the signal HAT is my personal attempt at an Harmonic Analysis Toolkit
Copyright Michele Carla' 2018 hat numpy.core.multiarray failed to import ? GCC: (GNU) 4.4.7 20120313 (Red Hat 4.4.7-23) zR x H AC
C <