diff --git a/.gitignore b/.gitignore index 779abc6..d739953 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,8 @@ *.o *.so -*.dll \ No newline at end of file +*.dll +*.exe +*.txt +*.dat +*.bin +*.csv \ No newline at end of file diff --git a/Makefile.dirs b/Makefile.dirs index e998cf9..850d1d8 100644 --- a/Makefile.dirs +++ b/Makefile.dirs @@ -1,38 +1,18 @@ CC = gcc INC_DIR = include +RELEASE_DIR = release/lib ifeq ($(OS),Windows_NT) DSPL_LIBNAME = libdspl.dll DEF_OS = WIN_OS LFLAGS = -lm - ifeq ($(PROCESSOR_ARCHITEW6432),AMD64) - RELEASE_DIR = release/lib/win64 - EXTLIB_DIR = ext/win64 - else - ifeq ($(PROCESSOR_ARCHITECTURE),AMD64) - RELEASE_DIR = release/lib/win64 - EXTLIB_DIR = ext/win64 - endif - ifeq ($(PROCESSOR_ARCHITECTURE),x86) - RELEASE_DIR = release/lib/win32 - EXTLIB_DIR = ext/win32 - endif - endif else UNAME_S := $(shell uname -s) UNAME_P := $(shell uname -p) ifeq ($(UNAME_S),Linux) DSPL_LIBNAME = libdspl.so DEF_OS = LINUX_OS - LFLAGS = -lm -ldl - ifeq ($(UNAME_P),x86_64) - RELEASE_DIR = release/lib/linux64 - EXTLIB_DIR = ext/linux64 - endif - ifneq ($(filter %86,$(UNAME_P)),) - RELEASE_DIR = release/lib/linux32 - EXTLIB_DIR = ext/linux32 - endif + LFLAGS = -lm -ldl endif endif diff --git a/Makefile.dspl b/Makefile.dspl index 11fae30..9e7dffe 100644 --- a/Makefile.dspl +++ b/Makefile.dspl @@ -21,10 +21,10 @@ all: $(RELEASE_DIR)/$(LIB_NAME) $(RELEASE_DIR)/$(LIB_NAME): $(OBJ_FILES) - $(CC) -shared -o $(RELEASE_DIR)/$(LIB_NAME) $(OBJ_FILES) -L$(EXTLIB_DIR) -lfftw3 -lblas -lm + $(CC) -shared -o $(RELEASE_DIR)/$(LIB_NAME) $(OBJ_FILES) -lm $(OBJ_DIR)/%.o:$(SRC_DIR)/%.c - $(CC) $(CFLAGS) $< -o $@ -L$(EXTLIB_DIR) -lfftw3 -lblas -llapack -lgfortran -lm + $(CC) $(CFLAGS) $< -o $@ -lm clean: rm -f $(OBJ_DIR)/*.o diff --git a/dspl/src/array.c b/dspl/src/array.c deleted file mode 100644 index 55a8ca6..0000000 --- a/dspl/src/array.c +++ /dev/null @@ -1,100 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - - -#include -#include -#include -#include "dspl.h" - -/************************************************************************************************** -Concntenate arrays -***************************************************************************************************/ -int concat(void* a, size_t na, void *b, size_t nb, void* c) -{ - if(!a || !b || !c || c == b) - return ERROR_PTR; - if(na < 1 || nb < 1) - return ERROR_SIZE; - - if(c != a) - memcpy(c, a, na); - - memcpy(c+na, b, nb); - return RES_OK; -} - - - - - - -/************************************************************************************************** -Flip real array in place -***************************************************************************************************/ -int flipip(double* x, int n) -{ - int k; - double tmp; - if(!x) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - - for(k = 0; k < n/2; k++) - { - tmp = x[k]; - x[k] = x[n-1-k]; - x[n-1-k] = tmp; - } - return RES_OK; -} - - - -/************************************************************************************************** -Flip complex array in place -***************************************************************************************************/ -int flipip_cmplx(complex_t* x, int n) -{ - int k; - complex_t tmp; - if(!x) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - - for(k = 0; k < n/2; k++) - { - RE(tmp) = RE(x[k]); - RE(x[k]) = RE(x[n-1-k]); - RE(x[n-1-k]) = RE(tmp); - - IM(tmp) = IM(x[k]); - IM(x[k]) = IM(x[n-1-k]); - IM(x[n-1-k]) = IM(tmp); - } - return RES_OK; -} - - - - - diff --git a/dspl/src/blas_lev1.c b/dspl/src/blas_lev1.c deleted file mode 100644 index 161ec5f..0000000 --- a/dspl/src/blas_lev1.c +++ /dev/null @@ -1,58 +0,0 @@ -/* -* Copyright (c) 2015-2018 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of libdspl-2.0. -* -* is free software: you can redistribute it and/or modify -* it under the terms of the GNU Lesser General Public License as published by -* the Free Software Foundation, either version 3 of the License, or -* (at your option) any later version. -* -* DSPL 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 Lesser General Public License -* along with Foobar. If not, see . -*/ - - -#include -#include -#include -#include "dspl.h" - - - -void dscal_(int*, double*, double*, int*); - - - - -/************************************************************************************************** -DSCAL scales a vector by a constant -***************************************************************************************************/ -int DSPL_API blas_dscal(int n, double a, double* x, int incx) -{ - if(!x) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - if(incx < 1) - return ERROR_INC_SIZE; - - dscal_(&n, &a, x, &incx); - - return RES_OK; -} - - - - - - - - - diff --git a/dspl/src/cheby.c b/dspl/src/cheby.c index b1105b9..7693968 100644 --- a/dspl/src/cheby.c +++ b/dspl/src/cheby.c @@ -24,7 +24,6 @@ #include "dspl.h" - /************************************************************************************************** Chebyshev polynomials of the first kind ***************************************************************************************************/ @@ -119,7 +118,4 @@ int DSPL_API cheby_poly2(double* x, int n, int ord, double* y) } } return RES_OK; - - - } diff --git a/dspl/src/complex.c b/dspl/src/complex.c deleted file mode 100644 index d63fff9..0000000 --- a/dspl/src/complex.c +++ /dev/null @@ -1,236 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - - -#include -#include -#include "dspl.h" - - - -/************************************************************************************************** -Acos complex -***************************************************************************************************/ -int acos_cmplx(complex_t* x, int n, complex_t *y) -{ - int k, res; - double pi2 = 0.5 * M_PI; - - res = asin_cmplx(x, n, y); - if(res != RES_OK) - return res; - - for(k = 0; k < n; k++) - { - RE(y[k]) = pi2 - RE(y[k]); - IM(y[k]) = - IM(y[k]); - } - return RES_OK; -} - - - - -/************************************************************************************************** -Asin complex -***************************************************************************************************/ -int asin_cmplx(complex_t* x, int n, complex_t *y) -{ - int k; - complex_t tmp; - if(!x || !y) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - - for(k = 0; k < n; k++) - { - RE(tmp) = 1.0 - CMRE(x[k], x[k]); // 1-x[k]^2 - IM(tmp) = - CMIM(x[k], x[k]); // 1-x[k]^2 - sqrt_cmplx(&tmp, 1, y+k); // sqrt(1 - x[k]^2) - RE(y[k]) -= IM(x[k]); // j * x[k] + sqrt(1 - x[k]^2) - IM(y[k]) += RE(x[k]); // j * x[k] + sqrt(1 - x[k]^2) - log_cmplx(y+k, 1, &tmp); // log( j * x[k] + sqrt(1 - x[k]^2) ) - RE(y[k]) = IM(tmp); // -j * log( j * x[k] + sqrt(1 - x[k]^2) ) - IM(y[k]) = -RE(tmp); // -j * log( j * x[k] + sqrt(1 - x[k]^2) ) - } - return RES_OK; -} - - - -/************************************************************************************************** -convert double array to a complex array -***************************************************************************************************/ -int re2cmplx(double* x, int n, complex_t *y) -{ - int k; - if(!x || !y) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - - for(k = 0; k < n; k++) - { - RE(y[k]) = x[k]; - IM(y[k]) = 0.0; - } - return RES_OK; -} - - - - - -/************************************************************************************************** -convert complex array to a re and im arrays -***************************************************************************************************/ -int cmplx2re(complex_t* x, int n, double *re, double *im) -{ - int k; - if(!x) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - - if(re) - { - for(k = 0; k < n; k++) - re[k] = RE(x[k]); - } - if(im) - { - for(k = 0; k < n; k++) - im[k] = IM(x[k]); - } - return RES_OK; -} - - - - - - -/************************************************************************************************** -Complex cosine -***************************************************************************************************/ -int cos_cmplx(complex_t* x, int n, complex_t *y) -{ - int k; - double ep, em, sx, cx; - if(!x || !y) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - - for(k = 0; k < n; k++) - { - ep = exp( IM(x[k])); - em = exp(-IM(x[k])); - sx = 0.5 * sin(RE(x[k])); - cx = 0.5 * cos(RE(x[k])); - RE(y[k]) = cx * (em + ep); - IM(y[k]) = sx * (em - ep); - } - return RES_OK; -} - - - - -/************************************************************************************************** -Complex cosine -***************************************************************************************************/ -int sin_cmplx(complex_t* x, int n, complex_t *y) -{ - int k; - double ep, em, sx, cx; - if(!x || !y) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - - for(k = 0; k < n; k++) - { - ep = exp( IM(x[k])); - em = exp(-IM(x[k])); - sx = 0.5 * sin(RE(x[k])); - cx = 0.5 * cos(RE(x[k])); - RE(y[k]) = sx * (em + ep); - IM(y[k]) = cx * (ep - em); - } - return RES_OK; -} - - - - -/************************************************************************************************** -Logarithm complex -***************************************************************************************************/ -int log_cmplx(complex_t* x, int n, complex_t *y) -{ - int k; - if(!x || !y) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - - for(k = 0; k < n; k++) - { - RE(y[k]) = 0.5 * log(ABSSQR(x[k])); - IM(y[k]) = atan2(IM(x[k]), RE(x[k])); - } - return RES_OK; -} - - - - -/************************************************************************************************** -SQRT complex -***************************************************************************************************/ -int sqrt_cmplx(complex_t* x, int n, complex_t *y) -{ - int k; - double r, zr; - complex_t t; - if(!x || !y) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - - for(k = 0; k < n; k++) - { - r = ABS(x[k]); - RE(t) = RE(x[k]) + r; - IM(t) = IM(x[k]); - zr = 1.0 / ABS(t); - r = sqrt(r); - RE(y[k]) = RE(t) * zr * r; - IM(y[k]) = IM(t) * zr * r; - - } - return RES_OK; -} - - - - diff --git a/dspl/src/composing.c b/dspl/src/composing.c deleted file mode 100644 index 50bc14a..0000000 --- a/dspl/src/composing.c +++ /dev/null @@ -1,122 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - - -#include -#include -#include "dspl.h" - - -/************************************************************************************************** -Composition of rational functions -***************************************************************************************************/ -int rcompos(double *b, double *a, int n, double *c, double *d, int p, double *beta, double *alpha) -{ - int k2, i, k, pn, pd, ln, ld, k2s, nk2s; - double *num = NULL, *den = NULL, *ndn = NULL, *ndd = NULL; - int res; - - if (!a || !b || !c || !d || !beta || !alpha) - return ERROR_PTR; - - if(n < 0 || p < 0) - return ERROR_SIZE; - - - k2 = (n*p)+1; - k2s = k2*sizeof(double); /* alpha and beta size */ - nk2s = (n+1)*k2*sizeof(double); /* num, den, ndn and ndd size */ - - num = (double*)malloc(nk2s); - den = (double*)malloc(nk2s); - ndn = (double*)malloc(nk2s); - ndd = (double*)malloc(nk2s); - - memset(num, 0, nk2s); - memset(den, 0, nk2s); - memset(ndn, 0, nk2s); - memset(ndd, 0, nk2s); - - - num[0] = den[0] = 1.0; - pn = 0; - ln = 1; - for(i = 1; i < n+1; i++) - { - res = conv(num+pn, ln, c, p+1, num+pn+k2); - if(res!=RES_OK) - goto exit_label; - res = conv(den+pn, ln, d, p+1, den+pn+k2); - if(res!=RES_OK) - goto exit_label; - pn += k2; - ln += p; - } - - pn = 0; - pd = n*k2; - ln = 1; - ld = k2; - - for (i = 0; i < n+1; i++) - { - res = conv(num + pn, ln, den + pd, ld, ndn + i*k2); - if(res!=RES_OK) - goto exit_label; - ln += p; - ld -= p; - pn += k2; - pd -= k2; - } - - for (i = 0; i < n+1; i++) - { - for (k = 0; k < k2; k++) - { - ndd[i*k2 + k] = ndn[i*k2 + k] * a[i]; - ndn[i*k2 + k] *= b[i]; - } - } - - memset(alpha, 0, k2s); - memset(beta, 0, k2s); - - for (k = 0; k < k2; k++) - { - for (i = 0; i < n+1; i++) - { - beta[k] += ndn[i*k2 + k]; - alpha[k] += ndd[i*k2 + k]; - } - } - - res = RES_OK; -exit_label: - if(num) - free(num); - if(den) - free(den); - if(ndn) - free(ndn); - if(ndd) - free(ndd); - - return res; -} diff --git a/dspl/src/conv.c b/dspl/src/conv.c index 0b999f3..3063b47 100644 --- a/dspl/src/conv.c +++ b/dspl/src/conv.c @@ -1,11 +1,11 @@ /* -* Copyright (c) 2015-2017 Sergey Bakhurin +* Copyright (c) 2015-2018 Sergey Bakhurin * Digital Signal Processing Library [http://dsplib.org] * -* This file is part of DSPL. +* This file is part of libdspl-2.0. * * is free software: you can redistribute it and/or modify -* it under the terms of the GNU General Public License as published by +* it under the terms of the GNU Lesser General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * @@ -14,7 +14,7 @@ * 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 +* You should have received a copy of the GNU Lesser General Public License * along with Foobar. If not, see . */ @@ -29,7 +29,7 @@ /************************************************************************************************** Real vectors linear convolution ***************************************************************************************************/ -int conv(double* a, int na, double* b, int nb, double* c) +int DSPL_API conv(double* a, int na, double* b, int nb, double* c) { int k; int n; @@ -72,7 +72,7 @@ int conv(double* a, int na, double* b, int nb, double* c) /************************************************************************************************** Complex vectors linear convolution ***************************************************************************************************/ -int conv_cmplx(complex_t* a, int na, complex_t* b, int nb, complex_t* c) +int DSPL_API conv_cmplx(complex_t* a, int na, complex_t* b, int nb, complex_t* c) { int k; int n; @@ -119,7 +119,7 @@ int conv_cmplx(complex_t* a, int na, complex_t* b, int nb, complex_t* c) /************************************************************************************************** IIR FILTER for real vector **************************************************************************************************/ -int filter_iir(double* b, double* a, int ord, double* x, int n, double* y) +int DSPL_API filter_iir(double* b, double* a, int ord, double* x, int n, double* y) { double* buf = NULL; double* an = NULL; diff --git a/dspl/src/ellipj.c b/dspl/src/ellipj.c deleted file mode 100644 index 2c74c81..0000000 --- a/dspl/src/ellipj.c +++ /dev/null @@ -1,494 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ -#include -#include -#include -#include "dspl.h" - - - -/************************************************************************************************** -inverse cd function -***************************************************************************************************/ -int ellip_acd(double* w, int n, double k, double* u) -{ - double lnd[ELLIP_ITER], t; - int i, m; - - if(!u || !w) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - if(k < 0.0 || k>= 1.0) - return ERROR_ELLIP_MODULE; - - ellip_landen(k,ELLIP_ITER, lnd); - - - for(m = 0; m < n; m++) - { - u[m] = w[m]; - for(i = 1; i < ELLIP_ITER; i++) - { - t = lnd[i-1]*u[m]; - t *= t; - t = 1.0 + sqrt(1.0 - t); - u[m] = 2.0 * u[m] / (t+t*lnd[i]); - } - u[m] = 2.0 * acos(u[m]) / M_PI; - } - return RES_OK; -} - - - - - -/************************************************************************************************** -inverse cd function -***************************************************************************************************/ -int ellip_acd_cmplx(complex_t* w, int n, double k, complex_t* u) -{ - double lnd[ELLIP_ITER], t; - complex_t tmp0, tmp1; - int i, m; - - if(!u || !w) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - if(k < 0.0 || k>= 1.0) - return ERROR_ELLIP_MODULE; - - ellip_landen(k,ELLIP_ITER, lnd); - - - for(m = 0; m < n; m++) - { - RE(u[m]) = RE(w[m]); - IM(u[m]) = IM(w[m]); - for(i = 1; i < ELLIP_ITER; i++) - { - RE(tmp0) = lnd[i-1]*RE(u[m]); - IM(tmp0) = lnd[i-1]*IM(u[m]); - RE(tmp1) = 1.0 - CMRE(tmp0, tmp0); - IM(tmp1) = - CMIM(tmp0, tmp0); - - sqrt_cmplx(&tmp1, 1, &tmp0); - RE(tmp0) += 1.0; - - RE(tmp1) = RE(tmp0) * (1.0 + lnd[i]); - IM(tmp1) = IM(tmp0) * (1.0 + lnd[i]); - - t = 2.0 / ABSSQR(tmp1); - - RE(tmp0) = t * CMCONJRE(tmp1, u[m]); - IM(tmp0) = t * CMCONJIM(tmp1, u[m]); - - RE(u[m]) = RE(tmp0); - IM(u[m]) = IM(tmp0); - - } - acos_cmplx(&tmp0, 1, u+m); - t = 2.0 / M_PI; - RE(u[m]) *= t; - IM(u[m]) *= t; - } - return RES_OK; -} - - - - - -/************************************************************************************************** -inverse sn function -***************************************************************************************************/ -int ellip_asn(double* w, int n, double k, double* u) -{ - double lnd[ELLIP_ITER], t; - int i, m; - - if(!u || !w) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - if(k < 0.0 || k>= 1.0) - return ERROR_ELLIP_MODULE; - - ellip_landen(k,ELLIP_ITER, lnd); - - - for(m = 0; m < n; m++) - { - u[m] = w[m]; - for(i = 1; i < ELLIP_ITER; i++) - { - t = lnd[i-1]*u[m]; - t *= t; - t = 1.0 + sqrt(1.0 - t); - u[m] = 2.0 * u[m] / (t+t*lnd[i]); - } - u[m] = 2.0 * asin(u[m]) / M_PI; - } - return RES_OK; -} - - - - - -/************************************************************************************************** -inverse sn function -***************************************************************************************************/ -int ellip_asn_cmplx(complex_t* w, int n, double k, complex_t* u) -{ - double lnd[ELLIP_ITER], t; - complex_t tmp0, tmp1; - int i, m; - - if(!u || !w) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - if(k < 0.0 || k>= 1.0) - return ERROR_ELLIP_MODULE; - - ellip_landen(k,ELLIP_ITER, lnd); - - - for(m = 0; m < n; m++) - { - RE(u[m]) = RE(w[m]); - IM(u[m]) = IM(w[m]); - for(i = 1; i < ELLIP_ITER; i++) - { - RE(tmp0) = lnd[i-1]*RE(u[m]); - IM(tmp0) = lnd[i-1]*IM(u[m]); - RE(tmp1) = 1.0 - CMRE(tmp0, tmp0); - IM(tmp1) = - CMIM(tmp0, tmp0); - - sqrt_cmplx(&tmp1, 1, &tmp0); - RE(tmp0) += 1.0; - - RE(tmp1) = RE(tmp0) * (1.0 + lnd[i]); - IM(tmp1) = IM(tmp0) * (1.0 + lnd[i]); - - t = 2.0 / ABSSQR(tmp1); - - RE(tmp0) = t * CMCONJRE(tmp1, u[m]); - IM(tmp0) = t * CMCONJIM(tmp1, u[m]); - - RE(u[m]) = RE(tmp0); - IM(u[m]) = IM(tmp0); - - } - asin_cmplx(&tmp0, 1, u+m); - t = 2.0 / M_PI; - RE(u[m]) *= t; - IM(u[m]) *= t; - } - return RES_OK; -} - - - - - - -/************************************************************************************************** -Elliptic cd function -***************************************************************************************************/ -int ellip_cd(double* u, int n, double k, double* y) -{ - double lnd[ELLIP_ITER]; - int i, m; - - if(!u || !y) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - if(k < 0.0 || k>= 1.0) - return ERROR_ELLIP_MODULE; - - ellip_landen(k,ELLIP_ITER, lnd); - - - for(m = 0; m < n; m++) - { - y[m] = cos(u[m] * M_PI * 0.5); - for(i = ELLIP_ITER-1; i>0; i--) - { - y[m] = (1.0 + lnd[i]) / (1.0 / y[m] + lnd[i]*y[m]); - } - } - return RES_OK; -} - - - - - - -/************************************************************************************************** -Elliptic cd function -***************************************************************************************************/ -int ellip_cd_cmplx(complex_t* u, int n, double k, complex_t* y) -{ - double lnd[ELLIP_ITER], t; - int i, m; - complex_t tmp; - - if(!u || !y) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - if(k < 0.0 || k>= 1.0) - return ERROR_ELLIP_MODULE; - - ellip_landen(k,ELLIP_ITER, lnd); - - - for(m = 0; m < n; m++) - { - RE(tmp) = RE(u[m]) * M_PI * 0.5; - IM(tmp) = IM(u[m]) * M_PI * 0.5; - - cos_cmplx(&tmp, 1, y+m); - - for(i = ELLIP_ITER-1; i>0; i--) - { - t = 1.0 / ABSSQR(y[m]); - - RE(tmp) = RE(y[m]) * t + RE(y[m]) * lnd[i]; - IM(tmp) = -IM(y[m]) * t + IM(y[m]) * lnd[i]; - - t = (1.0 + lnd[i]) / ABSSQR(tmp); - - RE(y[m]) = RE(tmp) * t; - IM(y[m]) = -IM(tmp) * t; - - } - } - return RES_OK; -} - - - - - - -/************************************************************************************************** -Landen transform -***************************************************************************************************/ -int ellip_landen(double k, int n, double* y) -{ - int i; - y[0] = k; - - if(!y) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - if(k < 0.0 || k>= 1.0) - return ERROR_ELLIP_MODULE; - - for(i = 1; i < n; i++) - { - y[i] = y[i-1] / (1.0 + sqrt(1.0 - y[i-1] * y[i-1])); - y[i] *= y[i]; - } - - return RES_OK; -} - - - - - -/************************************************************************************************** -Elliptic modular equation -***************************************************************************************************/ -int ellip_modulareq(double rp, double rs, int ord, double *k) -{ - double ep, es, ke, kp, t, sn; - int i, L, r; - - if(rp < 0 || rp == 0) - return ERROR_FILTER_RP; - if(rs < 0 || rs == 0) - return ERROR_FILTER_RS; - if(ord < 1) - return ERROR_FILTER_ORD; - if(!k) - return ERROR_PTR; - - - ep = sqrt(pow(10.0, rp*0.1)-1.0); - es = sqrt(pow(10.0, rs*0.1)-1.0); - - ke = ep/es; - - ke = sqrt(1.0 - ke*ke); - - r = ord % 2; - L = (ord-r)/2; - - kp = 1.0; - for(i = 0; i < L; i++) - { - t = (double)(2*i+1) / (double)ord; - ellip_sn(&t, 1, ke, &sn); - sn*=sn; - kp *= sn*sn; - } - - kp *= pow(ke, (double)ord); - *k = sqrt(1.0 - kp*kp); - - return RES_OK; - -} - - - -/************************************************************************************************** -Elliptic rational function -***************************************************************************************************/ -int ellip_rat(double* w, int n, int ord, double k, double* u) -{ - double t, xi, w2, xi2, k2; - int i, m, r, L; - - if(!u || !w) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - if(k < 0.0 || k>= 1.0) - return ERROR_ELLIP_MODULE; - - r = ord%2; - L = (ord-r)/2; - - if(r) - memcpy(u, w, n*sizeof(double)); - else - { - for(m = 0; m < n; m++) - { - u[m] = 1.0; - } - } - - k2 = k*k; - for(i = 0; i < L; i++) - { - t = (double)(2*i+1) / (double)ord; - ellip_cd(&t, 1, k, &xi); - xi2 = xi*xi; - for(m = 0; m < n; m++) - { - w2 = w[m]*w[m]; - u[m] *= (w2 - xi2) / (1.0 - w2 * k2 * xi2); - u[m] *= (1.0 - k2*xi2) / (1.0 - xi2); - } - } - return RES_OK; -} - - - - - - -/************************************************************************************************** -Elliptic sn function -***************************************************************************************************/ -int ellip_sn(double* u, int n, double k, double* y) -{ - double lnd[ELLIP_ITER]; - int i, m; - - if(!u || !y) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - if(k < 0.0 || k>= 1.0) - return ERROR_ELLIP_MODULE; - - ellip_landen(k,ELLIP_ITER, lnd); - - - for(m = 0; m < n; m++) - { - y[m] = sin(u[m] * M_PI * 0.5); - for(i = ELLIP_ITER-1; i>0; i--) - { - y[m] = (1.0 + lnd[i]) / (1.0 / y[m] + lnd[i]*y[m]); - } - } - return RES_OK; -} - - -/************************************************************************************************** -Elliptic sn function -***************************************************************************************************/ -int ellip_sn_cmplx(complex_t* u, int n, double k, complex_t* y) -{ - double lnd[ELLIP_ITER], t; - int i, m; - complex_t tmp; - - if(!u || !y) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - if(k < 0.0 || k>= 1.0) - return ERROR_ELLIP_MODULE; - - ellip_landen(k,ELLIP_ITER, lnd); - - - for(m = 0; m < n; m++) - { - RE(tmp) = RE(u[m]) * M_PI * 0.5; - IM(tmp) = IM(u[m]) * M_PI * 0.5; - - sin_cmplx(&tmp, 1, y+m); - - for(i = ELLIP_ITER-1; i>0; i--) - { - t = 1.0 / ABSSQR(y[m]); - - RE(tmp) = RE(y[m]) * t + RE(y[m]) * lnd[i]; - IM(tmp) = -IM(y[m]) * t + IM(y[m]) * lnd[i]; - - t = (1.0 + lnd[i]) / ABSSQR(tmp); - - RE(y[m]) = RE(tmp) * t; - IM(y[m]) = -IM(tmp) * t; - - } - } - return RES_OK; -} - - diff --git a/dspl/src/fft.c b/dspl/src/fft.c deleted file mode 100644 index 0b4adc7..0000000 --- a/dspl/src/fft.c +++ /dev/null @@ -1,167 +0,0 @@ -/* -* Copyright (c) 2015-2018 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of libdspl-2.0. -* -* is free software: you can redistribute it and/or modify -* it under the terms of the GNU Lesser General Public License as published by -* the Free Software Foundation, either version 3 of the License, or -* (at your option) any later version. -* -* DSPL 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 Lesser General Public License -* along with Foobar. If not, see . -*/ - -#include -#include -#include -#include "dspl.h" -#include "fftw3.h" - - - -/************************************************************************************************** -Create FFT object -***************************************************************************************************/ -int DSPL_API fft_create( fft_t *pfft, int n) -{ - if(!pfft) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - - printf("1\n"); - if(pfft->in) - { - if(pfft->size != n) - pfft->in = (complex_t*) realloc(pfft->in, n*sizeof(complex_t)); - } - else - pfft->in = (complex_t*) malloc(n*sizeof(complex_t)); - - if(pfft->out) - { - if(pfft->size != n) - pfft->out = (complex_t*) realloc(pfft->out, n*sizeof(complex_t)); - } - else - pfft->out = (complex_t*) malloc(n*sizeof(complex_t)); - - pfft->size = n; - printf("2\n"); - - if(pfft->pfftw) - fftw_destroy_plan(pfft->pfftw); - - printf("3\n"); - pfft->pfftw = (void*)fftw_plan_dft_1d(n, pfft->in, pfft->out, FFTW_FORWARD, FFTW_ESTIMATE); - printf("4\n"); - - if(!pfft->pfftw) - return ERROR_FFT_CREATE; - - return RES_OK; -} - - - -/************************************************************************************************** -Destroy FFT object -***************************************************************************************************/ -void DSPL_API fft_destroy(fft_t *pfft) -{ - if(pfft) - return; - - if(pfft->pfftw) - fftw_destroy_plan(pfft->pfftw); - - if(pfft->in) - free(pfft->in); - - if(pfft->out) - free(pfft->out); - - pfft->size = 0; -} - - - -/************************************************************************************************** -Destroy FFT object -***************************************************************************************************/ -int DSPL_API fft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y) -{ - if(!x || !y || !pfft) - return ERROR_PTR; - - if(n<1) - return ERROR_SIZE; - - if(n != pfft->size) - fft_create(pfft, n); - - memcpy(pfft->in, x, n*sizeof(complex_t)); - - fftw_execute(pfft->pfftw); - - memcpy(y, pfft->out, n*sizeof(complex_t)); - return RES_OK; - -} - - - - -/************************************************************************************************** -FFT shifting -***************************************************************************************************/ -int DSPL_API fft_shift(double* x, int n, double* y) -{ - int n2, r; - int k; - double tmp; - double *buf; - - if(!x || !y) - return ERROR_PTR; - - if(n<1) - return ERROR_SIZE; - - r = n%2; - if(!r) - { - n2 = n>>1; - for(k = 0; k < n2; k++) - { - tmp = x[k]; - y[k] = x[k+n2]; - y[k+n2] = tmp; - } - } - else - { - n2 = (n-1) >> 1; - buf = (double*) malloc(n2*sizeof(double)); - memcpy(buf, x, n2*sizeof(double)); - memcpy(y, x+n2, (n2+1)*sizeof(double)); - memcpy(y+n2+1, buf, n2*sizeof(double)); - free(buf); - } - return RES_OK; -} - - - - - - - - diff --git a/dspl/src/fillarray.c b/dspl/src/fillarray.c index 655cb9d..70d7195 100644 --- a/dspl/src/fillarray.c +++ b/dspl/src/fillarray.c @@ -1,11 +1,11 @@ /* -* Copyright (c) 2015-2017 Sergey Bakhurin +* Copyright (c) 2015-2018 Sergey Bakhurin * Digital Signal Processing Library [http://dsplib.org] * -* This file is part of DSPL. +* This file is part of libdspl-2.0. * * is free software: you can redistribute it and/or modify -* it under the terms of the GNU General Public License as published by +* it under the terms of the GNU Lesser General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * @@ -14,7 +14,7 @@ * 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 +* You should have received a copy of the GNU Lesser General Public License * along with Foobar. If not, see . */ @@ -28,7 +28,7 @@ /************************************************************************************************** Linspace array filling ***************************************************************************************************/ -int linspace(double x0, double x1, int n, int type, double* x) +int DSPL_API linspace(double x0, double x1, int n, int type, double* x) { double dx; int k; @@ -65,7 +65,7 @@ int linspace(double x0, double x1, int n, int type, double* x) /************************************************************************************************** Logspace array filling ***************************************************************************************************/ -int logspace(double x0, double x1, int n, int type, double* x) +int DSPL_API logspace(double x0, double x1, int n, int type, double* x) { double mx, a, b; int k; diff --git a/dspl/src/filter_an.c b/dspl/src/filter_an.c deleted file mode 100644 index e244b51..0000000 --- a/dspl/src/filter_an.c +++ /dev/null @@ -1,171 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - -#include -#include -#include -#include "dspl.h" - - - -/************************************************************************************************** -Complex frequency response of an analog filter H(s) -***************************************************************************************************/ -int freqs(double* b, double* a, int ord, double* w, int n, complex_t *h) -{ - - complex_t jw; - complex_t *bc = NULL; - complex_t *ac = NULL; - complex_t num, den; - double mag; - int k; - int res; - - if(!b || !a || !w || !h) - return ERROR_PTR; - if(ord<0) - return ERROR_FILTER_ORD; - if(n<1) - return ERROR_SIZE; - - RE(jw) = 0.0; - - bc = (complex_t*) malloc((ord+1) * sizeof(complex_t)); - res = re2cmplx(b, ord+1, bc); - if( res!=RES_OK ) - goto exit_label; - - ac = (complex_t*) malloc((ord+1) * sizeof(complex_t)); - res = re2cmplx(a, ord+1, ac); - if( res!=RES_OK ) - goto exit_label; - - for(k = 0; k < n; k++) - { - IM(jw) = w[k]; - res = polyval_cmplx(bc, ord, &jw, 1, &num); - if(res != RES_OK) - goto exit_label; - res = polyval_cmplx(ac, ord, &jw, 1, &den); - if(res != RES_OK) - goto exit_label; - mag = ABSSQR(den); - if(mag == 0.0) - { - res = ERROR_DIV_ZERO; - goto exit_label; - } - mag = 1.0 / mag; - RE(h[k]) = CMCONJRE(num, den) * mag; - IM(h[k]) = CMCONJIM(num, den) * mag; - - } - res = RES_OK; -exit_label: - if(bc) - free(bc); - if(ac) - free(ac); - return res; -} - - - - - - - -/************************************************************************************************** -Complex frequency response of a digital filter H(z) -**************************************************************************************************/ -int freqz(double* b, double* a, int ord, double* w, int n, complex_t *h) -{ - - complex_t jw; - complex_t *bc = NULL; - complex_t *ac = NULL; - complex_t num, den; - double mag; - int k; - int res; - - if(!b || !w || !h) - return ERROR_PTR; - if(ord<0) - return ERROR_FILTER_ORD; - if(n<1) - return ERROR_SIZE; - - - bc = (complex_t*) malloc((ord+1) * sizeof(complex_t)); - res = re2cmplx(b, ord+1, bc); - if( res!=RES_OK ) - goto exit_label; - - if(a) - { - // IIR filter if a != NULL - ac = (complex_t*) malloc((ord+1) * sizeof(complex_t)); - res = re2cmplx(a, ord+1, ac); - if( res!=RES_OK ) - goto exit_label; - for(k = 0; k < n; k++) - { - RE(jw) = cos(w[k]); - IM(jw) = -sin(w[k]); - res = polyval_cmplx(bc, ord, &jw, 1, &num); - if(res != RES_OK) - goto exit_label; - res = polyval_cmplx(ac, ord, &jw, 1, &den); - if(res != RES_OK) - goto exit_label; - mag = ABSSQR(den); - if(mag == 0.0) - { - res = ERROR_DIV_ZERO; - goto exit_label; - } - mag = 1.0 / mag; - RE(h[k]) = CMCONJRE(num, den) * mag; - IM(h[k]) = CMCONJIM(num, den) * mag; - } - } - else - { - // FIR filter if a == NULL - for(k = 0; k < n; k++) - { - RE(jw) = cos(w[k]); - IM(jw) = -sin(w[k]); - res = polyval_cmplx(bc, ord, &jw, 1, h+k); - if(res != RES_OK) - goto exit_label; - } - } - - res = RES_OK; -exit_label: - if(bc) - free(bc); - if(ac) - free(ac); - return res; -} diff --git a/dspl/src/filter_ap.c b/dspl/src/filter_ap.c deleted file mode 100644 index 7f2612f..0000000 --- a/dspl/src/filter_ap.c +++ /dev/null @@ -1,332 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - -#include -#include -#include -#include "dspl.h" - - - - -int butter_ap_zp(int ord, double rp, complex_t *z, int* nz, complex_t *p, int* np) -{ - double alpha; - double theta; - double ep; - int r; - int L; - int ind = 0, k; - - if(rp < 0 || rp == 0) - return ERROR_FILTER_RP; - if(ord < 1) - return ERROR_FILTER_ORD; - if(!z || !p || !nz || !np) - return ERROR_PTR; - - ep = sqrt(pow(10.0, rp*0.1) - 1.0); - r = ord % 2; - L = (int)((ord-r)/2); - - alpha = pow(ep, -1.0/(double)ord); - if(r) - { - RE(p[ind]) = -alpha; - IM(p[ind]) = 0.0; - ind++; - } - for(k = 0; k < L; k++) - { - theta = M_PI*(double)(2*k + 1)/(double)(2*ord); - RE(p[ind]) = RE(p[ind+1]) = -alpha * sin(theta); - IM(p[ind]) = alpha * cos(theta); - IM(p[ind+1]) = -alpha * cos(theta); - ind+=2; - } - return RES_OK; -} - - - - -/************************************************************************************************** -Analog Normalized Butterworth filter -***************************************************************************************************/ -int butter_ap(double Rp, int ord, double* b, double* a) -{ - double ep; - double p[3] = {0.0, 0.0, 1.0}; - double alpha; - double teta; - double *acc = NULL; - int r; - int L; - int n; - int k; - int res; - - if(Rp < 0 || Rp == 0) - return ERROR_FILTER_RP; - if(ord < 1) - return ERROR_FILTER_ORD; - if(!b || !a) - return ERROR_PTR; - - acc = (double*)malloc((ord+1)*sizeof(double)); - - memset(acc, 0, (ord+1)*sizeof(double)); - memset(a, 0, (ord+1)*sizeof(double)); - memset(b, 0, (ord+1)*sizeof(double)); - - - ep = sqrt(pow(10.0, Rp*0.1) - 1.0); - r = ord % 2; - L = (int)((ord-r)/2); - - b[0] = 1.0/ep; - - alpha = pow(ep, -1.0/(double)ord); - /* first pole according to filter order */ - if(r) - { - /* we have one real pole if filter order is odd */ - a[0] = alpha; - a[1] = 1.0; - k = 2; - } - else - { - /* all poles are complex if filter order is even */ - a[0] = 1.0; - k = 1; - } - - /* coeff calculation */ - for(n = 0; n < L; n++) - { - memcpy(acc, a, k*sizeof(double)); - teta = M_PI*(double)(2*n + 1)/(double)(2*ord); - p[0] = alpha * alpha; - p[1] = 2.0 * alpha * sin(teta); - - res = conv(p, 3, acc, k, a); - if(res!=RES_OK) - goto exit_label; - - k+=2; - } - -exit_label: - if(acc) - free(acc); - return res; -} - - - -/************************************************************************************************** -Analog Normalized Chebyshev type 1 filter -***************************************************************************************************/ -int cheby1_ap(double Rp, int ord, double* b, double* a) -{ - double ep; - double alpha; - double beta; - double shbeta; - double chbeta; - double sigma; - double omega; - double gain; - - - double p[3] = {0.0, 0.0, 1.0}; - - double *acc = NULL; - int r; - int L; - int n; - int k; - int res; - - if(Rp < 0 || Rp == 0) - return ERROR_FILTER_RP; - if(ord < 1) - return ERROR_FILTER_ORD; - if(!b || !a) - return ERROR_PTR; - - acc = (double*)malloc((ord+1)*sizeof(double)); - - memset(acc, 0, (ord+1)*sizeof(double)); - memset(a, 0, (ord+1)*sizeof(double)); - memset(b, 0, (ord+1)*sizeof(double)); - - ep = sqrt(pow(10.0, Rp*0.1) - 1.0); - r = ord % 2; - L = (int)((ord-r)/2); - - beta = asinh(1.0/ep)/(double)ord; - chbeta = cosh(beta); - shbeta = sinh(beta); - - /* first pole according to filter order */ - if(r) - { - /* we have one real pole if filter order is odd */ - a[0] = -shbeta; - a[1] = 1.0; - k = 2; - gain = shbeta; - } - else - { - /* all poles are complex if filter order is even */ - a[0] = 1.0; - k = 1; - gain = 1.0/sqrt(1.0 + ep*ep); - } - - - /* coeff calculation */ - for(n = 0; n < L; n++) - { - memcpy(acc, a, k*sizeof(double)); - alpha = M_PI*(double)(2*n + 1)/(double)(2*ord); - sigma = -sin(alpha)*shbeta; - omega = cos(alpha)*chbeta; - - p[0] = sigma * sigma + omega * omega; - gain*=p[0]; - p[1] = -2.0 * sigma; - - res = conv(p, 3, acc, k, a); - if(res!=RES_OK) - goto exit_label; - k+=2; - } - b[0] = gain; - - res = RES_OK; -exit_label: - if(acc) - free(acc); - return res; -} - - - -/************************************************************************************************** -Analog Normalized Chebyshev type II filter -***************************************************************************************************/ -int cheby2_ap(double Rs, int ord, double *b, double *a) -{ - double es, *acc = NULL, *bcc = NULL; - double alpha, beta, sigma, omega, sh, ch, so2; - double p[3] = {0.0, 0.0, 1.0}; - double q[3] = {0.0, 0.0, 1.0}; - double gain; - int r, L, n, kp; - int res; - - if(Rs < 0 || Rs == 0) - return ERROR_FILTER_RS; - if(ord < 1) - return ERROR_FILTER_ORD; - if(!b || !a) - return ERROR_PTR; - - acc = (double*)malloc((ord+1)*sizeof(double)); - bcc = (double*)malloc((ord+1)*sizeof(double)); - - memset(acc, 0, (ord+1)*sizeof(double)); - memset(bcc, 0, (ord+1)*sizeof(double)); - memset(a, 0, (ord+1)*sizeof(double)); - memset(b, 0, (ord+1)*sizeof(double)); - - es = sqrt(pow(10.0, Rs*0.1) - 1.0); - r = ord % 2; - L = (int)((ord-r)/2); - - beta = asinh(es)/(double)ord; - - sh = sinh(beta); - ch = cosh(beta); - - b[0] = 1.0; - /* first pole according to filter order */ - if(r) - { - /* we have one real pole if filter order is odd */ - a[0] = 1.0/sh; - a[1] = 1.0; - kp = 2; - } - else - { - /* all poles are TCOMPLEX if filter order is even */ - a[0] = 1.0; - kp = 1; - } - - /* coeff calculation */ - for(n = 0; n < L; n++) - { - memcpy(acc, a, kp * sizeof(double)); - memcpy(bcc, b, kp * sizeof(double)); - - alpha = M_PI*(double)(2*n+1)/(double)(2*ord); - - sigma = sh*sin(alpha); - omega = ch*cos(alpha); - - so2 = 1.0 / (sigma*sigma+omega*omega); - - q[0] = 1.0/(cos(alpha)*cos(alpha)); - p[0] = so2; - p[1] = 2.0*sigma*so2; - - res = conv(p, 3, acc, kp, a); - if(res!=RES_OK) - goto exit_label; - - res = conv(q, 3, bcc, kp, b); - if(res!=RES_OK) - goto exit_label; - - kp+=2; - - } - gain = b[0] / a[0]; - for(n = 0; n < ord+1; n++) - b[n] /= gain; - -exit_label: - if(acc) - free(acc); - if(bcc) - free(bcc); - return res; - -} - - - - diff --git a/dspl/src/fourier_series.c b/dspl/src/fourier_series.c deleted file mode 100644 index c6eaf0d..0000000 --- a/dspl/src/fourier_series.c +++ /dev/null @@ -1,100 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - - -#include -#include -#include -#include "dspl.h" - - - -int fourier_series_dec(double* t, double* s, int nt, double period, int nw, double* w, complex_t* y) -{ - int k, m; - double dw = M_2PI / period; - complex_t e[2]; - - if(!t || !s || !w || !y) - return ERROR_PTR; - if(nt<1 || nw < 1) - return ERROR_SIZE; - if(period <= 0.0) - return ERROR_NEGATIVE; - - memset(y, 0 , nw*sizeof(complex_t)); - - for(k = 0; k < nw; k++) - { - w[k] = (k - nw/2) * dw; - RE(e[1]) = s[0] * cos(w[k] * t[0]); - IM(e[1]) = -s[0] * sin(w[k] * t[0]); - for(m = 1; m < nt; m++) - { - RE(e[0]) = RE(e[1]); - IM(e[0]) = IM(e[1]); - RE(e[1]) = s[m] * cos(w[k] * t[m]); - IM(e[1]) = - s[m] * sin(w[k] * t[m]); - RE(y[k]) += 0.5 * (RE(e[0]) + RE(e[1])) * (t[m] - t[m-1]); - IM(y[k]) += 0.5 * (IM(e[0]) + IM(e[1])) * (t[m] - t[m-1]); - } - RE(y[k]) /= period; - IM(y[k]) /= period; - } - - if(!(nw%2)) - RE(y[0]) = RE(y[1]) = 0.0; - - return RES_OK; -} - - - - - - -int fourier_series_rec(double* w, complex_t* s, int nw, double *t, int nt, complex_t* y) -{ - int k, m; - complex_t e; - - if(!t || !s || !w || !y) - return ERROR_PTR; - if(nt<1 || nw < 1) - return ERROR_SIZE; - - memset(y, 0, nt*sizeof(complex_t)); - - - for(k = 0; k < nw; k++) - { - for(m = 1; m < nt; m++) - { - RE(e) = cos(w[k] * t[m]); - IM(e) = sin(w[k] * t[m]); - - RE(y[k]) += CMRE(s[k], e); - IM(y[k]) += CMIM(s[k], e); - - } - } - return RES_OK; - -} diff --git a/dspl/src/goertzel.c b/dspl/src/goertzel.c deleted file mode 100644 index 94f7d54..0000000 --- a/dspl/src/goertzel.c +++ /dev/null @@ -1,110 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - - -#include -#include -#include "dspl.h" - - -/************************************************************************************************** -Goertzel algorithm for real vector -**************************************************************************************************/ -int goertzel(double *x, int n, int *ind, int k, complex_t *y) -{ - - int m, p; - double wR, wI; - double alpha; - double v[3]; - - if(!x || !y || !ind) - return ERROR_PTR; - - if(n < 1 || k < 1) - return ERROR_SIZE; - - for(p = 0; p < k; p++) - { - wR = cos(M_2PI * (double)ind[p] / (double)n); - wI = sin(M_2PI * (double)ind[p] / (double)n); - - alpha = 2.0 * wR; - v[0] = v[1] = v[2] = 0.0; - - for(m = 0; m < n; m++) - { - v[2] = v[1]; - v[1] = v[0]; - v[0] = x[m]+alpha*v[1] - v[2]; - } - RE(y[p]) = wR * v[0] - v[1]; - IM(y[p]) = wI * v[0]; - } - - return RES_OK; -} - - - - - -/************************************************************************************************** -Goertzel algorithm for complex vector -**************************************************************************************************/ -int goertzel_cmplx(complex_t *x, int n, int *ind, int k, complex_t *y) -{ - - int m, p; - complex_t w; - double alpha; - complex_t v[3]; - - if(!x || !y || !ind) - return ERROR_PTR; - - if(n < 1 || k < 1) - return ERROR_SIZE; - - for(p = 0; p < k; p++) - { - RE(w) = cos(M_2PI * (double)ind[p] / (double)n); - IM(w) = sin(M_2PI * (double)ind[p] / (double)n); - - alpha = 2.0 * RE(w); - memset(v, 0, 3*sizeof(complex_t)); - - for(m = 0; m < n; m++) - { - RE(v[2]) = RE(v[1]); - RE(v[1]) = RE(v[0]); - RE(v[0]) = RE(x[m]) + alpha * RE(v[1]) - RE(v[2]); - - IM(v[2]) = IM(v[1]); - IM(v[1]) = IM(v[0]); - IM(v[0]) = IM(x[m]) + alpha * IM(v[1]) - IM(v[2]); - } - - RE(y[p]) = CMRE(w, v[0]) - RE(v[1]); - IM(y[p]) = CMIM(w, v[0]) - IM(v[1]); - } - - return RES_OK; -} diff --git a/dspl/src/inout.c b/dspl/src/inout.c deleted file mode 100644 index 79f97c4..0000000 --- a/dspl/src/inout.c +++ /dev/null @@ -1,129 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - - -#include -#include -#include "dspl.h" - - -/************************************************************************************************** -Write a real array to the binary file "fn" -***************************************************************************************************/ -int writebin(void* x, int n, int dtype, char* fn) -{ - int k, res; - FILE* pFile = NULL; - - if(!x) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - if(!fn) - return ERROR_FNAME; - - pFile = fopen(fn, "wb"); - if(pFile == NULL) - return ERROR_FOPEN; - - - if(fwrite(&dtype, sizeof(int), 1, pFile) != 1) - { - res = ERROR_FWRITE_SIZE; - goto exit_label; - } - - - if(fwrite(&n, sizeof(int), 1, pFile) != 1) - { - res = ERROR_FWRITE_SIZE; - goto exit_label; - } - - k = 1; - if(fwrite(&k, sizeof(int), 1, pFile) != 1) - { - res = ERROR_FWRITE_SIZE; - goto exit_label; - }; - - switch(dtype) - { - case DAT_DOUBLE: - if(fwrite((double*)x, sizeof(double), n, pFile) != n) - { - res = ERROR_FWRITE_SIZE; - goto exit_label; - } - break; - case DAT_COMPLEX: - if(fwrite((complex_t*)x, sizeof(complex_t), n, pFile) != n) - { - res = ERROR_FWRITE_SIZE; - goto exit_label; - } - break; - default: - res = ERROR_DAT_TYPE; - goto exit_label; - } - res = RES_OK; -exit_label: - if(pFile) - fclose(pFile); - return res; -} - - - - - - -/************************************************************************************************** -Write a real arrays to the text file "fn" -***************************************************************************************************/ -int writetxt(double* x, double *y, int n, char* fn) -{ - int k; - FILE* pFile = NULL; - - if(!x) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - if(!fn) - return ERROR_FNAME; - - pFile = fopen(fn, "w+"); - if(pFile == NULL) - return ERROR_FOPEN; - - if(y) - for(k = 0; k < n; k++) - fprintf(pFile, "%+.12E\t%+.12E\n", x[k], y[k]); - else - for(k = 0; k < n; k++) - fprintf(pFile, "%+.12E\n", x[k]); - - fclose(pFile); - return RES_OK; -} - - diff --git a/dspl/src/math.c b/dspl/src/math.c deleted file mode 100644 index b8504ce..0000000 --- a/dspl/src/math.c +++ /dev/null @@ -1,33 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - - -#include -#include -#include "dspl.h" - - - -double dmod (double x, double y) -{ - if(y == 0.0) - return x; - return x - floor(x/y) * y; -} diff --git a/dspl/src/polyval.c b/dspl/src/polyval.c deleted file mode 100644 index a252404..0000000 --- a/dspl/src/polyval.c +++ /dev/null @@ -1,136 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - - -#include -#include "dspl.h" - - - - - -/************************************************************************************************** -Real polynomial evaluation -***************************************************************************************************/ -int polyval(double* a, int ord, double* x, int n, double* y) -{ - int k, m; - - if(!a || !x || !y) - return ERROR_PTR; - if(ord<0) - return ERROR_POLY_ORD; - if(n<1) - return ERROR_SIZE; - - for(k = 0; k < n; k++) - { - y[k] = a[ord]; - for(m = ord-1; m>-1; m--) - y[k] = y[k]*x[k] + a[m]; - } - return RES_OK; -} - - - - - -/************************************************************************************************** -Complex polynomial evaluation -***************************************************************************************************/ -int polyval_cmplx(complex_t* a, int ord, complex_t* x, int n, complex_t* y) -{ - int k, m; - complex_t t; - - if(!a || !x || !y) - return ERROR_PTR; - if(ord<0) - return ERROR_POLY_ORD; - if(n<1) - return ERROR_SIZE; - - for(k = 0; k < n; k++) - { - RE(y[k]) = RE(a[ord]); - IM(y[k]) = IM(a[ord]); - for(m = ord-1; m>-1; m--) - { - RE(t) = CMRE(y[k], x[k]); - IM(t) = CMIM(y[k], x[k]); - RE(y[k]) = RE(t) + RE(a[m]); - IM(y[k]) = IM(t) + IM(a[m]); - } - } - return RES_OK; -} - - - - - - - -/************************************************************************************************** -polynomial zeros to coeff reecalc -***************************************************************************************************/ -int poly_z2a(complex_t *z, int nz, int ord, complex_t *a) -{ - int k, ind, res; - complex_t x[2]; - - if(!z || !a) - return ERROR_PTR; - if(ord<0) - return ERROR_POLY_ORD; - if(nz<1 || nz > ord) - return ERROR_SIZE; - - memset(a, 0, (ord+1) * sizeof(complex_t)); - RE(a[0]) = 1.0; - IM(a[0]) = 0.0; - - RE(x[1]) = 1.0; - IM(x[1]) = 0.0; - - - ind = 1; - for(k=0; k < nz; k++) - { - RE(x[0]) = -RE(z[k]); - IM(x[0]) = -IM(z[k]); - res = conv_cmplx(a, ind, x, 2, a); - if(res!=RES_OK) - return res; - ind++; - } - return RES_OK; -} - - - - - - - - - - diff --git a/dspl/src/randgen.c b/dspl/src/randgen.c deleted file mode 100644 index 3f32339..0000000 --- a/dspl/src/randgen.c +++ /dev/null @@ -1,100 +0,0 @@ - - -#include -#include -#include - -#include "dspl.h" - - -#define DSPL_RAND_MOD_X1 2147483647 -#define DSPL_RAND_MOD_X2 2145483479 - - -/************************************************************************************************** -Uniform random numbers generator -***************************************************************************************************/ -int randu(double* x, int n) -{ - int k,m; - unsigned int x1[4], x2[4], y; - - if(!x) - return ERROR_PTR; - if(n<1) - return ERROR_SIZE; - - x1[1] = rand(); - x2[1] = rand(); - x1[2] = rand(); - x2[2] = rand(); - x1[3] = rand(); - x2[3] = rand(); - for(k = 0; k 0; m--) - { - x1[m] = x1[m-1]; - x2[m] = x2[m-1]; - } - - x[k] = (double)y/DSPL_RAND_MOD_X1; - } - - return RES_OK; -} - - - - - -/************************************************************************************************** -Gaussian random numbers generator -***************************************************************************************************/ -int randn(double* x, int n, double mu, double sigma) -{ - int k, m; - double x1[128], x2[128]; - int res; - if(!x) - return ERROR_PTR; - - if(n<1) - return ERROR_SIZE; - - if(sigma < 0.0) - return ERROR_RAND_SIGMA; - - k=0; - while(k < n) - { - res = randu(x1, 128); - if(res != RES_OK) - goto exit_label; - - res = randu(x2, 128); - if(res != RES_OK) - goto exit_label; - m = 0 ; - while(k. -*/ - - -#include -#include -#include -#include "dspl.h" - -#define DSPL_FARROW_LAGRANGE_COEFF 0.16666666666666666666666666666667 - - - -/************************************************************************************************** -Farrow resampler based on the cubic Lagrange polynomials -***************************************************************************************************/ -int farrow_lagrange(double *s, int n, double p, double q, double frd, double **y, int *ny) -{ - double a[4]; - double t, x, dt; - int ind, k, res; - double g[4]; - double *z; - - if(!s || !y) - return ERROR_PTR; - - if(n<1) - return ERROR_SIZE; - - if(p <= 0.0 || q <= 0.0) - return ERROR_RESAMPLE_RATIO; - - if(frd <= -1.0 || frd >= 1.0) - return ERROR_RESAMPLE_FRAC_DELAY; - - dt = q/p; - - if((*ny) != (int)((double)(n-1)/dt)+1 || !(*y)) - { - - *ny = (int)((double)(n-1)/dt)+1; - (*y) = (double*)realloc((*y), (*ny)*sizeof(double)); - } - - t = -frd; - k = 0; - while(k < (*ny)) - { - ind = floor(t)+1; - x = t - (double)ind; - ind-=2; - if(ind < 0) - { - memset(g, 0, 4*sizeof(double)); - if(ind > (-3)) - memcpy(g-ind, s, (4+ind)*sizeof(double)); - z = g; - } - else - { - if(ind < n-3) - z = s+ind; - else - { - memset(g, 0, 4*sizeof(double)); - if((n-ind)>0) - memcpy(g, s+ind, (n-ind)*sizeof(double)); - z = g; - } - } - a[0] = z[2]; - a[3] = DSPL_FARROW_LAGRANGE_COEFF*(z[3] -z[0]) + 0.5*(z[1] - z[2]); - a[1] = 0.5*(z[3] - z[1])-a[3]; - a[2] = z[3] - z[2] -a[3]-a[1]; - - res = polyval(a, 3, &x, 1, (*y)+k); - - if(res != RES_OK) - goto exit_label; - t+=dt; - k++; - } - -exit_label: - return res; -} - - - - - -/************************************************************************************************** -Farrow resampler based on the cubic splines -**************************************************************************************************/ -int farrow_spline(double *s, int n, double p, double q, double frd, double **y, int *ny) -{ - double a[4]; - double t, x, dt; - int ind, k, res; - double g[4]; - double *z; - - if(!s || !y) - return ERROR_PTR; - - if(n<1) - return ERROR_SIZE; - - if(p <= 0.0 || q <= 0.0) - return ERROR_RESAMPLE_RATIO; - - if(frd <= -1.0 || frd >= 1.0) - return ERROR_RESAMPLE_FRAC_DELAY; - - dt = q/p; - - if((*ny) != (int)((double)(n-1)/dt)+1 || !(*y)) - { - - *ny = (int)((double)(n-1)/dt)+1; - (*y) = (double*)realloc((*y), (*ny)*sizeof(double)); - } - - t = -frd; - k = 0; - while(k < (*ny)) - { - ind = floor(t)+1; - x = t - (double)ind; - ind-=2; - if(ind < 0) - { - memset(g, 0, 4*sizeof(double)); - if(ind > (-3)) - memcpy(g-ind, s, (4+ind)*sizeof(double)); - z = g; - } - else - { - if(ind < n-3) - z = s+ind; - else - { - memset(g, 0, 4*sizeof(double)); - if((n-ind)>0) - memcpy(g, s+ind, (n-ind)*sizeof(double)); - z = g; - } - } - a[0] = z[2]; - a[1] = 0.5*(z[3] - z[1]); - a[3] = 2.0*(z[1] - z[2]) + a[1] + 0.5*(z[2] - z[0]); - a[2] = z[1] - z[2] +a[3] + a[1]; - - res = polyval(a, 3, &x, 1, (*y)+k); - - if(res != RES_OK) - goto exit_label; - t+=dt; - k++; - } - -exit_label: - return res; -} - - - diff --git a/dspl/src/signals.c b/dspl/src/signals.c deleted file mode 100644 index ebe3c5c..0000000 --- a/dspl/src/signals.c +++ /dev/null @@ -1,56 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - - -#include -#include -#include -#include "dspl.h" - - -/************************************************************************************************** -Rectangle pulse signal -***************************************************************************************************/ -int signal_pimp(double* t, size_t n, double amp, double tau, double dt, double period, double* y) -{ - int k; - double ll, lr, p2, tp; - - if(!t || !y) - return ERROR_PTR; - if(n < 1) - return ERROR_SIZE; - if(tau < 0.0 || period < 0.0) - return ERROR_NEGATIVE; - - - ll = -0.5 * tau; - lr = 0.5 * tau; - p2 = period*0.5; - for(k = 0; k < n; k++) - { - tp = dmod(t[k] - dt + p2, period) - p2; - y[k] = (tp < ll || tp > lr) ? 0.0 : amp; - } - - return RES_OK; -} - - diff --git a/dspl/src/trapint.c b/dspl/src/trapint.c deleted file mode 100644 index aa37601..0000000 --- a/dspl/src/trapint.c +++ /dev/null @@ -1,70 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - - -#include -#include -#include "dspl.h" - - -/************************************************************************************************** -Numerical integration of real data using the trapezoidal method. -***************************************************************************************************/ -int trapint(double* x, double* y, int n, double* sum) -{ - int k; - - if(!x || !y) - return ERROR_PTR; - if(n<2) - return ERROR_SIZE; - *sum = 0.0; - - for(k = 1; k < n; k++) - *sum += 0.5 * (x[k] - x[k-1]) * (y[k] + y[k-1]); - - return RES_OK; -} - - - - -/************************************************************************************************** -Numerical integration of complex data using the trapezoidal method. -***************************************************************************************************/ -int trapint_cmplx(double* x, complex_t* y, int n, complex_t* sum) -{ - - int k; - double dx; - if(!x || !y) - return ERROR_PTR; - if(n<2) - return ERROR_SIZE; - RE(*sum) = IM(*sum) = 0.0; - - for(k = 1; k < n; k++) - { - dx = 0.5 * (x[k] - x[k-1]); - RE(*sum) += dx * (RE(y[k]) + RE(y[k-1])); - IM(*sum) += dx * (IM(y[k]) + IM(y[k-1])); - } - return RES_OK; -} diff --git a/dspl/src/win.c b/dspl/src/win.c deleted file mode 100644 index ff1bbb2..0000000 --- a/dspl/src/win.c +++ /dev/null @@ -1,520 +0,0 @@ -/* -* Copyright (c) 2015-2017 Sergey Bakhurin -* Digital Signal Processing Library [http://dsplib.org] -* -* This file is part of DSPL. -* -* 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. -* -* DSPL 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 Foobar. If not, see . -*/ - - -#include -#include -#include "dspl.h" - - - -int win_bartlett (double *w, int n, int win_type); -int win_bartlett_hann (double *w, int n, int win_type); -int win_blackman (double *w, int n, int win_type); -int win_blackman_harris (double *w, int n, int win_type); -int win_blackman_nuttall(double *w, int n, int win_type); -int win_cos (double *w, int n, int win_type); -int win_flat_top (double *w, int n, int win_type); -int win_gaussian (double *w, int n, int win_type, double sigma); -int win_hamming (double *w, int n, int win_type); -int win_hann (double *w, int n, int win_type); -int win_lanczos (double *w, int n, int win_type); -int win_nuttall (double *w, int n, int win_type); -int win_rect (double *w, int n); - - - - - -/************************************************************************************************** -Window function -***************************************************************************************************/ -int window(double* w, int n, int win_type, double param) -{ - - switch(win_type & DSPL_WIN_MASK) - { - case DSPL_WIN_BARTLETT : return win_bartlett (w, n, win_type); - case DSPL_WIN_BARTLETT_HANN : return win_bartlett_hann (w, n, win_type); break; - case DSPL_WIN_BLACKMAN : return win_blackman (w, n, win_type); break; - case DSPL_WIN_BLACKMAN_HARRIS : return win_blackman_harris (w, n, win_type); break; - case DSPL_WIN_BLACKMAN_NUTTALL : return win_blackman_nuttall (w, n, win_type); break; - case DSPL_WIN_FLAT_TOP : return win_flat_top (w, n, win_type); break; - case DSPL_WIN_GAUSSIAN : return win_gaussian (w, n, win_type, param); break; - case DSPL_WIN_HAMMING : return win_hamming (w, n, win_type); break; - case DSPL_WIN_HANN : return win_hann (w, n, win_type); break; - case DSPL_WIN_LANCZOS : return win_lanczos (w, n, win_type); break; - case DSPL_WIN_NUTTALL : return win_nuttall (w, n, win_type); break; - case DSPL_WIN_RECT : return win_rect (w, n); break; - case DSPL_WIN_COS : return win_cos (w, n, win_type); break; - default : return ERROR_WIN_TYPE; - } - return RES_OK; -} - - - -/************************************************************************************************** -Barlett window function -***************************************************************************************************/ -int win_bartlett(double *w, int n, int win_type) -{ - double x = 0.0; - int i; - - if(!w) - return ERROR_PTR; - if(n<2) - return ERROR_SIZE; - - switch(win_type & DSPL_WIN_SYM_MASK) - { - case DSPL_WIN_SYMMETRIC: x = (double)(n-1); break; - case DSPL_WIN_PERIODIC : x = (double)n; break; - default: return ERROR_WIN_SYM; - } - - for(i = 0; i < n; i++) - { - w[i] = 2.0 / x * (x * 0.5-fabs((double)i - x * 0.5)); - } - return RES_OK; -} - - - - - -/************************************************************************************************** -Barlett - Hann window function -***************************************************************************************************/ -int win_bartlett_hann(double *w, int n, int win_type) -{ - double y; - double x = 0.0; - int i; - - if(!w) - return ERROR_PTR; - if(n<2) - return ERROR_SIZE; - - switch(win_type & DSPL_WIN_SYM_MASK) - { - case DSPL_WIN_SYMMETRIC: x = 1.0/(double)(n-1); break; - case DSPL_WIN_PERIODIC : x = 1.0/(double)n; break; - default: return ERROR_WIN_SYM; - } - - y = 0.0; - for(i = 0; i #endif //WIN_OS @@ -35,82 +34,82 @@ #ifndef BUILD_LIB -p_blas_dscal blas_dscal ; - -p_cheby_poly1 cheby_poly1 ; -p_cheby_poly2 cheby_poly2 ; +p_cheby_poly1 cheby_poly1 ; +p_cheby_poly2 cheby_poly2 ; +p_conv conv ; +p_conv_cmplx conv_cmplx ; p_dft dft ; p_dft_cmplx dft_cmplx ; - -p_fft_create fft_create ; -p_fft_destroy fft_destroy ; -p_fft_cmplx fft_cmplx ; -p_fft_shift fft_shift ; - - +p_filter_iir filter_iir ; +p_linspace linspace ; +p_logspace logspace ; #endif //BUILD_LIB +#ifdef WIN_OS +#define LOAD_FUNC(fn) \ + fname = #fn;\ + fn = (p_##fn)GetProcAddress(handle, fname);\ + if(! fn) goto exit_label; +#endif + + + +#ifdef LINUX_OS +#define LOAD_FUNC(fn) \ + fname = #fn;\ + fn = (p_##fn)dlsym(handle, fname);\ + if ((error = dlerror()) != NULL) goto exit_label +#endif + void* dspl_load() { + char* fname; #ifdef WIN_OS - HINSTANCE handle; - char* fname; + HINSTANCE handle; handle = LoadLibrary(TEXT("libdspl.dll")); if (!handle) { printf("libdspl.dll loading ERROR!\n"); return NULL; } - - - fname = "blas_dscal"; - blas_dscal = (p_blas_dscal)GetProcAddress(handle, fname); - if(!blas_dscal) goto exit_label; + #endif //WIN_OS - fname = "cheby_poly1"; - cheby_poly1 = (p_cheby_poly1)GetProcAddress(handle, fname); - if(!cheby_poly1) goto exit_label; + #ifdef LINUX_OS + char* error; + void *handle; + // open the *.so + handle = dlopen ("./libdspl.so", RTLD_LAZY); + if (!handle) + { + printf("libdspl.so loading ERROR!\n"); + return NULL; + } + #endif //LINUX_OS - fname = "cheby_poly2"; - cheby_poly2 = (p_cheby_poly2)GetProcAddress(handle, fname); - if(!cheby_poly2) goto exit_label; - - fname = "dft"; - dft = (p_dft)GetProcAddress(handle, fname); - if(!dft) goto exit_label; - - fname = "dft_cmplx"; - dft_cmplx = (p_dft_cmplx)GetProcAddress(handle, fname); - if(!dft_cmplx) goto exit_label; - - - fname = "fft_create"; - fft_create = (p_fft_create)GetProcAddress(handle, fname); - if(!fft_create) goto exit_label; - - fname = "fft_destroy"; - fft_destroy = (p_fft_destroy)GetProcAddress(handle, fname); - if(!fft_destroy) goto exit_label; - - - fname = "fft_cmplx"; - fft_cmplx = (p_fft_cmplx)GetProcAddress(handle, fname); - if(!fft_cmplx) goto exit_label; - - fname = "fft_shift"; - fft_shift = (p_fft_shift)GetProcAddress(handle, fname); - if(!fft_shift) goto exit_label; - - - return (void*)handle; + + + LOAD_FUNC(cheby_poly1); + LOAD_FUNC(cheby_poly2); + LOAD_FUNC(conv); + LOAD_FUNC(conv_cmplx); + LOAD_FUNC(dft); + LOAD_FUNC(dft_cmplx); + LOAD_FUNC(filter_iir); + LOAD_FUNC(linspace); + LOAD_FUNC(logspace); + + + + #ifdef WIN_OS + return (void*)handle; exit_label: printf("function %s loading ERROR!\n", fname); if(handle) @@ -119,75 +118,14 @@ void* dspl_load() #endif //WIN_OS - - - - - #ifdef LINUX_OS - char* error; - void *handle; - char* fname; - - // open the *.so - handle = dlopen ("./libdspl.so", RTLD_LAZY); - if (!handle) - { - printf("libdspl.so loading ERROR!\n"); - return NULL; - } - - - fname = "blas_dscal"; - blas_dscal = (p_blas_dscal)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - fname = "cheby_poly1"; - cheby_poly1 = (p_cheby_poly1)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - fname = "cheby_poly2"; - cheby_poly2 = (p_cheby_poly2)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - fname = "dft"; - dft = (p_dft)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - fname = "dft_cmplx"; - dft_cmplx = (p_dft_cmplx)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - fname = "fft_create"; - fft_create = (p_fft_create)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - fname = "fft_destroy"; - fft_destroy = (p_fft_destroy)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - fname = "fft_cmplx"; - fft_cmplx = (p_fft_cmplx)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - fname = "fft_shift"; - fft_shift = (p_fft_shift)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - return handle; - + #ifdef LINUX_OS + return handle; exit_label: printf("function %s loading ERROR!\n", fname); if(handle) dlclose(handle); return NULL; - #endif //LINUX_OS - + #endif //LINUX_OS } diff --git a/include/dspl.h b/include/dspl.h index 3703878..5e87762 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -162,66 +162,45 @@ extern "C" { #endif -// Declare DSPL_API for Windows OS -#ifdef BUILD_LIB +#ifdef BUILD_LIB +// Declare DSPL_API for Windows OS #ifdef WIN_OS #define DSPL_API __declspec(dllexport) #endif // WIN_OS - +// Declare DSPL_API for LINUX OS #ifdef LINUX_OS #define DSPL_API #endif //LINUX_OS - - -int DSPL_API blas_dscal(int n, double a, double* x, int incx); - -int DSPL_API cheby_poly1(double* x, int n, int ord, double* y); -int DSPL_API cheby_poly2(double* x, int n, int ord, double* y); - -int DSPL_API dft (double* x, int n, complex_t *y); -int DSPL_API dft_cmplx (complex_t* x, int n, complex_t *y); - - -int DSPL_API fft_create ( fft_t *pfft, int n); -void DSPL_API fft_destroy (fft_t *pfft); -int DSPL_API fft_cmplx (complex_t *x, int n, fft_t* pfft, complex_t* y); -int DSPL_API fft_shift (double* x, int n, double* y); - - -#else //BUILD_LIB - - - -typedef int (*p_blas_dscal ) (int n, double a, double* x, int incx); - -typedef int (*p_cheby_poly1 ) (double* x, int n, int ord, double* y); -typedef int (*p_cheby_poly2 ) (double* x, int n, int ord, double* y); - -typedef int (*p_dft ) (double* x, int n, complex_t *y); -typedef int (*p_dft_cmplx ) (complex_t* x, int n, complex_t *y); -typedef int (*p_fft_create ) (fft_t *pfft, int n); -typedef void(*p_fft_destroy ) (fft_t *pfft); -typedef int (*p_fft_cmplx ) (complex_t *x, int n, fft_t* pfft, complex_t* y); -typedef int (*p_fft_shift ) (double* x, int n, double* y); - - - -extern p_blas_dscal blas_dscal ; - -extern p_cheby_poly1 cheby_poly1 ; -extern p_cheby_poly2 cheby_poly2 ; - -extern p_dft dft ; -extern p_dft_cmplx dft_cmplx ; - -extern p_fft_create fft_create ; -extern p_fft_destroy fft_destroy ; -extern p_fft_cmplx fft_cmplx ; -extern p_fft_shift fft_shift ; - #endif //BUILD_DLL +#define COMMA , + + +#ifdef BUILD_LIB +#define DECLARE_FUNC(type, fn, param)\ + type DSPL_API fn(param);\ + +#endif + +#ifndef BUILD_LIB +#define DECLARE_FUNC(type, fn, param)\ + typedef type (*p_##fn)(param);\ + extern p_##fn fn; + +#endif + + +DECLARE_FUNC(int, cheby_poly1, double* COMMA int COMMA int COMMA double*); +DECLARE_FUNC(int, cheby_poly2, double* COMMA int COMMA int COMMA double*); +DECLARE_FUNC(int, conv, double* COMMA int COMMA double* COMMA int COMMA double*); +DECLARE_FUNC(int, conv_cmplx, complex_t* COMMA int COMMA complex_t* COMMA int COMMA complex_t*); +DECLARE_FUNC(int, dft, double* COMMA int COMMA complex_t*); +DECLARE_FUNC(int, dft_cmplx, complex_t* COMMA int COMMA complex_t*); +DECLARE_FUNC(int, filter_iir, double* COMMA double* COMMA int COMMA double* COMMA int COMMA double*); +DECLARE_FUNC(int, linspace, double COMMA double COMMA int COMMA int COMMA double*); +DECLARE_FUNC(int, logspace, double COMMA double COMMA int COMMA int COMMA double*); + #ifdef __cplusplus } diff --git a/include/fftw3.h b/include/fftw3.h deleted file mode 100644 index bfde54e..0000000 --- a/include/fftw3.h +++ /dev/null @@ -1,514 +0,0 @@ -/* - * Copyright (c) 2003, 2007-14 Matteo Frigo - * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology - * - * The following statement of license applies *only* to this header file, - * and *not* to the other files distributed with FFTW or derived therefrom: - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS - * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED - * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY - * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE - * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS - * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, - * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ - -/***************************** NOTE TO USERS ********************************* - * - * THIS IS A HEADER FILE, NOT A MANUAL - * - * If you want to know how to use FFTW, please read the manual, - * online at http://www.fftw.org/doc/ and also included with FFTW. - * For a quick start, see the manual's tutorial section. - * - * (Reading header files to learn how to use a library is a habit - * stemming from code lacking a proper manual. Arguably, it's a - * *bad* habit in most cases, because header files can contain - * interfaces that are not part of the public, stable API.) - * - ****************************************************************************/ - -#ifndef FFTW3_H -#define FFTW3_H - -#include - -#ifdef __cplusplus -extern "C" -{ -#endif /* __cplusplus */ - -/* If is included, use the C99 complex type. Otherwise - define a type bit-compatible with C99 complex */ -#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) -# define FFTW_DEFINE_COMPLEX(R, C) typedef R _Complex C -#else -# define FFTW_DEFINE_COMPLEX(R, C) typedef R C[2] -#endif - -#define FFTW_CONCAT(prefix, name) prefix ## name -#define FFTW_MANGLE_DOUBLE(name) FFTW_CONCAT(fftw_, name) -#define FFTW_MANGLE_FLOAT(name) FFTW_CONCAT(fftwf_, name) -#define FFTW_MANGLE_LONG_DOUBLE(name) FFTW_CONCAT(fftwl_, name) -#define FFTW_MANGLE_QUAD(name) FFTW_CONCAT(fftwq_, name) - -/* IMPORTANT: for Windows compilers, you should add a line - #define FFTW_DLL - here and in kernel/ifftw.h if you are compiling/using FFTW as a - DLL, in order to do the proper importing/exporting, or - alternatively compile with -DFFTW_DLL or the equivalent - command-line flag. This is not necessary under MinGW/Cygwin, where - libtool does the imports/exports automatically. */ -#if defined(FFTW_DLL) && (defined(_WIN32) || defined(__WIN32__)) - /* annoying Windows syntax for shared-library declarations */ -# if defined(COMPILING_FFTW) /* defined in api.h when compiling FFTW */ -# define FFTW_EXTERN extern __declspec(dllexport) -# else /* user is calling FFTW; import symbol */ -# define FFTW_EXTERN extern __declspec(dllimport) -# endif -#else -# define FFTW_EXTERN extern -#endif - -/* specify calling convention (Windows only) */ -#if defined(_WIN32) || defined(__WIN32__) -# define FFTW_CDECL __cdecl -#else -# define FFTW_CDECL -#endif - -enum fftw_r2r_kind_do_not_use_me { - FFTW_R2HC=0, FFTW_HC2R=1, FFTW_DHT=2, - FFTW_REDFT00=3, FFTW_REDFT01=4, FFTW_REDFT10=5, FFTW_REDFT11=6, - FFTW_RODFT00=7, FFTW_RODFT01=8, FFTW_RODFT10=9, FFTW_RODFT11=10 -}; - -struct fftw_iodim_do_not_use_me { - int n; /* dimension size */ - int is; /* input stride */ - int os; /* output stride */ -}; - -#include /* for ptrdiff_t */ -struct fftw_iodim64_do_not_use_me { - ptrdiff_t n; /* dimension size */ - ptrdiff_t is; /* input stride */ - ptrdiff_t os; /* output stride */ -}; - -typedef void (FFTW_CDECL *fftw_write_char_func_do_not_use_me)(char c, void *); -typedef int (FFTW_CDECL *fftw_read_char_func_do_not_use_me)(void *); - -/* - huge second-order macro that defines prototypes for all API - functions. We expand this macro for each supported precision - - X: name-mangling macro - R: real data type - C: complex data type -*/ - -#define FFTW_DEFINE_API(X, R, C) \ - \ -FFTW_DEFINE_COMPLEX(R, C); \ - \ -typedef struct X(plan_s) *X(plan); \ - \ -typedef struct fftw_iodim_do_not_use_me X(iodim); \ -typedef struct fftw_iodim64_do_not_use_me X(iodim64); \ - \ -typedef enum fftw_r2r_kind_do_not_use_me X(r2r_kind); \ - \ -typedef fftw_write_char_func_do_not_use_me X(write_char_func); \ -typedef fftw_read_char_func_do_not_use_me X(read_char_func); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(execute)(const X(plan) p); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft)(int rank, const int *n, \ - C *in, C *out, int sign, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft_1d)(int n, C *in, C *out, int sign, \ - unsigned flags); \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft_2d)(int n0, int n1, \ - C *in, C *out, int sign, unsigned flags); \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft_3d)(int n0, int n1, int n2, \ - C *in, C *out, int sign, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_many_dft)(int rank, const int *n, \ - int howmany, \ - C *in, const int *inembed, \ - int istride, int idist, \ - C *out, const int *onembed, \ - int ostride, int odist, \ - int sign, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru_dft)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - C *in, C *out, \ - int sign, unsigned flags); \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru_split_dft)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *ri, R *ii, R *ro, R *io, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru64_dft)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - C *in, C *out, \ - int sign, unsigned flags); \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru64_split_dft)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *ri, R *ii, R *ro, R *io, \ - unsigned flags); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(execute_dft)(const X(plan) p, C *in, C *out); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(execute_split_dft)(const X(plan) p, R *ri, R *ii, \ - R *ro, R *io); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_many_dft_r2c)(int rank, const int *n, \ - int howmany, \ - R *in, const int *inembed, \ - int istride, int idist, \ - C *out, const int *onembed, \ - int ostride, int odist, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft_r2c)(int rank, const int *n, \ - R *in, C *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft_r2c_1d)(int n,R *in,C *out,unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft_r2c_2d)(int n0, int n1, \ - R *in, C *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft_r2c_3d)(int n0, int n1, \ - int n2, \ - R *in, C *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_many_dft_c2r)(int rank, const int *n, \ - int howmany, \ - C *in, const int *inembed, \ - int istride, int idist, \ - R *out, const int *onembed, \ - int ostride, int odist, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft_c2r)(int rank, const int *n, \ - C *in, R *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft_c2r_1d)(int n,C *in,R *out,unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft_c2r_2d)(int n0, int n1, \ - C *in, R *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_dft_c2r_3d)(int n0, int n1, \ - int n2, \ - C *in, R *out, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru_dft_r2c)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *in, C *out, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru_dft_c2r)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - C *in, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru_split_dft_r2c)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *in, R *ro, R *io, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru_split_dft_c2r)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *ri, R *ii, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru64_dft_r2c)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *in, C *out, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru64_dft_c2r)(int rank, \ - const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - C *in, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru64_split_dft_r2c)(int rank, const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *in, R *ro, R *io, \ - unsigned flags); \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru64_split_dft_c2r)(int rank, const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *ri, R *ii, R *out, \ - unsigned flags); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(execute_dft_r2c)(const X(plan) p, R *in, C *out); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(execute_dft_c2r)(const X(plan) p, C *in, R *out); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(execute_split_dft_r2c)(const X(plan) p, \ - R *in, R *ro, R *io); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(execute_split_dft_c2r)(const X(plan) p, \ - R *ri, R *ii, R *out); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_many_r2r)(int rank, const int *n, \ - int howmany, \ - R *in, const int *inembed, \ - int istride, int idist, \ - R *out, const int *onembed, \ - int ostride, int odist, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_r2r)(int rank, const int *n, R *in, R *out, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_r2r_1d)(int n, R *in, R *out, \ - X(r2r_kind) kind, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_r2r_2d)(int n0, int n1, R *in, R *out, \ - X(r2r_kind) kind0, X(r2r_kind) kind1, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_r2r_3d)(int n0, int n1, int n2, \ - R *in, R *out, X(r2r_kind) kind0, \ - X(r2r_kind) kind1, X(r2r_kind) kind2, \ - unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru_r2r)(int rank, const X(iodim) *dims, \ - int howmany_rank, \ - const X(iodim) *howmany_dims, \ - R *in, R *out, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN X(plan) \ -FFTW_CDECL X(plan_guru64_r2r)(int rank, const X(iodim64) *dims, \ - int howmany_rank, \ - const X(iodim64) *howmany_dims, \ - R *in, R *out, \ - const X(r2r_kind) *kind, unsigned flags); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(execute_r2r)(const X(plan) p, R *in, R *out); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(destroy_plan)(X(plan) p); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(forget_wisdom)(void); \ -FFTW_EXTERN void \ -FFTW_CDECL X(cleanup)(void); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(set_timelimit)(double t); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(plan_with_nthreads)(int nthreads); \ - \ -FFTW_EXTERN int \ -FFTW_CDECL X(init_threads)(void); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(cleanup_threads)(void); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(make_planner_thread_safe)(void); \ - \ -FFTW_EXTERN int \ -FFTW_CDECL X(export_wisdom_to_filename)(const char *filename); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(export_wisdom_to_file)(FILE *output_file); \ - \ -FFTW_EXTERN char * \ -FFTW_CDECL X(export_wisdom_to_string)(void); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(export_wisdom)(X(write_char_func) write_char, \ - void *data); \ -FFTW_EXTERN int \ -FFTW_CDECL X(import_system_wisdom)(void); \ - \ -FFTW_EXTERN int \ -FFTW_CDECL X(import_wisdom_from_filename)(const char *filename); \ - \ -FFTW_EXTERN int \ -FFTW_CDECL X(import_wisdom_from_file)(FILE *input_file); \ - \ -FFTW_EXTERN int \ -FFTW_CDECL X(import_wisdom_from_string)(const char *input_string); \ - \ -FFTW_EXTERN int \ -FFTW_CDECL X(import_wisdom)(X(read_char_func) read_char, void *data); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(fprint_plan)(const X(plan) p, FILE *output_file); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(print_plan)(const X(plan) p); \ - \ -FFTW_EXTERN char * \ -FFTW_CDECL X(sprint_plan)(const X(plan) p); \ - \ -FFTW_EXTERN void * \ -FFTW_CDECL X(malloc)(size_t n); \ - \ -FFTW_EXTERN R * \ -FFTW_CDECL X(alloc_real)(size_t n); \ -FFTW_EXTERN C * \ -FFTW_CDECL X(alloc_complex)(size_t n); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(free)(void *p); \ - \ -FFTW_EXTERN void \ -FFTW_CDECL X(flops)(const X(plan) p, \ - double *add, double *mul, double *fmas); \ -FFTW_EXTERN double \ -FFTW_CDECL X(estimate_cost)(const X(plan) p); \ - \ -FFTW_EXTERN double \ -FFTW_CDECL X(cost)(const X(plan) p); \ - \ -FFTW_EXTERN int \ -FFTW_CDECL X(alignment_of)(R *p); \ - \ -FFTW_EXTERN const char X(version)[]; \ -FFTW_EXTERN const char X(cc)[]; \ -FFTW_EXTERN const char X(codelet_optim)[]; - - -/* end of FFTW_DEFINE_API macro */ - -FFTW_DEFINE_API(FFTW_MANGLE_DOUBLE, double, fftw_complex) -FFTW_DEFINE_API(FFTW_MANGLE_FLOAT, float, fftwf_complex) -FFTW_DEFINE_API(FFTW_MANGLE_LONG_DOUBLE, long double, fftwl_complex) - -/* __float128 (quad precision) is a gcc extension on i386, x86_64, and ia64 - for gcc >= 4.6 (compiled in FFTW with --enable-quad-precision) */ -#if (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)) \ - && !(defined(__ICC) || defined(__INTEL_COMPILER) || defined(__CUDACC__) || defined(__PGI)) \ - && (defined(__i386__) || defined(__x86_64__) || defined(__ia64__)) -# if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) -/* note: __float128 is a typedef, which is not supported with the _Complex - keyword in gcc, so instead we use this ugly __attribute__ version. - However, we can't simply pass the __attribute__ version to - FFTW_DEFINE_API because the __attribute__ confuses gcc in pointer - types. Hence redefining FFTW_DEFINE_COMPLEX. Ugh. */ -# undef FFTW_DEFINE_COMPLEX -# define FFTW_DEFINE_COMPLEX(R, C) typedef _Complex float __attribute__((mode(TC))) C -# endif -FFTW_DEFINE_API(FFTW_MANGLE_QUAD, __float128, fftwq_complex) -#endif - -#define FFTW_FORWARD (-1) -#define FFTW_BACKWARD (+1) - -#define FFTW_NO_TIMELIMIT (-1.0) - -/* documented flags */ -#define FFTW_MEASURE (0U) -#define FFTW_DESTROY_INPUT (1U << 0) -#define FFTW_UNALIGNED (1U << 1) -#define FFTW_CONSERVE_MEMORY (1U << 2) -#define FFTW_EXHAUSTIVE (1U << 3) /* NO_EXHAUSTIVE is default */ -#define FFTW_PRESERVE_INPUT (1U << 4) /* cancels FFTW_DESTROY_INPUT */ -#define FFTW_PATIENT (1U << 5) /* IMPATIENT is default */ -#define FFTW_ESTIMATE (1U << 6) -#define FFTW_WISDOM_ONLY (1U << 21) - -/* undocumented beyond-guru flags */ -#define FFTW_ESTIMATE_PATIENT (1U << 7) -#define FFTW_BELIEVE_PCOST (1U << 8) -#define FFTW_NO_DFT_R2HC (1U << 9) -#define FFTW_NO_NONTHREADED (1U << 10) -#define FFTW_NO_BUFFERING (1U << 11) -#define FFTW_NO_INDIRECT_OP (1U << 12) -#define FFTW_ALLOW_LARGE_GENERIC (1U << 13) /* NO_LARGE_GENERIC is default */ -#define FFTW_NO_RANK_SPLITS (1U << 14) -#define FFTW_NO_VRANK_SPLITS (1U << 15) -#define FFTW_NO_VRECURSE (1U << 16) -#define FFTW_NO_SIMD (1U << 17) -#define FFTW_NO_SLOW (1U << 18) -#define FFTW_NO_FIXED_RADIX_LARGE_N (1U << 19) -#define FFTW_ALLOW_PRUNING (1U << 20) - -#ifdef __cplusplus -} /* extern "C" */ -#endif /* __cplusplus */ - -#endif /* FFTW3_H */ diff --git a/instruction.txt b/instruction.txt deleted file mode 100644 index 5fcc182..0000000 --- a/instruction.txt +++ /dev/null @@ -1,7 +0,0 @@ -Сборка FFTW для LINUX -Для того чтобы разрешить позиционно-независимый код необходимо FFTW3 конфигурировать с флагом --enable-shared - -./configure --enable-shared -make - -После этого в скрытой дириктории .lib появится so файл и a файл. Статический файл a можно ликновать с so и подключать в so библиотеку. diff --git a/ext/linux32/.gitignore b/release/doc/.gitignore similarity index 100% rename from ext/linux32/.gitignore rename to release/doc/.gitignore diff --git a/release/include/dspl.c b/release/include/dspl.c index 40cdabf..18fdd46 100644 --- a/release/include/dspl.c +++ b/release/include/dspl.c @@ -20,7 +20,6 @@ - #ifdef WIN_OS #include #endif //WIN_OS @@ -35,82 +34,82 @@ #ifndef BUILD_LIB -p_blas_dscal blas_dscal ; - -p_cheby_poly1 cheby_poly1 ; -p_cheby_poly2 cheby_poly2 ; +p_cheby_poly1 cheby_poly1 ; +p_cheby_poly2 cheby_poly2 ; +p_conv conv ; +p_conv_cmplx conv_cmplx ; p_dft dft ; p_dft_cmplx dft_cmplx ; - -p_fft_create fft_create ; -p_fft_destroy fft_destroy ; -p_fft_cmplx fft_cmplx ; -p_fft_shift fft_shift ; - - +p_filter_iir filter_iir ; +p_linspace linspace ; +p_logspace logspace ; #endif //BUILD_LIB +#ifdef WIN_OS +#define LOAD_FUNC(fn) \ + fname = #fn;\ + fn = (p_##fn)GetProcAddress(handle, fname);\ + if(! fn) goto exit_label; +#endif + + + +#ifdef LINUX_OS +#define LOAD_FUNC(fn) \ + fname = #fn;\ + fn = (p_##fn)dlsym(handle, fname);\ + if ((error = dlerror()) != NULL) goto exit_label +#endif + void* dspl_load() { + char* fname; #ifdef WIN_OS - HINSTANCE handle; - char* fname; + HINSTANCE handle; handle = LoadLibrary(TEXT("libdspl.dll")); if (!handle) { printf("libdspl.dll loading ERROR!\n"); return NULL; } - - - fname = "blas_dscal"; - blas_dscal = (p_blas_dscal)GetProcAddress(handle, fname); - if(!blas_dscal) goto exit_label; + #endif //WIN_OS - fname = "cheby_poly1"; - cheby_poly1 = (p_cheby_poly1)GetProcAddress(handle, fname); - if(!cheby_poly1) goto exit_label; + #ifdef LINUX_OS + char* error; + void *handle; + // open the *.so + handle = dlopen ("./libdspl.so", RTLD_LAZY); + if (!handle) + { + printf("libdspl.so loading ERROR!\n"); + return NULL; + } + #endif //LINUX_OS - fname = "cheby_poly2"; - cheby_poly2 = (p_cheby_poly2)GetProcAddress(handle, fname); - if(!cheby_poly2) goto exit_label; - - fname = "dft"; - dft = (p_dft)GetProcAddress(handle, fname); - if(!dft) goto exit_label; - - fname = "dft_cmplx"; - dft_cmplx = (p_dft_cmplx)GetProcAddress(handle, fname); - if(!dft_cmplx) goto exit_label; - - - fname = "fft_create"; - fft_create = (p_fft_create)GetProcAddress(handle, fname); - if(!fft_create) goto exit_label; - - fname = "fft_destroy"; - fft_destroy = (p_fft_destroy)GetProcAddress(handle, fname); - if(!fft_destroy) goto exit_label; - - - fname = "fft_cmplx"; - fft_cmplx = (p_fft_cmplx)GetProcAddress(handle, fname); - if(!fft_cmplx) goto exit_label; - - fname = "fft_shift"; - fft_shift = (p_fft_shift)GetProcAddress(handle, fname); - if(!fft_shift) goto exit_label; - - - return (void*)handle; + + + LOAD_FUNC(cheby_poly1); + LOAD_FUNC(cheby_poly2); + LOAD_FUNC(conv); + LOAD_FUNC(conv_cmplx); + LOAD_FUNC(dft); + LOAD_FUNC(dft_cmplx); + LOAD_FUNC(filter_iir); + LOAD_FUNC(linspace); + LOAD_FUNC(logspace); + + + + #ifdef WIN_OS + return (void*)handle; exit_label: printf("function %s loading ERROR!\n", fname); if(handle) @@ -119,75 +118,14 @@ void* dspl_load() #endif //WIN_OS - - - - - #ifdef LINUX_OS - char* error; - void *handle; - char* fname; - - // open the *.so - handle = dlopen ("./libdspl.so", RTLD_LAZY); - if (!handle) - { - printf("libdspl.so loading ERROR!\n"); - return NULL; - } - - - fname = "blas_dscal"; - blas_dscal = (p_blas_dscal)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - fname = "cheby_poly1"; - cheby_poly1 = (p_cheby_poly1)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - fname = "cheby_poly2"; - cheby_poly2 = (p_cheby_poly2)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - fname = "dft"; - dft = (p_dft)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - fname = "dft_cmplx"; - dft_cmplx = (p_dft_cmplx)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - fname = "fft_create"; - fft_create = (p_fft_create)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - fname = "fft_destroy"; - fft_destroy = (p_fft_destroy)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - fname = "fft_cmplx"; - fft_cmplx = (p_fft_cmplx)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - fname = "fft_shift"; - fft_shift = (p_fft_shift)dlsym(handle, fname); - if ((error = dlerror()) != NULL) goto exit_label; - - - return handle; - + #ifdef LINUX_OS + return handle; exit_label: printf("function %s loading ERROR!\n", fname); if(handle) dlclose(handle); return NULL; - #endif //LINUX_OS - + #endif //LINUX_OS } diff --git a/release/include/dspl.h b/release/include/dspl.h index 3703878..5e87762 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -162,66 +162,45 @@ extern "C" { #endif -// Declare DSPL_API for Windows OS -#ifdef BUILD_LIB +#ifdef BUILD_LIB +// Declare DSPL_API for Windows OS #ifdef WIN_OS #define DSPL_API __declspec(dllexport) #endif // WIN_OS - +// Declare DSPL_API for LINUX OS #ifdef LINUX_OS #define DSPL_API #endif //LINUX_OS - - -int DSPL_API blas_dscal(int n, double a, double* x, int incx); - -int DSPL_API cheby_poly1(double* x, int n, int ord, double* y); -int DSPL_API cheby_poly2(double* x, int n, int ord, double* y); - -int DSPL_API dft (double* x, int n, complex_t *y); -int DSPL_API dft_cmplx (complex_t* x, int n, complex_t *y); - - -int DSPL_API fft_create ( fft_t *pfft, int n); -void DSPL_API fft_destroy (fft_t *pfft); -int DSPL_API fft_cmplx (complex_t *x, int n, fft_t* pfft, complex_t* y); -int DSPL_API fft_shift (double* x, int n, double* y); - - -#else //BUILD_LIB - - - -typedef int (*p_blas_dscal ) (int n, double a, double* x, int incx); - -typedef int (*p_cheby_poly1 ) (double* x, int n, int ord, double* y); -typedef int (*p_cheby_poly2 ) (double* x, int n, int ord, double* y); - -typedef int (*p_dft ) (double* x, int n, complex_t *y); -typedef int (*p_dft_cmplx ) (complex_t* x, int n, complex_t *y); -typedef int (*p_fft_create ) (fft_t *pfft, int n); -typedef void(*p_fft_destroy ) (fft_t *pfft); -typedef int (*p_fft_cmplx ) (complex_t *x, int n, fft_t* pfft, complex_t* y); -typedef int (*p_fft_shift ) (double* x, int n, double* y); - - - -extern p_blas_dscal blas_dscal ; - -extern p_cheby_poly1 cheby_poly1 ; -extern p_cheby_poly2 cheby_poly2 ; - -extern p_dft dft ; -extern p_dft_cmplx dft_cmplx ; - -extern p_fft_create fft_create ; -extern p_fft_destroy fft_destroy ; -extern p_fft_cmplx fft_cmplx ; -extern p_fft_shift fft_shift ; - #endif //BUILD_DLL +#define COMMA , + + +#ifdef BUILD_LIB +#define DECLARE_FUNC(type, fn, param)\ + type DSPL_API fn(param);\ + +#endif + +#ifndef BUILD_LIB +#define DECLARE_FUNC(type, fn, param)\ + typedef type (*p_##fn)(param);\ + extern p_##fn fn; + +#endif + + +DECLARE_FUNC(int, cheby_poly1, double* COMMA int COMMA int COMMA double*); +DECLARE_FUNC(int, cheby_poly2, double* COMMA int COMMA int COMMA double*); +DECLARE_FUNC(int, conv, double* COMMA int COMMA double* COMMA int COMMA double*); +DECLARE_FUNC(int, conv_cmplx, complex_t* COMMA int COMMA complex_t* COMMA int COMMA complex_t*); +DECLARE_FUNC(int, dft, double* COMMA int COMMA complex_t*); +DECLARE_FUNC(int, dft_cmplx, complex_t* COMMA int COMMA complex_t*); +DECLARE_FUNC(int, filter_iir, double* COMMA double* COMMA int COMMA double* COMMA int COMMA double*); +DECLARE_FUNC(int, linspace, double COMMA double COMMA int COMMA int COMMA double*); +DECLARE_FUNC(int, logspace, double COMMA double COMMA int COMMA int COMMA double*); + #ifdef __cplusplus } diff --git a/ext/linux64/.gitignore b/release/lib/.gitignore similarity index 100% rename from ext/linux64/.gitignore rename to release/lib/.gitignore diff --git a/release/lib/linux32/.gitignore b/release/lib/linux32/.gitignore deleted file mode 100644 index e69de29..0000000 diff --git a/release/lib/linux64/.gitignore b/release/lib/linux64/.gitignore deleted file mode 100644 index e69de29..0000000 diff --git a/release/lib/win32/.gitignore b/release/lib/win32/.gitignore deleted file mode 100644 index e69de29..0000000 diff --git a/release/lib/win64/.gitignore b/release/lib/win64/.gitignore deleted file mode 100644 index e69de29..0000000 diff --git a/test/bin/blas_test.exe b/test/bin/blas_test.exe deleted file mode 100755 index 64e7489..0000000 Binary files a/test/bin/blas_test.exe and /dev/null differ diff --git a/test/bin/dft_test.exe b/test/bin/dft_test.exe index 7d206b6..2b3c751 100755 Binary files a/test/bin/dft_test.exe and b/test/bin/dft_test.exe differ diff --git a/test/src/blas_test.c b/test/src/blas_test.c deleted file mode 100644 index bb4fb1b..0000000 --- a/test/src/blas_test.c +++ /dev/null @@ -1,28 +0,0 @@ -#include -#include -#include -#include "dspl.h" - -#define N 16 -int main() -{ - void* handle; - handle = dspl_load(); - - double x[N]; - double a = 2.0; - - for(int k = 0; k < N; k++) - x[k] = (double)k; - - blas_dscal(N, a, x, 1); - - for(int k = 0; k < N; k++) - printf("x[%2d] = %9.3f \n", k, x[k]); - - // remember to free the resource - dspl_free(handle); - return 0; -} - - diff --git a/test/src/dft_test.c b/test/src/dft_test.c index 1053642..f6feb10 100644 --- a/test/src/dft_test.c +++ b/test/src/dft_test.c @@ -11,7 +11,7 @@ int main() complex_t x[N]; complex_t y[N]; - complex_t z[N]; + //complex_t z[N]; for(int k = 0; k < N; k++) @@ -22,17 +22,9 @@ int main() dft_cmplx(x,N,y); - fft_t pfft; - memset(&pfft, 0, sizeof(fft_t)); - // - fft_create(&pfft,N); - fft_cmplx(x, N, &pfft, z); - for(int k = 0; k < N; k++) - printf("y[%2d] = %9.3f%9.3f z[%2d] = %9.3f%9.3f \n", k, RE(y[k]), IM(y[k]), k, RE(z[k]), IM(z[k])); + printf("y[%2d] = %9.3f%9.3f\n", k, RE(y[k]), IM(y[k])); - - fft_destroy(&pfft); // remember to free the resource dspl_free(handle); return 0;