mirror
/
libdspl-2.0
kopia lustrzana https://github.com/Dsplib/libdspl-2.0
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

changed project structure. Added macro for autoloading

pull/2/head
Dsplib 2018-03-13 23:18:11 +03:00
rodzic b8afe65f96
commit 4f60fc9c4e
46 zmienionych plików z 197 dodań i 4059 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

7
.gitignore vendored
Wyświetl plik

@ -1,3 +1,8 @@
*.o
*.so
*.dll
*.dll
*.exe
*.txt
*.dat
*.bin
*.csv

24
Makefile.dirs
Wyświetl plik

@ -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

4
Makefile.dspl
Wyświetl plik

@ -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

100
dspl/src/array.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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;
}

58
dspl/src/blas_lev1.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#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;
}

4
dspl/src/cheby.c
Wyświetl plik

@ -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;
}

236
dspl/src/complex.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdio.h>
#include <stdlib.h>
#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;
}

122
dspl/src/composing.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <string.h>
#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;
}

14
dspl/src/conv.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
@ -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;

494
dspl/src/ellipj.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#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;
}

167
dspl/src/fft.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#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;
}

12
dspl/src/fillarray.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
@ -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;

171
dspl/src/filter_an.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#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;
}

332
dspl/src/filter_ap.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#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;
}

100
dspl/src/fourier_series.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#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;
}

110
dspl/src/goertzel.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <string.h>
#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;
}

129
dspl/src/inout.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdio.h>
#include <stdlib.h>
#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;
}

33
dspl/src/math.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdio.h>
#include <stdlib.h>
#include "dspl.h"
double dmod (double x, double y)
{
if(y == 0.0)
return x;
return x - floor(x/y) * y;
}

136
dspl/src/polyval.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#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;
}

100
dspl/src/randgen.c
Wyświetl plik

@ -1,100 +0,0 @@
#include <stdlib.h>
#include <math.h>
#include <time.h>
#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<n; k++)
{
x1[0] = (63308 * x1[2] - 183326*x1[3]) % DSPL_RAND_MOD_X1;
x2[0] = (86098 * x2[1] - 539608*x2[3]) % DSPL_RAND_MOD_X2;
y = (x1[0] - x2[0]) % DSPL_RAND_MOD_X1;
for(m = 3; m > 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<n && m < 128)
{
x[k] = sqrt(-2.0*log(x1[m]))*cos(M_2PI*x2[m])*sigma + mu;
k++;
m++;
if(k<n && m < 128)
{
x[k] = sqrt(-2.0*log(x1[m]))*sin(M_2PI*x2[m])*sigma + mu;
k++;
m++;
}
}
}
res = RES_OK;
exit_label:
return res;
}

186
dspl/src/resampling.c
Wyświetl plik

@ -1,186 +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 <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#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;
}

56
dspl/src/signals.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#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;
}

70
dspl/src/trapint.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdio.h>
#include <stdlib.h>
#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;
}

