From c79c129fec2bc98a1f501aae696dfb04d0298b4b Mon Sep 17 00:00:00 2001 From: Dsplib Date: Sat, 5 May 2018 19:51:32 +0300 Subject: [PATCH] before merging --- dox/ru/mainpage.dox | 26 ++------ dspl/dox/ru/complex.dox | 137 ++++++++++++++++++++++++++++++++++++++ dspl/dox/ru/randgen.dox | 5 +- dspl/src/dft.c | 45 ++++++++++++- dspl/src/fft.c | 87 ++++++++++++++++++++++++ dspl/src/filter_an.c | 63 ++++++++++++++++++ include/dspl.c | 8 +++ include/dspl.h | 6 +- release/include/dspl.c | 8 +++ release/include/dspl.h | 6 +- test/src/butter_ap_test.c | 4 +- 11 files changed, 369 insertions(+), 26 deletions(-) diff --git a/dox/ru/mainpage.dox b/dox/ru/mainpage.dox index e93582c..674d508 100644 --- a/dox/ru/mainpage.dox +++ b/dox/ru/mainpage.dox @@ -14,31 +14,19 @@ \ref FILTER_ANALYSIS_GROUP
\ref RESAMPLING_GROUP
\ref SPEC_MATH_GROUP
- \ref IN_OUT_GROUP
- + \ref IN_OUT_GROUP

- + Библиотека поддерживает работу с вещественными и комплексными типами входных данных.
+ Описание испльзуемых типов:
+ + \ref TYPES_GROUP
DSPL-2.0 библиотека с открытым исходным кодом, написанная на языке Си.
- DSPL-2.0 использует наиболее эффективные сторонние библиотеки:
- libfftw3 библиотека алгоритмов быстрого преобразования Фурье.
- liblapack библиотека алгоритмов линейной алгебры.

