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

added some functions

pull/2/head
Dsplib 2018-03-16 00:01:34 +03:00
rodzic 3daecefae2
commit 26c3c73692
17 zmienionych plików z 1256 dodań i 73 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

0
dspl/dox/ru/fourier_series.dox 100644
Wyświetl plik

21
dspl/dox/ru/polyval.dox
Wyświetl plik

@ -60,37 +60,24 @@
P_N(x) = a_0 + x \cdot (a_1 + x \cdot (a_2 + \cdot ( \ldots x \cdot (a_{N-1} + x\cdot a_N) \ldots )))
\f]
\param[in] aR Указатель на вектор реальной части коэффициентов полинома. <BR>
\param[in] a Указатель на вектор комплексных коэффициентов полинома. <BR>
Размер вектора `[ord+1 x 1]`.<BR>
Коэффициент `aR[0]` соответствует коэффициенту полинома \f$a_0\f$.<BR><BR>
Коэффициент `a[0]` соответствует коэффициенту полинома \f$a_0\f$.<BR><BR>
\param[in] aI Указатель на вектор мнимой части коэффициентов полинома. <BR>
Размер вектора `[ord+1 x 1]`.<BR>
Коэффициент `aI[0]` соответствует коэффициенту полинома \f$a_0\f$.<BR><BR>
\param[in] ord Порядок полинома \f$N\f$. <BR><BR>
\param[in] xR Указатель на вектор реальной части аргумента полинома. <BR>
Размер вектора `[n x 1]`.<BR>
Значения полинома будут расчитаны для всех
значений аргумента вектора `x`.<BR><BR>
\param[in] xI Указатель на вектор мнимой части аргумента полинома. <BR>
\param[in] x Указатель на вектор аргумента полинома. <BR>
Размер вектора `[n x 1]`.<BR>
Значения полинома будут расчитаны для всех
значений аргумента вектора `x`.<BR><BR>
\param[in] n Размер вектора агрумента полинома. <BR><BR>
\param[out] yR Указатель вектор реальной части значения
полинома для аргумента `x`. <BR>
\param[out] y Указатель вектор значения полинома для аргумента `x`. <BR>
Размер вектора `[n x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[out] yI Указатель вектор реальной части значения
полинома для аргумента `x`. <BR>
Размер вектора `[n x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\return
`RES_OK` Полином расчитан успешно. <BR>

63
dspl/dox/ru/randgen.dox
Wyświetl plik

@ -1,16 +1,59 @@
/**************************************************************************************************
Uniform random numbers generator
int randu(double* x, int n)
***************************************************************************************************/
/*! *************************************************************************************************
\ingroup SPEC_MATH_RAND_GEN_GROUP
\fn int randn(double* x, int n, double mu, double sigma)
\brief Генерация вектора нормально распределенных псевдослучайных чисел с
заданным математическим ожиданием и среднеквадратическим отклонением.
Генерация производится при помощи преобразования Бокса — Мюллера равномерно-распределенной
случайной величины в нормально распределенную.<BR>
\param[in,out] x Указатель на вектор нормальной распределенных случайных чисел. <BR>
Размер вектора `[n x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[in] n Размер вектора случайных чисел.<BR><BR>
\param[in] mu Математическое ожидание.<BR><BR>
\param[in] sigma Среднеквадратическое отклонение (СКО).<BR><BR>
\return
`RES_OK` Вектор случайных чисел сгенерирован успешно.<BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки".
\author
Бахурин Сергей.
www.dsplib.org
**************************************************************************************************** */
/*! *************************************************************************************************
\ingroup SPEC_MATH_RAND_GEN_GROUP
\fn int randu(double* x, int n);
\brief Генерация вектора равномерно-распределенных в интервале от 0 до 1 псевдослучайных чисел.
Генерация производится при помощи рекурсивного алгоритма L'Ecluyer. Период датчика порядка \f$10^56\f$.<BR>
\param[in,out] x Указатель на вектор случайных чисел. <BR>
Размер вектора `[n x 1]`.<BR>
Память должна быть выделена.<BR><BR>
\param[in] n Размер вектора случайных чисел.<BR><BR>
\return
`RES_OK` Вектор случайных чисел сгенерирован успешно.<BR>
В противном случае \ref ERROR_CODE_GROUP "код ошибки".
\author
Бахурин Сергей.
www.dsplib.org
**************************************************************************************************** */
/**************************************************************************************************
Gaussian random numbers generator
int randn(double* x, int n, double mu, double sigma)
***************************************************************************************************/

101
dspl/src/array.c 100644
Wyświetl plik

@ -0,0 +1,101 @@
/*
* 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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
/**************************************************************************************************
Concntenate arrays
***************************************************************************************************/
int DSPL_API 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 DSPL_API 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 DSPL_API 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;
}

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