520
dspl/src/win.c
Wyświetl plik

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <stdio.h>
#include <math.h>
#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<n; i++)
{
w[i] = 0.62 - 0.48*fabs(y-0.5)-0.38*cos(M_2PI*y);
y += x;
}
return RES_OK;
}
/**************************************************************************************************
Blackman window function
***************************************************************************************************/
int win_blackman(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 = M_2PI/(double)(n-1); break;
case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
default: return ERROR_WIN_SYM;
}
y = 0.0;
for(i = 0; i<n; i++)
{
w[i] = 0.42 - 0.5* cos(y)+0.08*cos(2.0*y);
y += x;
}
return RES_OK;
}
/**************************************************************************************************
Blackman - Harris window function
***************************************************************************************************/
int win_blackman_harris(double *w, int n, int win_type)
{
double y;
double x = 0.0;
double a0 = 0.35875;
double a1 = 0.48829;
double a2 = 0.14128;
double a3 = 0.01168;
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 = M_2PI/(double)(n-1); break;
case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
default: return ERROR_WIN_SYM;
}
y = 0.0;
for(i = 0; i<n; i++)
{
w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
y += x;
}
return RES_OK;
}
/**************************************************************************************************
Blackman - Nuttull window function
***************************************************************************************************/
int win_blackman_nuttall(double *w, int n, int win_type)
{
double y;
double x = 0.0;
double a0 = 0.3635819;
double a1 = 0.4891775;
double a2 = 0.1365995;
double a3 = 0.0106411;
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 = M_2PI/(double)(n-1); break;
case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
default: return ERROR_WIN_SYM;
}
y = 0.0;
for(i = 0; i<n; i++)
{
w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
y += x;
}
return RES_OK;
}
/**************************************************************************************************
Cosine window function
***************************************************************************************************/
int win_cos(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 = M_PI/(double)(n-1); break;
case DSPL_WIN_PERIODIC : x = M_PI/(double)n; break;
default: return ERROR_WIN_SYM;
}
y = 0.0;
for(i = 0; i<n; i++)
{
w[i] = sin(y);
y += x;
}
return RES_OK;
}
/**************************************************************************************************
Flat - Top window function
***************************************************************************************************/
int win_flat_top(double *w, int n, int win_type)
{
double y;
double x = 0.0;
double a0 = 1.0;
double a1 = 1.93;
double a2 = 1.29;
double a3 = 0.388;
double a4 = 0.032;
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 = M_2PI/(double)(n-1); break;
case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
default: return ERROR_WIN_SYM;
}
y = 0.0;
for(i = 0; i<n; i++)
{
w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y)+a4*cos(4.0*y);
y += x;
}
return RES_OK;
}
/**************************************************************************************************
Gaussian window function
***************************************************************************************************/
int win_gaussian(double *w, int n, int win_type, double alpha)
{
double x = 0.0;
double y;
double sigma;
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)*0.5; break;
case DSPL_WIN_PERIODIC : x = (double)(n)*0.5; break;
default: return ERROR_WIN_SYM;
}
sigma = alpha / x;
for(i = 0; i<n; i++)
{
y = ((double)i - x)*sigma;
w[i] = exp(-0.5*y*y);
}
return RES_OK;
}
/**************************************************************************************************
Hamming window function
***************************************************************************************************/
int win_hamming(double *w, int n, int win_type)
{
double x = 0.0;
double y;
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 = M_2PI/(double)(n-1); break;
case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
default: return ERROR_WIN_SYM;
}
y = 0.0;
for(i = 0; i < n; i++)
{
w[i] = 0.54-0.46*cos(y);
y += x;
}
return RES_OK;
}
/**************************************************************************************************
Hann window function
***************************************************************************************************/
int win_hann(double *w, int n, int win_type)
{
double x = 0.0;
double y;
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 = M_2PI/(double)(n-1); break;
case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
default: return ERROR_WIN_SYM;
}
y = 0.0;
for(i = 0; i < n; i++)
{
w[i] = 0.5*(1-cos(y));
y += x;
}
return RES_OK;
}
/**************************************************************************************************
Lanczos window function
***************************************************************************************************/
int win_lanczos(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 = M_2PI/(double)(n-1); break;
case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
default: return ERROR_WIN_SYM;
}
y = 0.0;
for(i = 0; i < n; i++)
{
if((y - M_PI)==0.0)
w[i] = 1.0;
else
w[i] = sin(y - M_PI)/(y - M_PI);
y += x;
}
return RES_OK;
}
/**************************************************************************************************
Nuttall window function
***************************************************************************************************/
int win_nuttall(double *w, int n, int win_type)
{
double y;
double x = 0.0;
double a0 = 0.355768;
double a1 = 0.487396;
double a2 = 0.144232;
double a3 = 0.012604;
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 = M_2PI/(double)(n-1); break;
case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
default: return ERROR_WIN_SYM;
}
y = 0.0;
for(i = 0; i < n; i++)
{
w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
y += x;
}
return RES_OK;
}
/**************************************************************************************************
Rectangle window function
***************************************************************************************************/
int win_rect(double *w, int n)
{
int i;
if(!w)
return ERROR_PTR;
if(n<2)
return ERROR_SIZE;
for(i = 0; i < n; i++)
w[i] = 1.0;
return RES_OK;
}

BIN
ext/linux64/libblas.a
Wyświetl plik

Plik binarny nie jest wyświetlany.

BIN
ext/linux64/libfftw3.a
Wyświetl plik

Plik binarny nie jest wyświetlany.

BIN
ext/linux64/liblapack.a
Wyświetl plik

Plik binarny nie jest wyświetlany.

0
ext/win32/.gitignore vendored
Wyświetl plik

0
ext/win64/.gitignore vendored
Wyświetl plik

BIN
ext/win64/libfftw3.a
Wyświetl plik

Plik binarny nie jest wyświetlany.

174
include/dspl.c
Wyświetl plik

@ -20,7 +20,6 @@
#ifdef WIN_OS
#include <windows.h>
#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
}

81
include/dspl.h
Wyświetl plik

@ -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
}

514
include/fftw3.h
Wyświetl plik

@ -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 <stdio.h>
#ifdef __cplusplus
extern "C"
{
#endif /* __cplusplus */
/* If <complex.h> 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 <stddef.h> /* 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 */

7
instruction.txt
Wyświetl plik

@ -1,7 +0,0 @@
Сборка FFTW для LINUX
Для того чтобы разрешить позиционно-независимый код необходимо FFTW3 конфигурировать с флагом --enable-shared
./configure --enable-shared
make
После этого в скрытой дириктории .lib появится so файл и a файл. Статический файл a можно ликновать с so и подключать в so библиотеку.

0
ext/linux32/.gitignore → release/doc/.gitignore vendored
Wyświetl plik

174
release/include/dspl.c
Wyświetl plik

@ -20,7 +20,6 @@
#ifdef WIN_OS
#include <windows.h>
#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
}

81
release/include/dspl.h
Wyświetl plik

@ -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
}

0
ext/linux64/.gitignore → release/lib/.gitignore vendored
Wyświetl plik

0
release/lib/linux32/.gitignore vendored
Wyświetl plik

0
release/lib/linux64/.gitignore vendored
Wyświetl plik

0
release/lib/win32/.gitignore vendored
Wyświetl plik

0
release/lib/win64/.gitignore vendored
Wyświetl plik

BIN
test/bin/blas_test.exe
Wyświetl plik

Plik binarny nie jest wyświetlany.

BIN
test/bin/dft_test.exe
Wyświetl plik

Plik binarny nie jest wyświetlany.

28
test/src/blas_test.c
Wyświetl plik

@ -1,28 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#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;
}

12
test/src/dft_test.c
Wyświetl plik

@ -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;