- - - DSPL-2.0 содержит собственные алгоритмы цифровой обработки сигналов, а также реализует интерфейс - для работы со сторонними библиотеками как это показано на структуре:
- - \image html dspl_structure.png - - В дополнение к этому, совместно с DSPL-2.0, релизована библиотека libplot, которая использует - ANSI C интерфейс с пакетом построения графиков - GNUPLOT. - - Pаспространяется под лицензией GPL v3 - Исходные коды библиотеки доступны на GitHub.
+ Pаспространяется под лицензией LGPL v3 + Исходные коды библиотеки доступны на GitHub.
Вы также можете внести свой вклад в развитие данной библиотеки. Присоединяйтесь! */ diff --git a/dspl/dox/ru/complex.dox b/dspl/dox/ru/complex.dox index af1ec83..3d65ac8 100644 --- a/dspl/dox/ru/complex.dox +++ b/dspl/dox/ru/complex.dox @@ -1,3 +1,140 @@ +/*! ************************************************************************************************* + \ingroup TYPES_GROUP + \typedef complex_t + \brief Описание комплексного типа данных. + + Комплексный тип данных в библиотеке DSPL-2.0 определен как массив из двух элементов типа `double`. + При этом первый элемент массива определяет реальную часть комплексного числа, а второй - мнимую. + + Например: + + \code + complex_t z; + z[0] = 1.0; + z[1] = -2.0; + \endcode + + Переменная `z = 1-2j`, где `j` - мнимая единица. + + Для удобства работы с комплексными числами реализованы специальные макросы: + \ref RE, \ref IM, \ref ABSSQR + + + + + +**************************************************************************************************** */ + + + + + + + + + +/*! ************************************************************************************************* + \ingroup TYPES_GROUP + \def ABSSQR(x) + \brief Макрос возвращает квадрат модуля комплексного числа `x`. + + Квадрата модуля комплексного числа \f$ x = a + j \cdot b \f$ равен: + + \f[ + |x|^2 = x \cdot x^* = a^2 + b^2. + \f] + + Например: + \code + complex_t z; + double y; + RE(z) = 1.0; + IM(z) = -2.0; + y = ABSSQR(z); + \endcode + + Переменная `z = 1-2j`, где `j` - мнимая единица, а переменная `y = 5`. + + + + +**************************************************************************************************** */ + + + + + + + + +/*! ************************************************************************************************* + \ingroup TYPES_GROUP + \def IM(x) + \brief Макрос определяющий мнимую часть комплексного числа. + + Например: + \code + complex_t z; + RE(z) = 1.0; + IM(z) = -2.0; + \endcode + + Переменная `z = 1-2j`, где `j` - мнимая единица. + + Аналогично, макрос можно использовать для получения мнимой части комплексного числа: + + \code + complex_t z = {3.0, -4.0}; + double r; + r = IM(z); + \endcode + В данном примере переменная `z = 3-4i`, а в переменой `r` будет храниться число -4. + + +**************************************************************************************************** */ + + + + + + + + + + +/*! ************************************************************************************************* + \ingroup TYPES_GROUP + \def RE(x) + \brief Макрос определяющий реальную часть комплексного числа. + + Например: + \code + complex_t z; + RE(z) = 1.0; + IM(z) = -2.0; + \endcode + + Переменная `z = 1-2j`, где `j` - мнимая единица. + + Аналогично, макрос можно использовать для получения реальной части комплексного числа: + + \code + complex_t z = {3.0, -4.0}; + double r; + r = RE(z); + \endcode + В данном примере переменная `z = 3-4i`, а в переменой `r` будет храниться число 3. + + +**************************************************************************************************** */ + + + + + + + + /*! ************************************************************************************************* \ingroup TYPES_GROUP diff --git a/dspl/dox/ru/randgen.dox b/dspl/dox/ru/randgen.dox index 610ca31..924fd80 100644 --- a/dspl/dox/ru/randgen.dox +++ b/dspl/dox/ru/randgen.dox @@ -32,12 +32,15 @@ **************************************************************************************************** */ + + + /*! ************************************************************************************************* \ingroup SPEC_MATH_RAND_GEN_GROUP \fn int randu(double* x, int n); \brief Генерация вектора равномерно-распределенных в интервале от 0 до 1 псевдослучайных чисел. - Генерация производится при помощи рекурсивного алгоритма L'Ecluyer. Период датчика порядка \f$10^56\f$.
+ Генерация производится при помощи рекурсивного алгоритма L'Ecluyer. Период датчика порядка \f$10^{56}\f$.
\param[in,out] x Указатель на вектор случайных чисел.
diff --git a/dspl/src/dft.c b/dspl/src/dft.c index ccf26ef..cf41843 100644 --- a/dspl/src/dft.c +++ b/dspl/src/dft.c @@ -49,8 +49,8 @@ int DSPL_API dft(double* x, int n, complex_t *y) RE(y[k]) = IM(y[k]) = 0.0; for(m = 0; m < n; m++) { - phi = -M_2PI * divn * (double)k * (double)m; - RE(y[k]) += x[m] * cos(phi); + phi = -M_2PI * divn * (double)k * (double)m; + RE(y[k]) += x[m] * cos(phi); IM(y[k]) += x[m] * sin(phi); } } @@ -92,3 +92,44 @@ int DSPL_API dft_cmplx(complex_t* x, int n, complex_t *y) } return RES_OK; } + + + + + +/************************************************************************************************** +Complex vector DFT +***************************************************************************************************/ +int DSPL_API idft_cmplx(complex_t* x, int n, complex_t *y) +{ + int k; + int m; + double divn; + double phi; + complex_t e; + + if(!x || !y) + return ERROR_PTR; + + if(n<1) + return ERROR_SIZE; + + divn = 1.0 / (double)n; + + for(k = 0; k < n; k++) + { + RE(y[k]) = IM(y[k]) = 0.0; + for(m = 0; m < n; m++) + { + phi = M_2PI * divn * (double)k * (double)m; + RE(e) = cos(phi); + IM(e) = sin(phi); + RE(y[k]) += CMRE(x[m], e); + IM(y[k]) += CMIM(x[m], e); + } + RE(y[k]) /= (double)n; + IM(y[k]) /= (double)n; + } + return RES_OK; +} + diff --git a/dspl/src/fft.c b/dspl/src/fft.c index 68ac16a..eb1886f 100644 --- a/dspl/src/fft.c +++ b/dspl/src/fft.c @@ -29,6 +29,47 @@ int fft_p2(int n); + +/************************************************************************************************** +COMPLEX vector IFFT +**************************************************************************************************/ +int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y) +{ + int err, k; + double norm; + + if(!x || !pfft || !y) + return ERROR_PTR; + if(n<1) + return ERROR_SIZE; + + + err = fft_create(pfft, n); + if(err != RES_OK) + return err; + + memcpy(pfft->t1, x, n*sizeof(complex_t)); + for(k = 0; k < n; k++) + IM(pfft->t1[k]) = -IM(pfft->t1[k]); + + err = fft_dit(pfft, n, y); + + if(err!=RES_OK) + return err; + + norm = 1.0 / (double)n; + for(k = 0; k < n; k++) + { + RE(y[k]) = RE(y[k])*norm; + IM(y[k]) = -IM(y[k])*norm; + } + return RES_OK; +} + + + + + /************************************************************************************************** Real vector FFT **************************************************************************************************/ @@ -335,6 +376,52 @@ int DSPL_API fft_shift(double* x, int n, double* y) +/************************************************************************************************** +FFT shifting for complex vector +***************************************************************************************************/ +int DSPL_API fft_shift_cmplx(complex_t* x, int n, complex_t* y) +{ + int n2, r; + int k; + complex_t tmp; + complex_t *buf; + + if(!x || !y) + return ERROR_PTR; + + if(n<1) + return ERROR_SIZE; + + r = n%2; + if(!r) + { + n2 = n>>1; + for(k = 0; k < n2; k++) + { + RE(tmp) = RE(x[k]); + IM(tmp) = IM(x[k]); + + RE(y[k]) = RE(x[k+n2]); + IM(y[k]) = IM(x[k+n2]); + + RE(y[k+n2]) = RE(tmp); + IM(y[k+n2]) = IM(tmp); + } + } + else + { + n2 = (n-1) >> 1; + buf = (complex_t*) malloc(n2*sizeof(complex_t)); + memcpy(buf, x, n2*sizeof(complex_t)); + memcpy(y, x+n2, (n2+1)*sizeof(complex_t)); + memcpy(y+n2+1, buf, n2*sizeof(complex_t)); + free(buf); + } + return RES_OK; +} + + + diff --git a/dspl/src/filter_an.c b/dspl/src/filter_an.c index 9b77eed..43c5d46 100644 --- a/dspl/src/filter_an.c +++ b/dspl/src/filter_an.c @@ -92,6 +92,69 @@ exit_label: +/************************************************************************************************** +impulse response of an analog filter H(s) +***************************************************************************************************/ +int DSPL_API freqs2time(double* b, double* a, int ord, double fs, int n, fft_t* pfft, double *t, double *h) +{ + double *w = NULL; + complex_t *hs = NULL; + complex_t *ht = NULL; + int err, k; + + if(!b || !a || !t || !h) + return ERROR_PTR; + if(ord<1) + return ERROR_FILTER_ORD; + if(n<1) + return ERROR_SIZE; + + w = (double*)malloc(n*sizeof(double)); + hs = (complex_t*)malloc(n*sizeof(complex_t)); + + + err = linspace(-fs*0.5, fs*0.5, n, DSPL_PERIODIC, w); + if(err != RES_OK) + goto exit_label; + + err = freqs(b, a, ord, w, n, hs); + if(err != RES_OK) + goto exit_label; + + err = fft_shift_cmplx(hs, n, hs); + if(err != RES_OK) + goto exit_label; + + ht = (complex_t*)malloc(n*sizeof(complex_t)); + + err = ifft_cmplx(hs, n, pfft, ht); + if(err != RES_OK) + { + err = idft_cmplx(hs, n, ht); + if(err != RES_OK) + goto exit_label; + } + + for(k = 0; k < n; k++) + { + t[k] = (double)k/fs; + h[k] = RE(ht[k]) * fs; + } + +exit_label: + if(w) + free(w); + if(hs) + free(hs); + if(ht) + free(ht); + return err; + +} + + + + /************************************************************************************************** Magnitude, phase response and group delay of an analog filter H(s) diff --git a/include/dspl.c b/include/dspl.c index ef2f307..c01180c 100644 --- a/include/dspl.c +++ b/include/dspl.c @@ -51,6 +51,8 @@ p_cos_cmplx cos_cmplx ; p_dft dft ; p_dft_cmplx dft_cmplx ; p_dmod dmod ; +p_idft_cmplx idft_cmplx ; +p_ifft_cmplx ifft_cmplx ; p_farrow_lagrange farrow_lagrange ; p_farrow_spline farrow_spline ; p_filter_iir filter_iir ; @@ -60,12 +62,14 @@ p_fft_cmplx fft_cmplx ; p_fft_create fft_create ; p_fft_free fft_free ; p_fft_shift fft_shift ; +p_fft_shift_cmplx fft_shift_cmplx ; p_flipip flipip ; p_flipip_cmplx flipip_cmplx ; p_fourier_series_dec fourier_series_dec ; p_fourier_series_rec fourier_series_rec ; p_freqs freqs ; p_freqs_resp freqs_resp ; +p_freqs2time freqs2time ; p_freqz freqz ; p_goertzel goertzel ; p_goertzel_cmplx goertzel_cmplx ; @@ -153,6 +157,8 @@ void* dspl_load() LOAD_FUNC(dft); LOAD_FUNC(dft_cmplx); LOAD_FUNC(dmod); + LOAD_FUNC(idft_cmplx); + LOAD_FUNC(ifft_cmplx); LOAD_FUNC(farrow_lagrange); LOAD_FUNC(farrow_spline); LOAD_FUNC(filter_iir); @@ -162,6 +168,7 @@ void* dspl_load() LOAD_FUNC(fft_create); LOAD_FUNC(fft_free); LOAD_FUNC(fft_shift); + LOAD_FUNC(fft_shift_cmplx); LOAD_FUNC(flipip); LOAD_FUNC(flipip_cmplx); LOAD_FUNC(fourier_series_dec); @@ -169,6 +176,7 @@ void* dspl_load() LOAD_FUNC(freqz); LOAD_FUNC(freqs); LOAD_FUNC(freqs_resp); + LOAD_FUNC(freqs2time); LOAD_FUNC(goertzel); LOAD_FUNC(goertzel_cmplx); LOAD_FUNC(linspace); diff --git a/include/dspl.h b/include/dspl.h index 1c0288e..88d6e9c 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -213,6 +213,8 @@ DECLARE_FUNC(int, cos_cmplx, complex_t* COMMA int COMMA comple 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, idft_cmplx, complex_t* COMMA int COMMA complex_t*); +DECLARE_FUNC(int, ifft_cmplx, complex_t* COMMA int COMMA fft_t* COMMA complex_t* ); 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*); @@ -221,13 +223,15 @@ DECLARE_FUNC(int, fft, double* COMMA int COMMA fft_t* DECLARE_FUNC(int, fft_cmplx, complex_t* COMMA int COMMA fft_t* COMMA complex_t* ); DECLARE_FUNC(int, fft_create, fft_t* COMMA int); DECLARE_FUNC(void, fft_free, fft_t*); -DECLARE_FUNC(int, fft_shift, double* COMMA int n COMMA double*); +DECLARE_FUNC(int, fft_shift, double* COMMA int COMMA double*); +DECLARE_FUNC(int, fft_shift_cmplx, complex_t* COMMA int COMMA complex_t*); 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, freqs, double* COMMA double* COMMA int COMMA double* COMMA int COMMA complex_t*); DECLARE_FUNC(int, freqs_resp, double* COMMA double* COMMA int COMMA double* COMMA int COMMA int COMMA double* COMMA double* COMMA double*); +DECLARE_FUNC(int, freqs2time, double* COMMA double* COMMA int COMMA double COMMA int COMMA fft_t* COMMA double* COMMA double*); DECLARE_FUNC(int, freqz, double* COMMA double* 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*); diff --git a/release/include/dspl.c b/release/include/dspl.c index ef2f307..c01180c 100644 --- a/release/include/dspl.c +++ b/release/include/dspl.c @@ -51,6 +51,8 @@ p_cos_cmplx cos_cmplx ; p_dft dft ; p_dft_cmplx dft_cmplx ; p_dmod dmod ; +p_idft_cmplx idft_cmplx ; +p_ifft_cmplx ifft_cmplx ; p_farrow_lagrange farrow_lagrange ; p_farrow_spline farrow_spline ; p_filter_iir filter_iir ; @@ -60,12 +62,14 @@ p_fft_cmplx fft_cmplx ; p_fft_create fft_create ; p_fft_free fft_free ; p_fft_shift fft_shift ; +p_fft_shift_cmplx fft_shift_cmplx ; p_flipip flipip ; p_flipip_cmplx flipip_cmplx ; p_fourier_series_dec fourier_series_dec ; p_fourier_series_rec fourier_series_rec ; p_freqs freqs ; p_freqs_resp freqs_resp ; +p_freqs2time freqs2time ; p_freqz freqz ; p_goertzel goertzel ; p_goertzel_cmplx goertzel_cmplx ; @@ -153,6 +157,8 @@ void* dspl_load() LOAD_FUNC(dft); LOAD_FUNC(dft_cmplx); LOAD_FUNC(dmod); + LOAD_FUNC(idft_cmplx); + LOAD_FUNC(ifft_cmplx); LOAD_FUNC(farrow_lagrange); LOAD_FUNC(farrow_spline); LOAD_FUNC(filter_iir); @@ -162,6 +168,7 @@ void* dspl_load() LOAD_FUNC(fft_create); LOAD_FUNC(fft_free); LOAD_FUNC(fft_shift); + LOAD_FUNC(fft_shift_cmplx); LOAD_FUNC(flipip); LOAD_FUNC(flipip_cmplx); LOAD_FUNC(fourier_series_dec); @@ -169,6 +176,7 @@ void* dspl_load() LOAD_FUNC(freqz); LOAD_FUNC(freqs); LOAD_FUNC(freqs_resp); + LOAD_FUNC(freqs2time); LOAD_FUNC(goertzel); LOAD_FUNC(goertzel_cmplx); LOAD_FUNC(linspace); diff --git a/release/include/dspl.h b/release/include/dspl.h index 1c0288e..88d6e9c 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -213,6 +213,8 @@ DECLARE_FUNC(int, cos_cmplx, complex_t* COMMA int COMMA comple 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, idft_cmplx, complex_t* COMMA int COMMA complex_t*); +DECLARE_FUNC(int, ifft_cmplx, complex_t* COMMA int COMMA fft_t* COMMA complex_t* ); 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*); @@ -221,13 +223,15 @@ DECLARE_FUNC(int, fft, double* COMMA int COMMA fft_t* DECLARE_FUNC(int, fft_cmplx, complex_t* COMMA int COMMA fft_t* COMMA complex_t* ); DECLARE_FUNC(int, fft_create, fft_t* COMMA int); DECLARE_FUNC(void, fft_free, fft_t*); -DECLARE_FUNC(int, fft_shift, double* COMMA int n COMMA double*); +DECLARE_FUNC(int, fft_shift, double* COMMA int COMMA double*); +DECLARE_FUNC(int, fft_shift_cmplx, complex_t* COMMA int COMMA complex_t*); 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, freqs, double* COMMA double* COMMA int COMMA double* COMMA int COMMA complex_t*); DECLARE_FUNC(int, freqs_resp, double* COMMA double* COMMA int COMMA double* COMMA int COMMA int COMMA double* COMMA double* COMMA double*); +DECLARE_FUNC(int, freqs2time, double* COMMA double* COMMA int COMMA double COMMA int COMMA fft_t* COMMA double* COMMA double*); DECLARE_FUNC(int, freqz, double* COMMA double* 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*); diff --git a/test/src/butter_ap_test.c b/test/src/butter_ap_test.c index fd888b9..513585a 100644 --- a/test/src/butter_ap_test.c +++ b/test/src/butter_ap_test.c @@ -3,7 +3,7 @@ #include #include "dspl.h" -#define ORD 4 +#define ORD 3 #define N 1000 @@ -13,7 +13,7 @@ int main() handle = dspl_load(); // Load DSPL function double a[ORD+1], b[ORD+1]; - double Rp = 3.0; + double Rp = 1.0; double w[N], mag[N], phi[N], tau[N];