@ -0,0 +1,236 @@
/*
* 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 DSPL_API 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 DSPL_API 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 DSPL_API 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 DSPL_API 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 DSPL_API 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 DSPL_API 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 DSPL_API 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 DSPL_API 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;
}

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

@ -0,0 +1,100 @@
/*
* 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"
int DSPL_API 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 DSPL_API 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;
}

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

@ -0,0 +1,129 @@
/*
* 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 DSPL_API 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 DSPL_API 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 100644
Wyświetl plik

@ -0,0 +1,33 @@
/*
* 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 DSPL_API dmod (double x, double y)
{
if(y == 0.0)
return x;
return x - floor(x/y) * y;
}

119
dspl/src/randgen.c 100644
Wyświetl plik

@ -0,0 +1,119 @@
/*
* 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 <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 DSPL_API 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 DSPL_API 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;
}

187
dspl/src/resampling.c 100644
Wyświetl plik

@ -0,0 +1,187 @@
/*
* 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"
#define DSPL_FARROW_LAGRANGE_COEFF 0.16666666666666666666666666666667
/**************************************************************************************************
Farrow resampler based on the cubic Lagrange polynomials
***************************************************************************************************/
int DSPL_API 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 DSPL_API 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;
}

60
dspl/src/signals.c 100644
Wyświetl plik

@ -0,0 +1,60 @@
/*
* 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"
/**************************************************************************************************
Rectangle pulse signal
***************************************************************************************************/
int DSPL_API 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 100644
Wyświetl plik

