I attached an "translation" of your code (wavelet2.py) for the numerical
part. It is supposed to do no transforms to the arrays for the wavelet
transform. This is only done for the real / halfcomplex fourier
transform so that the arrays are of complex or real type and the data is
reordered to make it more pythonic.
As far as I understand wavelets are real to real transforms. So the
wrapper is only checking if lengthes are valid, and if a workspace is
required to be allocated.
In CVS was again a bug for a precompile statement. (Which is fixed now).
If you use pygsl 0.3.2 I assume that this bug does not exist. If so, try
the CVS version. If you get an GSL_SANITY error, please put the attached
file core.c to src/transform/core.c.
It would be nice if you could extend the attached example to a real one,
(perhaps generate a test signal and transform it, so that one can see
what the wavelets are doing).
Hoping to hear from you!
Sincerely yours
Pierre
basically what I need to do is apply the wavelet transform (this is
obtain the wavelet coefficients) to a discrete signal obtained from a
sampled vibration signal, the results of this transformation will be
applied to an algorithm for the compute of the energies of the bands of
frecuency and then obtain a codebook that will be used to obtain a group
of secuencies to be applied to an artificial intelligence software, but
that is another story. I already calculated the wavelet transform of my
vibration signal, I use the following program
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_wavelet.h>
int main (int argc, char **argv)
{
int i, n = 512;
double *data = malloc(n * sizeof(double));
FILE *f = fopen(argv[1],"r");
for (i=0; i<n; i++)
{
fscanf(f, "%lg", &data[i]);
}
fclose (f);
{
gsl_wavelet *w = gsl_wavelet_alloc(gsl_wavelet_daubechies, 4);
gsl_wavelet_workspace *work = gsl_wavelet_workspace_alloc(n);
gsl_wavelet_transform_forward(w,data,1,n,work);
}
FILE *g = fopen("wavcoef.dat","w");
for (i=0; i<n; i++)
{
fprintf(g, "%g\n", data[i]);
}
fclose (g);
return 0;
, as you can see this is very simple, I just need to know if I can do
this with PyGSL and if the output of the python library is the same as
the GSL C code, I mean the output array with the wavelet coeficcients
has the same order, I think that I can use the example that you include
with PyGSL, if this works everything is ok.
I had no check the Python code yet but I have to do in this days. So if
your have some other advice I will really thank you.
Thank you very much.
On 10/2/06, * Pierre SCHNIZER* <[EMAIL PROTECTED]
<mailto:[EMAIL PROTECTED]>> wrote:
Fernando Vargas wrote:
#!/usr/bin/env python
#
# Example for a wavelet analysis
# Authors: Fernando Vargas, Pierre Schnizer
# Date: October 2006
#
# To use the same numerical array package as pygsl.
# I only use it here, so not to import a package not installed on the user
# side. Normaly I directly import the array package I use for all my scripts
# and programs.
# There are subtle differences between the packages, thus I would be careful
# which I use in my scripts.
from pygsl import _numobj as numx
from pygsl import wavelet
def transform(data):
n = len(data)
w = wavelet.daubechies(4)
# If no workspace is given, the wrapper will allocate one itself for the
# transform. I am not sure how much of an performance impact that is.
result = w.transform_forward(data);
# If you perfer you could have the workspace and the transform stored
# somewhere, and then perhaps your code is faster
# work = wavelet.wavelet_workspace(n)
# result = w.transform_forward(data, work)
return result
def run():
# Here should go some meaningful data sample to test the transform.
# Can you write one?
n = 512
data = numx.zeros((n,))
result = transform(data)
print result
if __name__ == '__main__':
run()
#ifdef _PyGSL_HAS_WAVELET
#include "wavelet.h"
#endif
/*
* Checks existing objects if they are of proper type, and if so check that
* they are big enough. If not allocate space of approbriate size.
*/
static int
PyGSL_transform_helpers_alloc(PyObject *s_o, PyObject *t_o, struct _pygsl_transform_help_rf_s * h, int n)
{
int check;
FUNC_MESS_BEGIN();
/* Set everything to zero first! */
h->free_space = 0;
h->free_table = 0;
h->table = NULL;
h->space = NULL;
if(h->func == 0)
GSL_ERROR("Functions not set!", GSL_EFAULT);
if(n <= 0)
GSL_ERROR("n<=0!", GSL_ESANITY);
DEBUG_MESS(3, "Allocating/Checking space for %d elements", n);
if(h->func->space_type != NOSPACE && s_o){
if(PyGSL_transform_space_check(s_o) && ((PyGSL_transform_space * )s_o)->type == h->func->space_type){
check = PyGSL_transform_space_get_n((PyGSL_transform_space * )s_o);
if(check == -1 || check < n)
GSL_ERROR("Work Space not big enough!", GSL_EINVAL);
h->space = ((PyGSL_transform_space * )s_o) ->space.v;
} else {
PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__ - 4);
GSL_ERROR("Need a pygsl transform space of proper type!", GSL_EINVAL);
}
}
if(h->func->space_type != NOSPACE && t_o){
if(PyGSL_transform_space_check(t_o) && ((PyGSL_transform_space * )t_o)->type == h->func->table_type){
check = PyGSL_transform_space_get_n((PyGSL_transform_space * )s_o);
if(check == -1 || check < n)
GSL_ERROR("Wave table not big enough!", GSL_EINVAL);
h->table = ((PyGSL_transform_space * )t_o) ->space.v;
} else {
PyGSL_add_traceback(module, filename, __FUNCTION__, __LINE__ - 4);
GSL_ERROR("Need a pygsl transform wave table of proper type!", GSL_EINVAL);
}
}
/* Check for the approbriate type and initialise it!*/
if((h->space) == NULL || (h->table) == NULL){
DEBUG_MESS(3, "func %p alloc table %p alloc space %p, space %p, table %p",
h->func, (void *) h->func->table_alloc,
(void *) h->func->space_alloc, h->space, h->table);
/* Store if I need to free these arrays */
if((h->space) == NULL && h->func->space_type != NOSPACE){
h->space = h->func->space_alloc(n);
h->free_space = 1;
}
if((h->table) == NULL && h->func->table_type != NOSPACE){
h->table = h->func->table_alloc(n);
h->free_table = 1;
}
if((h->space == NULL && h->func->space_type != NOSPACE) ||
(h->table == NULL && h->func->table_type != NOSPACE)){
return GSL_ENOMEM;
}
DEBUG_MESS(3, "Allocated space @ %p table @ %p", h->space, h->table);
}
FUNC_MESS_END();
return GSL_SUCCESS;
}
#define PyGSL_TRANSFORM_HELPERS_FREE(helpers) \
((helpers != NULL) && ((helpers->free_table != 0) && (helpers->free_space != 0)) ? \
PyGSL_transform_helpers_free(helpers) : 0)
static void
PyGSL_transform_helpers_free(struct _pygsl_transform_help_rf_s * h)
{
FUNC_MESS_BEGIN();
assert(h->func);
DEBUG_MESS(3, "func @ %p", h->func);
if( (h->free_table == 1) && (h->table != NULL)){
assert(h->table);
DEBUG_MESS(3, "Free Table %p with func %p", h->table,
(void *) h->func->table_free);
h->func->table_free(h->table);
h->table = NULL;
h->free_table = 0;
}
if((h->free_space == 1) && (h->space != NULL)){
assert(h->space);
DEBUG_MESS(3, "Free Space %p with func %p", h->space,
(void *) h->func->space_free);
h->func->space_free(h->space);
h->space = NULL;
h->free_space = 0;
}
FUNC_MESS_END();
}
#ifdef _PyGSL_HAS_WAVELET
/*
* The two D function. Currently only supports wavelet2d_matrix...
*/
static PyObject *
PyGSL_transform_2d_(PyObject *self, PyObject *args, pygsl_transform_help_s *helps)
{
PyObject *ret = NULL, *data=NULL, *s_o = NULL;
PyArrayObject *m=NULL;
PyGSL_wavelet *wavelet = NULL;
gsl_matrix_view mv;
const enum radix_mode radix2 = helps->info->radix2;
const enum PyArray_TYPES input_array_type=helps->info->input_array_type;
const int sizeoftype = sizeof(double);
int call_n, line=-1;
switch(radix2){
case WAVELET:
if(!PyGSL_WAVELET_CHECK(self)){
gsl_error("Should be a wavelet method!", filename, line, GSL_ESANITY);
line = __LINE__ - 1;
goto fail;
}
wavelet = (PyGSL_wavelet *) self;
break;
default:
line = __LINE__;
gsl_error("Unknown switch!", filename, line, GSL_ESANITY);
goto fail;
}
/* Currently only implemented for wavelet 2d transform */
if(!PyArg_ParseTuple(args, "O|OO", &data, &s_o, &ret)){
line = __LINE__ - 1;
goto fail;
}
m = PyGSL_matrix_check(data, -1, -1,
PyGSL_BUILD_ARRAY_INFO(PyGSL_CONTIGUOUS | PyGSL_INPUT_ARRAY, input_array_type, 1, 2),
NULL, NULL, NULL);
if(m == NULL)
goto fail;
mv = gsl_matrix_view_array((double *)m->data, m->dimensions[0], m->dimensions[1]);
call_n = m->dimensions[0] + m->dimensions[1];
if (PyGSL_transform_helpers_alloc(s_o, NULL, helps->helpers, call_n) != GSL_SUCCESS){
line = __LINE__ -1;
goto fail;
}
if(PyGSL_ERROR_FLAG(helps->transform.wavelet2d(wavelet->wavelet, &mv.matrix, helps->helpers->space)) != GSL_SUCCESS){
line = __LINE__ -1;
goto fail;
}
PyGSL_TRANSFORM_HELPERS_FREE(helps->helpers);
ret = (PyObject *) m;
return ret;
fail:
FUNC_MESS("Fail");
PyGSL_add_traceback(module, filename, __FUNCTION__, line);
PyGSL_TRANSFORM_HELPERS_FREE(helps->helpers);
Py_XDECREF(m); m = NULL;
FUNC_MESS("Fail End");
return NULL;
}
#endif /* _PyGSL_HAS_WAVELET */
/*
* Catch all for all one dimensional functions.
*/
static PyObject *
PyGSL_transform_(PyObject *self, PyObject *args, pygsl_transform_help_s *helps)
{
PyObject *ret = NULL, *s_o = NULL, *t_o = NULL, *data = NULL;
PyArrayObject *result = NULL, *a=NULL, *r=NULL;
#ifdef _PyGSL_HAS_WAVELET
PyGSL_wavelet *wavelet = NULL;
#endif
double eps=1e-6;
int line = -1;
int length=-1;
/*
* how will it be called and what array length will I need.
* In the HalfComplexReal and the reversed case, an array of double
* data is provided. But the real array is of CDOUBLE. This takes the
* computation in place into account and minimizes the necessary
* copies.
*/
PyGSL_array_index_t n=0, call_n=0, return_n=0, strides=0;
const enum PyArray_TYPES input_array_type=helps->info->input_array_type,
output_array_type=helps->info->output_array_type;
const enum transform_mode mode = helps->info->mode;
const enum radix_mode radix2 = helps->info->radix2;
/*
* Helper for STRIDE recalculation. numpy/numarray use strides as
* multiples of char. But GSL uses the basis as common denominator.
* Only complex transforms use complex numbers, so only for them
* n_type will be two, and it will be one for all others.
*/
const enum pygsl_packed_type n_type = helps->info->packed;
/*
* Double mode or float mode?
* complex double or complex float ?
*/
const int datatype = helps->info->datatype;
/* and its size */
const int sizeoftype = (datatype == MODE_DOUBLE) ? sizeof(double) : sizeof(float);
/*
* Normaly the element 0 of the return array is passed. In exceptional
* cases, e.g. RealHalfComplex for RADIX_FREE an offset is needed!
*/
const int call_offset= helps->info->data_offset;
void * vdata;
FUNC_MESS_BEGIN();
assert(args);
/*
* I never know when I will call the _free function. So better to define
* them here already as any jump can go to fail and will start to free
* the tables.
*/
if(helps->helpers){
helps->helpers->free_table=0;
helps->helpers->free_space=0;
helps->helpers->table=NULL;
helps->helpers->space=NULL;
}
assert(mode>0);
assert(radix2>0);
assert(datatype>0);
/*
* Parse the arguments. Depends on the type how many arguments we are going
* to support.
*/
switch(radix2){
#ifdef _PyGSL_HAS_WAVELET
case WAVELET:
FUNC_MESS("Wavelet");
/*
* Yes thats a method here.... UiUi not sure if it is not too
* much in one function
*/
if(!PyGSL_WAVELET_CHECK(self)){
gsl_error("Should be a wavelet method!", filename, line, GSL_ESANITY);
line = __LINE__ - 1;
goto fail;
}
wavelet = (PyGSL_wavelet *) self;
if(!PyArg_ParseTuple(args, "O|OO", &data, &s_o, &ret)){
line = __LINE__ - 1;
goto fail;
}
break;
#endif
case RADIX_FREE:
FUNC_MESS("Radix_Free");
assert(helps->helpers);
switch(mode){
case ComplexComplex :
FUNC_MESS("ComplexComplex");
if(!PyArg_ParseTuple(args, "O|OOO", &data, &s_o, &t_o, &ret)){
line = __LINE__ - 1;
goto fail;
}
break;
case RealHalfcomplex:
FUNC_MESS("Real");
if(!PyArg_ParseTuple(args, "O|OO",&data, &s_o, &t_o)){
line = __LINE__ - 1;
goto fail;
}
break;
case HalfComplexReal:
FUNC_MESS("Halfcomplex");
if(!PyArg_ParseTuple(args, "O|iOOd",&data, &length, &s_o, &t_o, &eps)){
line = __LINE__ - 1;
goto fail;
}
break;
default:
line = __LINE__;
gsl_error("Unknown mode!", filename, line, GSL_ESANITY);
goto fail;
} /* mode */
break;
case RADIX_TWO:
FUNC_MESS("Radix_TWO");
assert(helps->helpers==NULL);
switch(mode){
case ComplexComplex :
FUNC_MESS("ComplexComplex");
if(!PyArg_ParseTuple(args, "O|O", &data, &ret)){
line = __LINE__ - 1;
goto fail;
}
break;
case RealHalfcomplex:
case HalfComplexReal:
FUNC_MESS("Real-Halfcomplex");
if(!PyArg_ParseTuple(args, "O",&data)){
line = __LINE__ - 1;
goto fail;
}
break;
default:
line = __LINE__;
gsl_error("Unknown mode!", filename, line, GSL_ESANITY);
goto fail;
} /* mode */
break;
default:
line = __LINE__;
DEBUG_MESS(3, "Radix2 was %d, wavelet = %d", (int) radix2, (int) WAVELET);
gsl_error("Unknown radix!", filename, line, GSL_ESANITY);
goto fail;
}/* radix2 */
FUNC_MESS("mode");
assert(input_array_type >0);
assert(output_array_type >0);
assert(n_type >0);
assert(sizeoftype >0);
DEBUG_MESS(3, "Double?Float %d IN Array Type = %d OUT Array Type = %d", datatype, input_array_type, output_array_type);
DEBUG_MESS(3, "n_type = %d sizeofbasis = %d", n_type, sizeoftype);
{
const int io_type = ((ret == NULL) ? PyGSL_INPUT_ARRAY : PyGSL_IO_ARRAY) | PyGSL_NON_CONTIGUOUS;
a = PyGSL_vector_check(data, -1, PyGSL_BUILD_ARRAY_INFO(io_type, input_array_type, 1, 1), NULL, NULL);
if(a == NULL){
line = __LINE__ - 1;
goto fail;
}
}
n = a->dimensions[0];
/*
* Calculate the size of the return array and the call array
* Generate the return array
*/
FUNC_MESS("Return Array");
switch(radix2){
#ifdef _PyGSL_HAS_WAVELET
case WAVELET:
return_n = n;
call_n = n;
if((r = PyGSL_Shadow_array((PyObject *) ret, (PyObject *) a, datatype)) == NULL){
line = __LINE__ -1;
goto fail;
}
break;
#endif
case RADIX_FREE:
FUNC_MESS("Radix Free");
switch(mode){
case ComplexComplex:
FUNC_MESS("ComplexComplex");
return_n = n;
call_n = n;
if((r = PyGSL_Shadow_array((PyObject *) ret, (PyObject *) a, datatype)) == NULL){
line = __LINE__ -1;
goto fail;
}
break;
case RealHalfcomplex:
FUNC_MESS("RealHalfcomplex");
return_n = n / 2 + 1;
call_n = n;
if(((r = (PyArrayObject *) PyGSL_New_Array(1, &return_n, output_array_type)) == NULL) ||
(PyGSL_copy_real_to_complex(r, a, datatype) != GSL_SUCCESS)){
line = __LINE__ -2;
goto fail;
}
break;
case HalfComplexReal:
FUNC_MESS("HalfcomplexReal");
call_n = PyGSL_guess_halfcomplex_length(a, length, eps, datatype);
return_n = call_n;
if(((r = (PyArrayObject *) PyGSL_New_Array(1, &return_n, output_array_type)) == NULL) ||
(PyGSL_ERROR_FLAG(PyGSL_copy_halfcomplex_to_real(r, a, eps, datatype)) != GSL_SUCCESS)){
line = __LINE__ - 2;
goto fail;
}
break;
default:
line = __LINE__ -1;
goto fail;
}/* mode */
break;
case RADIX_TWO:
FUNC_MESS("Radix TWO");
return_n = n;
call_n = n;
if((r = (PyArrayObject *) PyGSL_Copy_Array(a)) == NULL){
line = __LINE__ -1;
goto fail;
}
break;
default:
line = __LINE__ -1;
goto fail;
} /* radix2 */
Py_DECREF(a);
a = NULL;
/* make sure that the ntype was set */
assert(n_type > 0);
DEBUG_MESS(2, "Type(r) = %d, r->nd = %d, r->dimensions[0] = %d, Strides r->strides[0] %d",
r->descr->type_num, r->nd, r->dimensions[0], r->strides[0]);
if(PyGSL_STRIDE_RECALC(r->strides[0], n_type * sizeoftype, &strides) != GSL_SUCCESS){
line = __LINE__ -1;
goto fail;
}
DEBUG_MESS(2, "Strides r->strides[0] %ld strides = %ld", (long)r->strides[0], (long)strides);
/* build the helpers if necessary */
switch(radix2){
#ifdef _PyGSL_HAS_WAVELET
case WAVELET:
#endif
case RADIX_FREE:
FUNC_MESS("Build Helpers");
if (PyGSL_transform_helpers_alloc(s_o, t_o, helps->helpers, call_n) != GSL_SUCCESS){
line = __LINE__ -1;
goto fail;
}
break;
default:
/* No helpers needed! */ ;
}/* radix2 */
#if DEBUG > 0
if(PyGSL_Check_Array_Length(r, call_n, datatype, n_type) != GSL_SUCCESS){
line = __LINE__ -1;
goto fail;
}
#endif
vdata = (void *) r->data;
if(call_offset!=0){
switch(datatype){
case MODE_DOUBLE: vdata = (void *) (((double *) r->data) + call_offset); break;
case MODE_FLOAT: vdata = (void *) (((float *) r->data) + call_offset); break;
}
DEBUG_MESS(3, "Original Pointer at %p new pointer at %p", r->data, vdata);
}
/* call the function */
switch(radix2){
case RADIX_FREE:
FUNC_MESS("Transform free length");
assert(helps->helpers->table);
assert(helps->helpers->space);
DEBUG_MESS(3, "vdata = %f, strides = %ld, call_n = %ld",
*((double *)(vdata)), (long) strides, (long) call_n);
if(PyGSL_ERROR_FLAG(helps->transform.free(vdata, strides, call_n,
helps->helpers->table,
helps->helpers->space)) != GSL_SUCCESS){
line = __LINE__ -1;
goto fail;
}
DEBUG_MESS(3, "Transformed: r->data[0] = %f, strides = %ld, call_n = %ld",
*((double *)(r->data)), (long)strides, (long)call_n);
break;
case RADIX_TWO:
FUNC_MESS("Tranform radix2");
if(PyGSL_ERROR_FLAG(helps->transform.radix2(vdata, strides, call_n)) != GSL_SUCCESS){
line = __LINE__ -1;
goto fail;
}
break;
#ifdef _PyGSL_HAS_WAVELET
case WAVELET:
FUNC_MESS("Tranform wavelet");
assert(wavelet);
assert(helps->helpers->space);
if(PyGSL_ERROR_FLAG(helps->transform.wavelet(wavelet->wavelet, (double *)vdata, strides, call_n,
helps->helpers->space)) != GSL_SUCCESS){
line = __LINE__ -1;
goto fail;
}
break;
#endif
default:
line = __LINE__ -1;
goto fail;
}
/* restore the array */
switch(radix2){
case RADIX_FREE:
switch(mode){
case RealHalfcomplex:
FUNC_MESS("halfcomplex unpack");
if(PyGSL_fft_halfcomplex_unpack(r, call_n, datatype) != GSL_SUCCESS){
line = __LINE__ -1;
goto fail;
}
break;
case RealReal:
case HalfComplexReal:
case ComplexComplex: break;
}
break;
default:
/* No helpers needed! */ ;
}
PyGSL_TRANSFORM_HELPERS_FREE(helps->helpers);
result = r;
FUNC_MESS_END();
return (PyObject *) result;
fail:
FUNC_MESS("Fail");
PyGSL_add_traceback(module, filename, __FUNCTION__, line);
PyGSL_TRANSFORM_HELPERS_FREE(helps->helpers);
Py_XDECREF(r); r = NULL;
FUNC_MESS("Fail End");
return NULL;
}
/*
* Local Variables:
* mode: C
* c-file-style: "python"
* End:
*/
-------------------------------------------------------------------------
Take Surveys. Earn Cash. Influence the Future of IT
Join SourceForge.net's Techsay panel and you'll get the chance to share your
opinions on IT & business topics through brief surveys -- and earn cash
http://www.techsay.com/default.php?page=join.php&p=sourceforge&CID=DEVDEV
_______________________________________________
pygsl-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/pygsl-discuss
