From 6b28f9391089bdad8588115ea2d8cca8eab1cb6b Mon Sep 17 00:00:00 2001 From: Dsplib Date: Tue, 18 Aug 2020 20:38:42 +0300 Subject: [PATCH] added function for analog and digital filters group delay calculation. Changes to be committed: modified: dspl/src/filter_an.c modified: dspl/src/gnuplot.c modified: examples/src/iir_bstop.c modified: examples/src/writetxt_3d_test.c modified: include/dspl.c modified: include/dspl.h --- dspl/src/filter_an.c | 139 +++++++++++++++++++++++--------- dspl/src/gnuplot.c | 95 +++++++++++++++++++++- examples/src/iir_bstop.c | 4 +- examples/src/writetxt_3d_test.c | 25 +++--- include/dspl.c | 4 +- include/dspl.h | 2 + 6 files changed, 213 insertions(+), 56 deletions(-) diff --git a/dspl/src/filter_an.c b/dspl/src/filter_an.c index 7837604..fa539ec 100644 --- a/dspl/src/filter_an.c +++ b/dspl/src/filter_an.c @@ -25,8 +25,104 @@ #include "dspl.h" +int DSPL_API group_delay(double* b, double* a, int ord, int flag, + double* w, int n, double* tau) +{ + double r, i, dr, di, br, bi, ar, ai, dbr, dbi, dar, dai, ar2ai2, ar2ai22; + int t, m; + + double *pa = NULL; + + if(!b || !w || !tau || (!a && (flag & DSPL_FLAG_ANALOG))) + return ERROR_PTR; + if(ord < 1) + return ERROR_FILTER_ORD; + if(n < 1) + return ERROR_SIZE; + if(a) + pa = a; + else + { + pa = (double*)malloc((ord+1) * sizeof(double)); + memset(pa, 0, (ord+1) * sizeof(double)); + pa[0] = 1.0; + } + + for(t = 0; t < n; t++) + { + br = bi = ar = ai = dbr = dbi = dar = dai = 0.0; + if(flag & DSPL_FLAG_ANALOG) + { + /* Br, Ar, Bi, Ai, Br', Ar', Bi', Ai' for analog filter */ + for(m = 0; m < ord+1; m+=4) + { + br += b[m] * pow(w[t], (double)m); + ar += pa[m] * pow(w[t], (double)m); + dbr += b[m] * (double) m * pow(w[t], (double)(m-1)); + dar += pa[m] * (double) m * pow(w[t], (double)(m-1)); + } + for(m = 2; m < ord+1; m+=4) + { + br -= b[m] * pow(w[t], (double)m); + ar -= pa[m] * pow(w[t], (double)m); + dbr -= b[m] * (double) m * pow(w[t], (double)(m-1)); + dar -= pa[m] * (double) m * pow(w[t], (double)(m-1)); + } + + for(m = 1; m < ord+1; m+=4) + { + bi += b[m] * pow(w[t], (double)m) ; + ai += pa[m] * pow(w[t], (double)m) ; + dbi += b[m] * (double) m * pow(w[t], (double)(m-1)) ; + dai += pa[m] * (double) m * pow(w[t], (double)(m-1)) ; + } + + for(m = 3; m < ord+1; m+=4) + { + bi -= b[m] * pow(w[t], (double)m) ; + ai -= pa[m] * pow(w[t], (double)m) ; + dbi -= b[m] * (double) m * pow(w[t], (double)(m-1)) ; + dai -= pa[m] * (double) m * pow(w[t], (double)(m-1)) ; + } + } + else + { + /* Br, Ar, Bi, Ai, Br', Ar', Bi', Ai' for digital filter */ + for(m = 0; m < ord+1; m++) + { + br += b[m] * cos(w[t]*(double)m); + bi -= b[m] * sin(w[t]*(double)m); + ar += pa[m] * cos(w[t]*(double)m); + ai -= pa[m] * sin(w[t]*(double)m); + + dbr -= b[m] *(double)m * sin(w[t]*(double)m); + dbi -= b[m] *(double)m * cos(w[t]*(double)m); + dar -= pa[m] *(double)m * sin(w[t]*(double)m); + dai -= pa[m] *(double)m * cos(w[t]*(double)m); + } + } + ar2ai2 = (ar * ar + ai * ai); + ar2ai22 = ar2ai2 * ar2ai2; + r = (br * ar + bi * ai) / ar2ai2; + i = (br * ai - bi * ar) / ar2ai2; + + dr = ((dbr * ar + dar * br + dbi * ai + dai * bi) * ar2ai2 - + (2.0 * ar * dar + 2.0 * ai * dai) * (br * ar + bi * ai))/ar2ai22; + + di = ((dbr * ai + dai * br - dbi * ar - dar * bi) * ar2ai2 - + (2.0 * ar * dar + 2.0 * ai * dai) * (br * ai - bi * ar))/ar2ai22; + + tau[t] = -(dr * i - di * r) / (r*r + i*i); + } + + if(pa != a) + free(pa); + + return RES_OK; +} + #ifdef DOXYGEN_ENGLISH /*! **************************************************************************** @@ -241,9 +337,9 @@ int DSPL_API filter_freq_resp(double* b, double* a, int ord, { int res, k, flag_analog; - complex_t *hc = NULL; - double *phi0 = NULL; - double *phi1 = NULL; + complex_t *hc = NULL; + double *phi0 = NULL; + double *phi1 = NULL; double *w0 = NULL; double *w1 = NULL; @@ -296,34 +392,8 @@ int DSPL_API filter_freq_resp(double* b, double* a, int ord, if(tau) - { - phi0 = (double*) malloc(n*sizeof(double)); - phi1 = (double*) malloc(n*sizeof(double)); - w0 = (double*) malloc(n*sizeof(double)); - w1 = (double*) malloc(n*sizeof(double)); + res = group_delay(b, a, ord, flag, w, n, tau); - w0[0] = w[0] - (w[1] - w[0])*0.02; - w1[0] = w[0] + (w[1] - w[0])*0.02; - - for(k = 1; k < n; k++) - { - w0[k] = w[k] - (w[k] - w[k-1])*0.02; - w1[k] = w[k] + (w[k] - w[k-1])*0.02; - } - res = filter_freq_resp(b, a, ord, w0, n, - DSPL_FLAG_UNWRAP | flag_analog, - NULL, phi0, NULL); - if(res != RES_OK) - goto exit_label; - res = filter_freq_resp(b, a, ord, w1, n, - DSPL_FLAG_UNWRAP | flag_analog, - NULL, phi1, NULL); - if(res != RES_OK) - goto exit_label; - - for(k = 0; k < n; k++) - tau[k] = (phi0[k] - phi1[k])/(w1[k] - w0[k]); - } exit_label: @@ -585,11 +655,6 @@ int DSPL_API freqs_cmplx(double* b, double* a, int ord, } - - - - - #ifdef DOXYGEN_ENGLISH #endif @@ -734,7 +799,6 @@ H(z) = \frac {\sum_{k = 0}^{N} b_k z^{-k}} Комплексный коэффициент передачи рассчитывается путем подстановки \f$z = e^{j \omega} \f$. \n - \param[in] b Указатель на вектор коэффициентов числителя передаточной функции \f$H(z)\f$. \n @@ -754,7 +818,7 @@ H(z) = \frac {\sum_{k = 0}^{N} b_k z^{-k}} для которого будет рассчитан комплексный коэффициент передачи \f$ H \left(e^{j \omega} \right)\f$. \n Размер вектора `[n x 1]`. \n \n - + \param[in] n Размер вектора нормированной циклической частоты `w`. \n \n @@ -765,7 +829,6 @@ H(z) = \frac {\sum_{k = 0}^{N} b_k z^{-k}} Размер вектора `[n x 1]`. \n Память должна быть выделена. \n \n - \return `RES_OK` Комплексный коэффициент передачи рассчитан успешно. \n В противном случае \ref ERROR_CODE_GROUP "код ошибки". \n diff --git a/dspl/src/gnuplot.c b/dspl/src/gnuplot.c index be014d9..e131269 100644 --- a/dspl/src/gnuplot.c +++ b/dspl/src/gnuplot.c @@ -22,8 +22,8 @@ #include "dspl.h" #define GNUPLOT_NO 1 -#define GNUPLOT_WIN 2 -#define GNUPLOT_PNG 3 +#define GNUPLOT_WIN 3 +#define GNUPLOT_PNG 4 @@ -219,6 +219,8 @@ int DSPL_API gnuplot_create(int argc, char* argv[], if(!strcmp(argv[1], "--plotpng")) state = GNUPLOT_PNG; } + if(!hplot) + return ERROR_PTR; switch(state) { @@ -267,6 +269,7 @@ int DSPL_API gnuplot_create(int argc, char* argv[], + #ifdef DOXYGEN_ENGLISH /*! **************************************************************************** \ingroup PLOT_GROUP @@ -400,5 +403,93 @@ void DSPL_API gnuplot_cmd(void* h, char* cmd) +#ifdef DOXYGEN_ENGLISH +/*! **************************************************************************** +\ingroup PLOT_GROUP +\fn int gnuplot_open(void** hplot) +\brief Open GNUPLOT program. + +This function opens the GNUPLOT package. +After calling this function, the handle of the GNUPLOT +will be written to the address `hplot` and it becomes possible to send GNUPLOT +commands. + +\note From a system point of view, `hplot` is a pointer to an open file +in which you can write commands for execution by the GNUPLOT package. + +\param[in, out] hplot +Pointer to the handle address of the GNUPLOT package. \n +This pointer is required to send GNUPLOT commands. \n +\n + +\return +`RES_OK` if function is calculated successfully. \n +Else \ref ERROR_CODE_GROUP "code error". +The `hplot` pointer sets in `NULL` if function returns error. \n + +GNUPLOT handle must be closed by \ref gnuplot_close after plotting.\n + +An example of plotting sine and cosine is given in the following listing: + +\note +The difference between `gnuplot_open` and` gnuplot_create` is that +`gnuplot_create` processes the program execution parameters and +creates a GNUPLOT terminal. \n +The gnuplot_open function opens the GNUPLOT handle to be able to +send commands, regardless of program execution parameters, +but does not create terminals. + +\author Sergey Bakhurin www.dsplib.org +***************************************************************************** */ +#endif +#ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup PLOT_GROUP +\fn int gnuplot_open(void** hplot) +\brief Открыть пакет GNUPLOT. + +Данная функция открывает пакет GNUPLOT. +После вызова данной функции по адресу `hplot` будет записан +handle GNUPLOT и появляется возможность посылать GNUPLOT команды. + +\note С точки зрения системы, `hplot` является указателем на открытый файл, + в который можно записывать команды для исполнения пакетом GNUPLOT. + +\param[in, out] hplot +Указатель на адрес хэндла пакета GNUPLOT. \n +По данному адресу будет записан указатель на текщий график. Данный указатель +необходим для посылки команд GNUPLOT для построения графика. \n +\n + +\return +`RES_OK` --- функция выполнена успешно. \n +В противном случае \ref ERROR_CODE_GROUP "код ошибки". \n +При возникновении ошибки по адресу `hplot` будет записан `NULL`. \n + +После построения графиков необходимо закрыть хэндл GNUPLOT функцией +\ref gnuplot_close. \n + +\note +Отличие функции `gnuplot_open` от `gnuplot_create` заключается в том, +что `gnuplot_create` обрабатывает параметры выполнения программы +и создает терминал GNUPLOT.\n +Функция `gnuplot_open` открывает GNUPLOT хэндл для возможности посылки команд, +вне зависимости от параметров выполнения программы, но не создает терминалов. + +\author Бахурин Сергей www.dsplib.org +***************************************************************************** */ +#endif +int DSPL_API gnuplot_open(void** hplot) +{ + if(!hplot) + return ERROR_PTR; + *hplot = popen("gnuplot -p", "w"); + + if(!(*hplot)) + return ERROR_GNUPLOT_CREATE; + return RES_OK; +} + + diff --git a/examples/src/iir_bstop.c b/examples/src/iir_bstop.c index 9fe371b..3ee2bdd 100644 --- a/examples/src/iir_bstop.c +++ b/examples/src/iir_bstop.c @@ -7,7 +7,7 @@ #define ORD 14 // Frequency response vector size -#define N 1001 +#define N 5001 int main(int argc, char* argv[]) @@ -27,7 +27,7 @@ int main(int argc, char* argv[]) // Calculate Chebyshev type 2 digital bandstop filter int res = iir(rp, rs, ORD, 0.3, 0.7, - DSPL_FILTER_CHEBY2 | DSPL_FILTER_BSTOP, b, a); + DSPL_FILTER_ELLIP | DSPL_FILTER_BSTOP, b, a); if(res != RES_OK) printf("error code = 0x%8x\n", res); diff --git a/examples/src/writetxt_3d_test.c b/examples/src/writetxt_3d_test.c index b8f9c7a..fed3f96 100644 --- a/examples/src/writetxt_3d_test.c +++ b/examples/src/writetxt_3d_test.c @@ -2,31 +2,31 @@ #include #include #include - + #include "dspl.h" - + #define NX 20 #define NY 30 - + int main(int argc, char* argv[]) { void* hdspl; /* DSPL handle */ void* hplot; /* GNUPLOT handles */ - + hdspl = dspl_load(); /* Load DSPL function */ - + double x[NX]; double y[NY]; double z[NX * NY]; int n, m; int err; - + /* x vector from -2 to 2 */ linspace(-2.0, 2.0, NX, DSPL_SYMMETRIC, x); - + /* y vector from -2.5 to 2.5 */ linspace(-2.5, 2.5, NY, DSPL_SYMMETRIC, y); - + /* z(x,y) = x * exp(-x^2 - y^2) */ for(n = 0; n < NX; n++) { @@ -35,11 +35,11 @@ int main(int argc, char* argv[]) z[n + m*NX] = x[n]*exp(-x[n]*x[n] - y[m]*y[m]); } } - + /* Save to files "dat/data3d.txt" */ err = writetxt_3d(x, NX, y, NY, z, "dat/data3d.txt"); printf("writetxt_3d error 0x%8x\n", err); - + /* plotting 3d surface by GNUPLOT */ /* Create window 0 */ err = gnuplot_create(argc, argv, 560, 480, "img/writetxt_3d.png", &hplot); @@ -50,8 +50,7 @@ int main(int argc, char* argv[]) gnuplot_cmd(hplot, "set ylabel 'y'"); gnuplot_cmd(hplot, "splot 'dat/data3d.txt' with lines"); gnuplot_close(hplot); - + dspl_free(hdspl); /* free dspl handle */ return 0; -} - +} \ No newline at end of file diff --git a/include/dspl.c b/include/dspl.c index 117e57b..ba0ab98 100644 --- a/include/dspl.c +++ b/include/dspl.c @@ -110,6 +110,7 @@ p_freqz freqz ; p_gnuplot_close gnuplot_close ; p_gnuplot_cmd gnuplot_cmd ; p_gnuplot_create gnuplot_create ; +p_gnuplot_open gnuplot_open ; p_goertzel goertzel ; p_goertzel_cmplx goertzel_cmplx ; @@ -409,7 +410,8 @@ void* dspl_load() LOAD_FUNC(gnuplot_close); LOAD_FUNC(gnuplot_cmd); - LOAD_FUNC(gnuplot_create); + LOAD_FUNC(gnuplot_create); + LOAD_FUNC(gnuplot_open); LOAD_FUNC(goertzel); LOAD_FUNC(goertzel_cmplx); diff --git a/include/dspl.h b/include/dspl.h index 1f0a1fe..55b1901 100644 --- a/include/dspl.h +++ b/include/dspl.h @@ -1055,6 +1055,8 @@ DECLARE_FUNC(int, gnuplot_create, int argc COMMA char* fn_png COMMA void** hplot); /*----------------------------------------------------------------------------*/ +DECLARE_FUNC(int, gnuplot_open, void** hplot); +/*----------------------------------------------------------------------------*/ DECLARE_FUNC(int, goertzel, double* COMMA int COMMA int*