diff --git a/dspl/dox/ru/math.dox b/dspl/dox/ru/math.dox new file mode 100644 index 0000000..0bd0496 --- /dev/null +++ b/dspl/dox/ru/math.dox @@ -0,0 +1,58 @@ + +/*! **************************************************************************** +\ingroup SPEC_MATH_COMMON_GROUP +\fn int sinc(double* x, int n, double a, double* y) +\brief Функция \f$ \textrm{sinc}(x,a) = \frac{\sin(ax)}{ax}\f$ + +Функция рассчитывает занчения значения функции для вещественного вектора `x`. +
+ +\param[in] x Указатель на вектор переменной \f$ x \f$.
+ Размер вектора `[n x 1]`.
+ Память должна быть выделена.

+ +\param[in] n Размер входного вектора `x`.

+ +\param[in] a Параметр функции + \f$ \textrm{sinc}(x,a) = \frac{\sin(ax)}{ax}\f$ + +\param[out] y Указатель на вектор значений функции.
+ Размер вектора `[n x 1]`.
+ Память должна быть выделена.

+ + +\return + `RES_OK` Расчет произведен успешно.
+ В противном случае + \ref ERROR_CODE_GROUP "код ошибки".
+ + +Пример использования функции `sinc`: + +\include sinc_test.c + + +

+В каталоге `dat` будут созданы три файла:
+ +
+sinc_test_1.0.txt	  параметр a = 1.0
+sinc_test_pi.txt	  параметр a = pi
+sinc_test_2pi.txt	  параметр a = 2*pi
+
+ +Кроме того программа GNUPLOT произведет построение графика +по сохраненным в файлах данным: + +График: `img/sinc_test.png` +\image html sinc_test.png + +Скрипт GNUPLOT для построения графиков из текстовых файлов: +\include sinc_test.plt + +\author + Бахурин Сергей + www.dsplib.org + +***************************************************************************** */ + diff --git a/dspl/src/math.c b/dspl/src/math.c index ef1e13d..d0d9baf 100755 --- a/dspl/src/math.c +++ b/dspl/src/math.c @@ -41,7 +41,7 @@ double DSPL_API dmod (double x, double y) /******************************************************************************* sinc(x) = sin(pi*x)/(pi*x) *******************************************************************************/ -int DSPL_API sinc(double* x, int n, double* y) +int DSPL_API sinc(double* x, int n, double a, double* y) { int k; @@ -51,7 +51,7 @@ int DSPL_API sinc(double* x, int n, double* y) return ERROR_SIZE; for(k = 0; k < n; k++) - y[k] = (x[k]==0.0) ? 1.0 : sin(M_PI*x[k])/(M_PI*x[k]); + y[k] = (x[k]==0.0) ? 1.0 : sin(a*x[k])/(a*x[k]); return RES_OK; diff --git a/dspl/src/matrix.c b/dspl/src/matrix.c new file mode 100644 index 0000000..2cbd03c --- /dev/null +++ b/dspl/src/matrix.c @@ -0,0 +1,194 @@ +/* +* 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 "dspl.h" + + +void transpose(double* a, int n, int m, double* b); +void transpose_cmplx(complex_t* a, int n, int m, complex_t* b); +void transpose_hermite(complex_t* a, int n, int m, complex_t* b); + + + + +/******************************************************************************* +Real matrx transpose +*******************************************************************************/ +int DSPL_API matrix_create(matrix_t* a, int n, int m, int type) +{ + if(!a) + return ERROR_PTR; + if(n < 1 || m < 1) + return ERROR_MATRIX_SIZE; + + if(a->dat) + { + a->dat = (type & DAT_MASK) ? + (void*) realloc(a->dat, n*m*sizeof(complex_t)): + (void*) realloc(a->dat, n*m*sizeof(double)); + } + else + { + a->dat = (type & DAT_MASK) ? + (void*) malloc(n*m*sizeof(complex_t)): + (void*) malloc(n*m*sizeof(double)); + } + + a->n = n; + a->m = m; + a->type = type; + + return RES_OK; +} + + + + + + + + + + + +/******************************************************************************* +Real matrx transpose +*******************************************************************************/ +void transpose(double* a, int n, int m, double* b) +{ + int p, q, i, j, aind, bind; + + for(p = 0; p < n - DSPL_MATRIX_BLOCK; p+=DSPL_MATRIX_BLOCK) + { + for(q = 0; q < m - DSPL_MATRIX_BLOCK; q+=DSPL_MATRIX_BLOCK) + { + for(i = 0; i < DSPL_MATRIX_BLOCK; i++) + { + for(j = 0; j < DSPL_MATRIX_BLOCK; j++) + { + aind = (q+j) * n + p + i; + bind = (p+i) * m + q + j; + b[bind] = a[aind]; + } + } + } + } + for(i = p; i < n; i++) + for(j = 0; j < m; j++) + b[i*m + j] = a[j*n+i]; + + for(i = 0; i < p; i++) + for(j = q; j < m; j++) + b[i*m + j] = a[j*n+i]; + +} + + + + +/******************************************************************************* +Complex matrx transpose +*******************************************************************************/ +void transpose_cmplx(complex_t* a, int n, int m, complex_t* b) +{ + int p, q, i, j, aind, bind; + + for(p = 0; p < n - DSPL_MATRIX_BLOCK; p+=DSPL_MATRIX_BLOCK) + { + for(q = 0; q < m - DSPL_MATRIX_BLOCK; q+=DSPL_MATRIX_BLOCK) + { + for(i = 0; i < DSPL_MATRIX_BLOCK; i++) + { + for(j = 0; j < DSPL_MATRIX_BLOCK; j++) + { + aind = (q+j) * n + p + i; + bind = (p+i) * m + q + j; + RE(b[bind]) = RE(a[aind]); + IM(b[bind]) = IM(a[aind]); + } + } + } + } + for(i = p; i < n; i++) + { + for(j = 0; j < m; j++) + { + RE(b[i*m + j]) = RE(a[j*n+i]); + IM(b[i*m + j]) = IM(a[j*n+i]); + } + } + + for(i = 0; i < p; i++) + { + for(j = q; j < m; j++) + { + RE(b[i*m + j]) = RE(a[j*n+i]); + IM(b[i*m + j]) = IM(a[j*n+i]); + } + } +} + + + + +/******************************************************************************* +Hermite matrx transpose +*******************************************************************************/ +void transpose_hermite(complex_t* a, int n, int m, complex_t* b) +{ + int p, q, i, j, aind, bind; + + for(p = 0; p < n - DSPL_MATRIX_BLOCK; p+=DSPL_MATRIX_BLOCK) + { + for(q = 0; q < m - DSPL_MATRIX_BLOCK; q+=DSPL_MATRIX_BLOCK) + { + for(i = 0; i < DSPL_MATRIX_BLOCK; i++) + { + for(j = 0; j < DSPL_MATRIX_BLOCK; j++) + { + aind = (q+j) * n + p + i; + bind = (p+i) * m + q + j; + RE(b[bind]) = RE(a[aind]); + IM(b[bind]) = -IM(a[aind]); + } + } + } + } + for(i = p; i < n; i++) + { + for(j = 0; j < m; j++) + { + RE(b[i*m + j]) = RE(a[j*n+i]); + IM(b[i*m + j]) = IM(a[j*n+i]); + } + } + + for(i = 0; i < p; i++) + { + for(j = q; j < m; j++) + { + RE(b[i*m + j]) = RE(a[j*n+i]); + IM(b[i*m + j]) = -IM(a[j*n+i]); + } + } +} + diff --git a/include/dspl.c b/include/dspl.c index ae002ef..ef249ed 100755 --- a/include/dspl.c +++ b/include/dspl.c @@ -97,6 +97,7 @@ p_logspace logspace ; p_low2bp low2bp ; p_low2high low2high ; p_low2low low2low ; +p_matrix_create matrix_create ; p_poly_z2a_cmplx poly_z2a_cmplx ; p_polyval polyval ; p_polyval_cmplx polyval_cmplx ; @@ -229,6 +230,7 @@ void* dspl_load() LOAD_FUNC(low2bp); LOAD_FUNC(low2high); LOAD_FUNC(low2low); + LOAD_FUNC(matrix_create); LOAD_FUNC(poly_z2a_cmplx); LOAD_FUNC(polyval); LOAD_FUNC(polyval_cmplx); diff --git a/include/dspl.h b/include/dspl.h index 71e3e21..42d1a3f 100755 --- a/include/dspl.h +++ b/include/dspl.h @@ -43,15 +43,23 @@ typedef double complex_t[2]; typedef struct { - complex_t* w; - complex_t* t0; - complex_t* t1; - int n; - int p2; + complex_t* w; + complex_t* t0; + complex_t* t1; + int n; + int p2; } fft_t; +typedef struct +{ + void* dat; + int n; + int m; + int type; +} matrix_t; + #define RE(x) (x[0]) #define IM(x) (x[1]) @@ -99,6 +107,7 @@ typedef struct /* K 0x11xxxxxx*/ /* L 0x12xxxxxx*/ /* M 0x13xxxxxx*/ +#define ERROR_MATRIX_SIZE 0x13011909 /* N 0x14xxxxxx*/ #define ERROR_NEGATIVE 0x14050701 /* O 0x15xxxxxx*/ @@ -124,11 +133,11 @@ typedef struct /* Y 0x25xxxxxx*/ /* Z 0x26xxxxxx*/ +#define DAT_MASK 0x00000001 +#define DAT_DOUBLE 0x00000000 +#define DAT_COMPLEX 0x00000001 -#define DAT_DOUBLE 0 -#define DAT_COMPLEX 1 - - +#define DSPL_MATRIX_BLOCK 32 #define DSPL_SYMMETRIC 0x00000000 @@ -547,6 +556,11 @@ DECLARE_FUNC(int, low2low, double* b COMMA double* beta COMMA double* alpha); //------------------------------------------------------------------------------ +DECLARE_FUNC(int, matrix_create, matrix_t* a + COMMA int n + COMMA int m + COMMA int type); +//------------------------------------------------------------------------------ DECLARE_FUNC(int, poly_z2a_cmplx, complex_t* COMMA int COMMA int @@ -606,6 +620,7 @@ DECLARE_FUNC(int, sin_cmplx, complex_t* //------------------------------------------------------------------------------ DECLARE_FUNC(int, sinc, double* x COMMA int n + COMMA double a COMMA double* y); //------------------------------------------------------------------------------ DECLARE_FUNC(int, sqrt_cmplx, complex_t* diff --git a/release/include/dspl.c b/release/include/dspl.c index ae002ef..ef249ed 100755 --- a/release/include/dspl.c +++ b/release/include/dspl.c @@ -97,6 +97,7 @@ p_logspace logspace ; p_low2bp low2bp ; p_low2high low2high ; p_low2low low2low ; +p_matrix_create matrix_create ; p_poly_z2a_cmplx poly_z2a_cmplx ; p_polyval polyval ; p_polyval_cmplx polyval_cmplx ; @@ -229,6 +230,7 @@ void* dspl_load() LOAD_FUNC(low2bp); LOAD_FUNC(low2high); LOAD_FUNC(low2low); + LOAD_FUNC(matrix_create); LOAD_FUNC(poly_z2a_cmplx); LOAD_FUNC(polyval); LOAD_FUNC(polyval_cmplx); diff --git a/release/include/dspl.h b/release/include/dspl.h index 71e3e21..42d1a3f 100755 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -43,15 +43,23 @@ typedef double complex_t[2]; typedef struct { - complex_t* w; - complex_t* t0; - complex_t* t1; - int n; - int p2; + complex_t* w; + complex_t* t0; + complex_t* t1; + int n; + int p2; } fft_t; +typedef struct +{ + void* dat; + int n; + int m; + int type; +} matrix_t; + #define RE(x) (x[0]) #define IM(x) (x[1]) @@ -99,6 +107,7 @@ typedef struct /* K 0x11xxxxxx*/ /* L 0x12xxxxxx*/ /* M 0x13xxxxxx*/ +#define ERROR_MATRIX_SIZE 0x13011909 /* N 0x14xxxxxx*/ #define ERROR_NEGATIVE 0x14050701 /* O 0x15xxxxxx*/ @@ -124,11 +133,11 @@ typedef struct /* Y 0x25xxxxxx*/ /* Z 0x26xxxxxx*/ +#define DAT_MASK 0x00000001 +#define DAT_DOUBLE 0x00000000 +#define DAT_COMPLEX 0x00000001 -#define DAT_DOUBLE 0 -#define DAT_COMPLEX 1 - - +#define DSPL_MATRIX_BLOCK 32 #define DSPL_SYMMETRIC 0x00000000 @@ -547,6 +556,11 @@ DECLARE_FUNC(int, low2low, double* b COMMA double* beta COMMA double* alpha); //------------------------------------------------------------------------------ +DECLARE_FUNC(int, matrix_create, matrix_t* a + COMMA int n + COMMA int m + COMMA int type); +//------------------------------------------------------------------------------ DECLARE_FUNC(int, poly_z2a_cmplx, complex_t* COMMA int COMMA int @@ -606,6 +620,7 @@ DECLARE_FUNC(int, sin_cmplx, complex_t* //------------------------------------------------------------------------------ DECLARE_FUNC(int, sinc, double* x COMMA int n + COMMA double a COMMA double* y); //------------------------------------------------------------------------------ DECLARE_FUNC(int, sqrt_cmplx, complex_t* diff --git a/test/bin/gnuplot/sinc_test.plt b/test/bin/gnuplot/sinc_test.plt new file mode 100644 index 0000000..dfef9c1 --- /dev/null +++ b/test/bin/gnuplot/sinc_test.plt @@ -0,0 +1,13 @@ +set grid +set xlabel "x" + +set lmargin at screen 0.10 + +set terminal pngcairo size 560,280 enhanced font 'Verdana,8' +set output 'img/sinc_test.png' +set ylabel "sinc(x,a)" +set yrange [-0.25:1.1] +plot 'dat/sinc_test_1.0.txt' with lines title "a = 1.0", \ +'dat/sinc_test_pi.txt' with lines title "a = pi", \ +'dat/sinc_test_2pi.txt' with lines title "a = 2pi" + diff --git a/test/src/sinc_test.c b/test/src/sinc_test.c new file mode 100644 index 0000000..b805558 --- /dev/null +++ b/test/src/sinc_test.c @@ -0,0 +1,32 @@ +#include +#include +#include +#include "dspl.h" + +#define N 1000 + + +int main() +{ + void* handle; // DSPL handle + handle = dspl_load(); // Load DSPL function + + double x[N], y[N]; + + linspace(-10.0, 10.0, N , DSPL_SYMMETRIC, x); + sinc(x, N, 1.0, y); + writetxt(x, y, N, "dat/sinc_test_1.0.txt"); + + sinc(x, N, M_PI, y); + writetxt(x, y, N, "dat/sinc_test_pi.txt"); + + sinc(x, N, 2.0*M_PI, y); + writetxt(x, y, N, "dat/sinc_test_2pi.txt"); + + dspl_free(handle); // free dspl handle + + // выполнить скрипт GNUPLOT для построения графиков + // по рассчитанным данным + return system("gnuplot gnuplot/sinc_test.plt");; +} +