From 26c3c73692c0edbc6e8d885ecf25e9531dd600b4 Mon Sep 17 00:00:00 2001 From: Dsplib Date: Fri, 16 Mar 2018 00:01:34 +0300 Subject: [PATCH] added some functions --- dspl/dox/ru/fourier_series.dox | 0 dspl/dox/ru/polyval.dox | 21 +-- dspl/dox/ru/randgen.dox | 63 +++++++-- dspl/src/array.c | 101 ++++++++++++++ dspl/src/complex.c | 236 +++++++++++++++++++++++++++++++++ dspl/src/fourier_series.c | 100 ++++++++++++++ dspl/src/inout.c | 129 ++++++++++++++++++ dspl/src/math.c | 33 +++++ dspl/src/randgen.c | 119 +++++++++++++++++ dspl/src/resampling.c | 187 ++++++++++++++++++++++++++ dspl/src/signals.c | 60 +++++++++ dspl/src/trapint.c | 70 ++++++++++ include/dspl.c | 77 ++++++++--- include/dspl.h | 28 +++- release/include/dspl.c | 77 ++++++++--- release/include/dspl.h | 28 +++- test/bin/dft_test.exe | Bin 13568 -> 14344 bytes 17 files changed, 1256 insertions(+), 73 deletions(-) create mode 100644 dspl/dox/ru/fourier_series.dox create mode 100644 dspl/src/array.c create mode 100644 dspl/src/complex.c create mode 100644 dspl/src/fourier_series.c create mode 100644 dspl/src/inout.c create mode 100644 dspl/src/math.c create mode 100644 dspl/src/randgen.c create mode 100644 dspl/src/resampling.c create mode 100644 dspl/src/signals.c create mode 100644 dspl/src/trapint.c diff --git a/dspl/dox/ru/fourier_series.dox b/dspl/dox/ru/fourier_series.dox new file mode 100644 index 0000000..e69de29 diff --git a/dspl/dox/ru/polyval.dox b/dspl/dox/ru/polyval.dox index 648bf03..251dbf0 100644 --- a/dspl/dox/ru/polyval.dox +++ b/dspl/dox/ru/polyval.dox @@ -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 Указатель на вектор реальной части коэффициентов полинома.
+ \param[in] a Указатель на вектор комплексных коэффициентов полинома.
Размер вектора `[ord+1 x 1]`.
- Коэффициент `aR[0]` соответствует коэффициенту полинома \f$a_0\f$.

+ Коэффициент `a[0]` соответствует коэффициенту полинома \f$a_0\f$.

- \param[in] aI Указатель на вектор мнимой части коэффициентов полинома.
- Размер вектора `[ord+1 x 1]`.
- Коэффициент `aI[0]` соответствует коэффициенту полинома \f$a_0\f$.

\param[in] ord Порядок полинома \f$N\f$.

- \param[in] xR Указатель на вектор реальной части аргумента полинома.
- Размер вектора `[n x 1]`.
- Значения полинома будут расчитаны для всех - значений аргумента вектора `x`.

- - \param[in] xI Указатель на вектор мнимой части аргумента полинома.
+ \param[in] x Указатель на вектор аргумента полинома.
Размер вектора `[n x 1]`.
Значения полинома будут расчитаны для всех значений аргумента вектора `x`.

\param[in] n Размер вектора агрумента полинома.

- \param[out] yR Указатель вектор реальной части значения - полинома для аргумента `x`.
+ \param[out] y Указатель вектор значения полинома для аргумента `x`.
Размер вектора `[n x 1]`.
Память должна быть выделена.

- \param[out] yI Указатель вектор реальной части значения - полинома для аргумента `x`.
- Размер вектора `[n x 1]`.
- Память должна быть выделена.

\return `RES_OK` Полином расчитан успешно.
diff --git a/dspl/dox/ru/randgen.dox b/dspl/dox/ru/randgen.dox index 1383def..610ca31 100644 --- a/dspl/dox/ru/randgen.dox +++ b/dspl/dox/ru/randgen.dox @@ -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 Генерация вектора нормально распределенных псевдослучайных чисел с + заданным математическим ожиданием и среднеквадратическим отклонением. + + Генерация производится при помощи преобразования Бокса — Мюллера равномерно-распределенной + случайной величины в нормально распределенную.
+ + + + \param[in,out] x Указатель на вектор нормальной распределенных случайных чисел.
+ Размер вектора `[n x 1]`.
+ Память должна быть выделена.

+ + \param[in] n Размер вектора случайных чисел.

+ + \param[in] mu Математическое ожидание.

+ + \param[in] sigma Среднеквадратическое отклонение (СКО).

+ + \return + `RES_OK` Вектор случайных чисел сгенерирован успешно.
+ В противном случае \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$.
+ + + \param[in,out] x Указатель на вектор случайных чисел.
+ Размер вектора `[n x 1]`.
+ Память должна быть выделена.

+ + \param[in] n Размер вектора случайных чисел.