@ -0,0 +1,70 @@
/*
* 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 DSPL_API 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 DSPL_API 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;
}

77
include/dspl.c
Wyświetl plik

@ -35,20 +35,42 @@
#ifndef BUILD_LIB
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_filter_iir filter_iir ;
p_goertzel goertzel ;
p_goertzel_cmplx goertzel_cmplx ;
p_linspace linspace ;
p_logspace logspace ;
p_polyval polyval ;
p_polyval_cmplx polyval_cmplx ;
p_acos_cmplx acos_cmplx ;
p_asin_cmplx asin_cmplx ;
p_cheby_poly1 cheby_poly1 ;
p_cheby_poly2 cheby_poly2 ;
p_cmplx2re cmplx2re ;
p_concat concat ;
p_conv conv ;
p_conv_cmplx conv_cmplx ;
p_cos_cmplx cos_cmplx ;
p_dft dft ;
p_dft_cmplx dft_cmplx ;
p_dmod dmod ;
p_farrow_lagrange farrow_lagrange ;
p_farrow_spline farrow_spline ;
p_filter_iir filter_iir ;
p_flipip flipip ;
p_flipip_cmplx flipip_cmplx ;
p_fourier_series_dec fourier_series_dec ;
p_fourier_series_rec fourier_series_rec ;
p_goertzel goertzel ;
p_goertzel_cmplx goertzel_cmplx ;
p_linspace linspace ;
p_log_cmplx log_cmplx ;
p_logspace logspace ;
p_polyval polyval ;
p_polyval_cmplx polyval_cmplx ;
p_randn randn ;
p_randu randu ;
p_re2cmplx re2cmplx ;
p_signal_pimp signal_pimp ;
p_sin_cmplx sin_cmplx ;
p_sqrt_cmplx sqrt_cmplx ;
p_trapint trapint ;
p_trapint_cmplx trapint_cmplx ;
p_writebin writebin ;
p_writetxt writetxt ;
#endif //BUILD_LIB
@ -100,23 +122,42 @@ void* dspl_load()
}
#endif //LINUX_OS
LOAD_FUNC(acos_cmplx);
LOAD_FUNC(asin_cmplx);
LOAD_FUNC(cheby_poly1);
LOAD_FUNC(cheby_poly2);
LOAD_FUNC(cmplx2re);
LOAD_FUNC(concat);
LOAD_FUNC(conv);
LOAD_FUNC(conv_cmplx);
LOAD_FUNC(cos_cmplx);
LOAD_FUNC(dft);
LOAD_FUNC(dft_cmplx);
LOAD_FUNC(dmod);
LOAD_FUNC(farrow_lagrange);
LOAD_FUNC(farrow_spline);
LOAD_FUNC(filter_iir);
LOAD_FUNC(flipip);
LOAD_FUNC(flipip_cmplx);
LOAD_FUNC(fourier_series_dec);
LOAD_FUNC(fourier_series_rec);
LOAD_FUNC(goertzel);
LOAD_FUNC(goertzel_cmplx);
LOAD_FUNC(linspace);
LOAD_FUNC(log_cmplx);
LOAD_FUNC(logspace);
LOAD_FUNC(polyval);
LOAD_FUNC(polyval_cmplx);
LOAD_FUNC(randn);
LOAD_FUNC(randu);
LOAD_FUNC(re2cmplx);
LOAD_FUNC(signal_pimp);
LOAD_FUNC(sin_cmplx);
LOAD_FUNC(sqrt_cmplx);
LOAD_FUNC(trapint);
LOAD_FUNC(trapint_cmplx);
LOAD_FUNC(writebin);
LOAD_FUNC(writetxt);
#ifdef WIN_OS

28
include/dspl.h
Wyświetl plik

@ -191,24 +191,42 @@ extern "C" {
#endif
DECLARE_FUNC(int, acos_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, asin_cmplx, complex_t* COMMA int COMMA complex_t*);
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, cmplx2re, complex_t* COMMA int COMMA double* COMMA double*);
DECLARE_FUNC(int, concat, void* COMMA size_t COMMA void* COMMA size_t COMMA void*);
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, cos_cmplx, 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(double,dmod, double COMMA double);
DECLARE_FUNC(int, farrow_lagrange, double* COMMA int COMMA double COMMA double COMMA double COMMA double** COMMA int*);
DECLARE_FUNC(int, farrow_spline, double* COMMA int COMMA double COMMA double COMMA double COMMA double** COMMA int*);
DECLARE_FUNC(int, filter_iir, double* COMMA double* COMMA int COMMA double* COMMA int COMMA double*);
DECLARE_FUNC(int, flipip, double* COMMA int);
DECLARE_FUNC(int, flipip_cmplx, complex_t* COMMA int);
DECLARE_FUNC(int, fourier_series_dec, double* COMMA double* COMMA int COMMA double COMMA int COMMA double* COMMA complex_t*);
DECLARE_FUNC(int, fourier_series_rec, double* COMMA complex_t* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel, double* COMMA int COMMA int* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel_cmplx, complex_t* COMMA int COMMA int* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, linspace, double COMMA double COMMA int COMMA int COMMA double*);
DECLARE_FUNC(int, log_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, logspace, double COMMA double COMMA int COMMA int COMMA double*);
DECLARE_FUNC(int, polyval, double* COMMA int COMMA double* COMMA int COMMA double*);
DECLARE_FUNC(int, polyval_cmplx, complex_t* COMMA int COMMA complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, randn, double* COMMA int COMMA double COMMA double);
DECLARE_FUNC(int, randu, double* COMMA int);
DECLARE_FUNC(int, re2cmplx, double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, signal_pimp, double* COMMA size_t COMMA double COMMA double COMMA double COMMA double COMMA double*);
DECLARE_FUNC(int, sin_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, sqrt_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, trapint, double* COMMA double* COMMA int COMMA double* sum);
DECLARE_FUNC(int, trapint_cmplx, double* COMMA complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, writebin, void* COMMA int COMMA int COMMA char*);
DECLARE_FUNC(int, writetxt, double* COMMA double* COMMA int COMMA char*);
#ifdef __cplusplus
}

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

@ -35,20 +35,42 @@
#ifndef BUILD_LIB
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_filter_iir filter_iir ;
p_goertzel goertzel ;
p_goertzel_cmplx goertzel_cmplx ;
p_linspace linspace ;
p_logspace logspace ;
p_polyval polyval ;
p_polyval_cmplx polyval_cmplx ;
p_acos_cmplx acos_cmplx ;
p_asin_cmplx asin_cmplx ;
p_cheby_poly1 cheby_poly1 ;
p_cheby_poly2 cheby_poly2 ;
p_cmplx2re cmplx2re ;
p_concat concat ;
p_conv conv ;
p_conv_cmplx conv_cmplx ;
p_cos_cmplx cos_cmplx ;
p_dft dft ;
p_dft_cmplx dft_cmplx ;
p_dmod dmod ;
p_farrow_lagrange farrow_lagrange ;
p_farrow_spline farrow_spline ;
p_filter_iir filter_iir ;
p_flipip flipip ;
p_flipip_cmplx flipip_cmplx ;
p_fourier_series_dec fourier_series_dec ;
p_fourier_series_rec fourier_series_rec ;
p_goertzel goertzel ;
p_goertzel_cmplx goertzel_cmplx ;
p_linspace linspace ;
p_log_cmplx log_cmplx ;
p_logspace logspace ;
p_polyval polyval ;
p_polyval_cmplx polyval_cmplx ;
p_randn randn ;
p_randu randu ;
p_re2cmplx re2cmplx ;
p_signal_pimp signal_pimp ;
p_sin_cmplx sin_cmplx ;
p_sqrt_cmplx sqrt_cmplx ;
p_trapint trapint ;
p_trapint_cmplx trapint_cmplx ;
p_writebin writebin ;
p_writetxt writetxt ;
#endif //BUILD_LIB
@ -100,23 +122,42 @@ void* dspl_load()
}
#endif //LINUX_OS
LOAD_FUNC(acos_cmplx);
LOAD_FUNC(asin_cmplx);
LOAD_FUNC(cheby_poly1);
LOAD_FUNC(cheby_poly2);
LOAD_FUNC(cmplx2re);
LOAD_FUNC(concat);
LOAD_FUNC(conv);
LOAD_FUNC(conv_cmplx);
LOAD_FUNC(cos_cmplx);
LOAD_FUNC(dft);
LOAD_FUNC(dft_cmplx);
LOAD_FUNC(dmod);
LOAD_FUNC(farrow_lagrange);
LOAD_FUNC(farrow_spline);
LOAD_FUNC(filter_iir);
LOAD_FUNC(flipip);
LOAD_FUNC(flipip_cmplx);
LOAD_FUNC(fourier_series_dec);
LOAD_FUNC(fourier_series_rec);
LOAD_FUNC(goertzel);
LOAD_FUNC(goertzel_cmplx);
LOAD_FUNC(linspace);
LOAD_FUNC(log_cmplx);
LOAD_FUNC(logspace);
LOAD_FUNC(polyval);
LOAD_FUNC(polyval_cmplx);
LOAD_FUNC(randn);
LOAD_FUNC(randu);
LOAD_FUNC(re2cmplx);
LOAD_FUNC(signal_pimp);
LOAD_FUNC(sin_cmplx);
LOAD_FUNC(sqrt_cmplx);
LOAD_FUNC(trapint);
LOAD_FUNC(trapint_cmplx);
LOAD_FUNC(writebin);
LOAD_FUNC(writetxt);
#ifdef WIN_OS

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

@ -191,24 +191,42 @@ extern "C" {
#endif
DECLARE_FUNC(int, acos_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, asin_cmplx, complex_t* COMMA int COMMA complex_t*);
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, cmplx2re, complex_t* COMMA int COMMA double* COMMA double*);
DECLARE_FUNC(int, concat, void* COMMA size_t COMMA void* COMMA size_t COMMA void*);
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, cos_cmplx, 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(double,dmod, double COMMA double);
DECLARE_FUNC(int, farrow_lagrange, double* COMMA int COMMA double COMMA double COMMA double COMMA double** COMMA int*);
DECLARE_FUNC(int, farrow_spline, double* COMMA int COMMA double COMMA double COMMA double COMMA double** COMMA int*);
DECLARE_FUNC(int, filter_iir, double* COMMA double* COMMA int COMMA double* COMMA int COMMA double*);
DECLARE_FUNC(int, flipip, double* COMMA int);
DECLARE_FUNC(int, flipip_cmplx, complex_t* COMMA int);
DECLARE_FUNC(int, fourier_series_dec, double* COMMA double* COMMA int COMMA double COMMA int COMMA double* COMMA complex_t*);
DECLARE_FUNC(int, fourier_series_rec, double* COMMA complex_t* COMMA int COMMA double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel, double* COMMA int COMMA int* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, goertzel_cmplx, complex_t* COMMA int COMMA int* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, linspace, double COMMA double COMMA int COMMA int COMMA double*);
DECLARE_FUNC(int, log_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, logspace, double COMMA double COMMA int COMMA int COMMA double*);
DECLARE_FUNC(int, polyval, double* COMMA int COMMA double* COMMA int COMMA double*);
DECLARE_FUNC(int, polyval_cmplx, complex_t* COMMA int COMMA complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, randn, double* COMMA int COMMA double COMMA double);
DECLARE_FUNC(int, randu, double* COMMA int);
DECLARE_FUNC(int, re2cmplx, double* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, signal_pimp, double* COMMA size_t COMMA double COMMA double COMMA double COMMA double COMMA double*);
DECLARE_FUNC(int, sin_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, sqrt_cmplx, complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, trapint, double* COMMA double* COMMA int COMMA double* sum);
DECLARE_FUNC(int, trapint_cmplx, double* COMMA complex_t* COMMA int COMMA complex_t*);
DECLARE_FUNC(int, writebin, void* COMMA int COMMA int COMMA char*);
DECLARE_FUNC(int, writetxt, double* COMMA double* COMMA int COMMA char*);
#ifdef __cplusplus
}

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

Plik binarny nie jest wyświetlany.