+ + \return + `RES_OK` Вектор случайных чисел сгенерирован успешно.
+ В противном случае \ref ERROR_CODE_GROUP "код ошибки". + + \author + Бахурин Сергей. + www.dsplib.org + +**************************************************************************************************** */ -/************************************************************************************************** -Gaussian random numbers generator - -int randn(double* x, int n, double mu, double sigma) -***************************************************************************************************/ diff --git a/dspl/src/array.c b/dspl/src/array.c new file mode 100644 index 0000000..f1032f3 --- /dev/null +++ b/dspl/src/array.c @@ -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 . +*/ + + +#include +#include +#include +#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; +} + + + + + diff --git a/dspl/src/complex.c b/dspl/src/complex.c new file mode 100644 index 0000000..f6f7517 --- /dev/null +++ b/dspl/src/complex.c @@ -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 . +*/ + + +#include +#include +#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; +} + + + + diff --git a/dspl/src/fourier_series.c b/dspl/src/fourier_series.c new file mode 100644 index 0000000..9c9b87f --- /dev/null +++ b/dspl/src/fourier_series.c @@ -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 . +*/ + + +#include +#include +#include +#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; + +} diff --git a/dspl/src/inout.c b/dspl/src/inout.c new file mode 100644 index 0000000..6de0363 --- /dev/null +++ b/dspl/src/inout.c @@ -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 . +*/ + + +#include +#include +#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; +} + + diff --git a/dspl/src/math.c b/dspl/src/math.c new file mode 100644 index 0000000..e814b02 --- /dev/null +++ b/dspl/src/math.c @@ -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 . +*/ + + +#include +#include +#include "dspl.h" + + + +double DSPL_API dmod (double x, double y) +{ + if(y == 0.0) + return x; + return x - floor(x/y) * y; +} diff --git a/dspl/src/randgen.c b/dspl/src/randgen.c new file mode 100644 index 0000000..571ea35 --- /dev/null +++ b/dspl/src/randgen.c @@ -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 . +*/ + + +#include +#include +#include + +#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 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. +*/ + + + +#include +#include +#include +#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; +} + + + diff --git a/dspl/src/signals.c b/dspl/src/signals.c new file mode 100644 index 0000000..e3e352d --- /dev/null +++ b/dspl/src/signals.c @@ -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 . +*/ + + + +#include +#include +#include +#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; +} + + diff --git a/dspl/src/trapint.c b/dspl/src/trapint.c new file mode 100644 index 0000000..898df8f --- /dev/null +++ b/dspl/src/trapint.c @@ -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 . +*/ + + +#include +#include +#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; +} diff --git a/include/dspl.c b/include/dspl.c index 5fd8953..38e1d06 100644 --- a/include/dspl.c +++ b/include/dspl.c @@ -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 diff --git a/include/dspl.h b/include/dspl.h index 01403ef..adf2c99 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -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 } diff --git a/release/include/dspl.c b/release/include/dspl.c index 5fd8953..38e1d06 100644 --- a/release/include/dspl.c +++ b/release/include/dspl.c @@ -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 diff --git a/release/include/dspl.h b/release/include/dspl.h index 01403ef..adf2c99 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -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 } diff --git a/test/bin/dft_test.exe b/test/bin/dft_test.exe index 2cb8e662b4783bec01e49126b113fa972b311dcd..de5dff4c02f017e6af0ebb50058697ccc8f81c7f 100755 GIT binary patch delta 4170 zcmZve4Nz3q702&BK3p;ET^1Cs8a@7;G^yklp$zg+WH?Wy>A*3vv?9!vA~oT<_J(}ctFqGQk@#X*P1 z>{Z`14}$4d0JE9Z)ehK2{V{pqnWD`R{~U5Wh)w`Yt~;Kbkq z3aT(TN`VK1k11%t;1ddb7<@)S7=x2$&|4CRE14Ph&xm7uNpH!}3<2YVN<5~-$CY?e ziKie=&V8c5$#_v5`mFCu;|uO#?^&QfqQT#e_buKX<1lBj$~sTGl_}O)J+qGk;yrTM z$HgnwOR;qh4@|Vmq*k;J}%*q=-+ zHa3bT*WX0$AxihU$#ZDPuMJpjbXYugrUdvOYnpkA;Knam`idLvweYmSQ7NO z{yt(qkXS77JIQq*_bzd0a{X20wo|&j63fL=e?k~~T;Gb=Zo)Q6Y(HiT5r!Vu8xVVq zusVqy$cq8zB@RulFG21t;ucEoZtA@(j|(_fi$eFm``*oV2I$Mw(QnK3>l z>=TL2N9=^eVxwqs9Upq*Byn#@?r1J>zoK;LaUC~-ah9;{5^F|GCk#EV<0dc?gl&-6 zF~puF3_Y&HCWs}CDdLt(ZYgs26NetxaT6GK2)jpOFCk_p3_Y&nCNR?BwSe38!gX^= zo<-~;_F+lT<2r5vBa5(@A-N+C;C@dWnq0?CVB``vAhAKjULg!UuHz;!ZrSCiVToCA z)F@%-aUC~-QAFA5C3eN(jIqSy#G%P`+yq7iaf>DQc@Cs2BMd#R<0dew3A_1kbFP0t zEQx)XD|%eVO<;Hkn~>N`h@B=Zh9=i>6BzZx9gy6s$Q_bgv~U#Ni+g};>wm8a`iU~I z@MY`mJB@E#*5=|xy;ug=#mB(98}`>eQ~RCTzGUtC+NamYZ>J>&rib}Ge#=~cOQ^f4WplW*hxxmMU9u}e*@ZULsv`l`66$L4 zN5v1bi`rW^hgw;iKN1OTZR+&5NBmvwfmGZb?hJN86zq%!B2B?ygtc`B!@)2$G)7x! zOC$(>cL191rq)0UQ=~|sMI6j7&7>>&G~4+UjMLi1@P7-71ghoLb_d(L{GCnV;N~#v zel8NF0iqFqIM@|s)RMY25{w3#gI&z*M0=uQD!Z;87SXaG?ZG-|XTfjeYBJdmZR5X^ z$q{H1*OEyew$sVv6tp9-LFPl%#rjtT?E~;~=)o$?#PhbrTKt-L+2+!mSH(xRGHvjx zxMHhWJq0<&q1>`sn|+6-ZOF1&ADNT({S0Vp=hWG(d0_2_-yl3`h@}p^rgz}ie^Yeu zM>9si)`k!Bzt7&gs6m_H!|N6qY2TiJ8m}H{%+4miR?lIn@2HvQq%jAZbPR7;Y}2&z zJke3ITZ7N68zrtB{DQy(KQZvQUmPmUx2zGbmF2rPgXm-htkLN+z0_X|EJqG2WId$A zN0WJUh;ya+1uu{rSKMBaQm~3SEPTa=ozp|xEjY5-6r+b?6hLcr4uF<-AQ^v(JO@i zzWi7{xl@Y!H0jqBy^(ZQ>@#mp>LzUh-TpNE7POyq!r61cW1hIx zV@efd{Dtz&SMt10H_Tb=vzcgWfD3q8D{sJQILt<6GMQ;|n5j%T%#R=q9&=C^pIUNb zrO?AWUakzdg7hjyUqku>iv9%Y9z|bIx>wPgrS4)4im?$4T$u)1;nX`O0=g3~`SBO= zRJe2MWanveDve6Uy>ck#Q}lsU#=s#WHY&v1q=yy#1JZjG{WH?ziarUtx!eU;v4{74 zw#U2?=9?CJe-wwx^D}OzxGV=ec6eKM77u?}{)mMuB}}*~7S>a>3U7&adpv8_)~|lJ zsoGWTTA8XOvUFe#X$rLZqyCgm73OXo6)hDtwW&J8!lBM>&qJO1|DyHZN+t{JjPlDm zL+xUyVx>bVjZiJcL`7tlnX5V?>MEaHt{_l7Q$te0l=GtT%}OeMR#};jrCh8^zjs+% zs&dXKqAeQ(&D)yr(yRZsybn)T_F0DG%Y13td8;^I<;qIEp(nT~bk{girDu|xH-!C0 z6U}Z{R;uzqX?fXQlQjunHh4hcIpfaMEal>yyJ+uj@ZiWvJ3*9JMg delta 2221 zcmZvde@qi+7{}kcQY{5~^r#kuNg>4*r441ZEc%C65LPxe9RgFg1d2kGbg{T_f{B>(2>9Zs&xd>u&xUPetxc(PdDAE`uZwTn2iKyKV~l zJ?`m{ULf*)EZ)7ZCMkAt~T5KVSd<7aoJ6IIUg-0*%bm zD3RD+d71IV53;GfifjCg<8XWgyn7FvB_&tP-(KhM^~fpF`{b!`3L+(;?)l8Hc72wjy_!ak&b2c@ba*3`0){ z&piO_B*Uivlb_QpVt?ZpYD7;6|Ag2G!@g0lC}J07EEFF>QwSeN?rX+LN3Bkss@`N`LqE-qm@}&S?_n+3351k27!C*500+gccn!qt7BUdgG|9t4|@7MaG=qva-l+Ivga^NpF& zhNUshvyE^f)lEaPN$X?dMw0sK&7Y$oBI!0nHy0atJ}%Hs>wX?SinCUS5zjj|ECE|3 zeX+!xHx89sNFLc}X_cd?+XmXT%*?gZn;zqKB1lU|1o7_d;C^2)d z%PYK3Pr@#rFtT4o9!I9`B69~-?k_Aptcu@&c=ls*o}?bgc?Rp^v~>DObDjvr8s|O1 zG?rlzEwY(?E5Q|sflOpba=e|nHS(+_pw~ zBimHDexD+ocaVT;;;Oa(V(|;AcoH*piW*jHx+Lc>p9mPOfmtHSAU0N?Y{IR88cWT3 zC(OdVE?J`DiW-jI*A(Oib_HLPq-W~us;XNm92Jh0q@{hUFX*G+mAMzE6v(l+m`zT7ILB9;