kopia lustrzana https://github.com/Dsplib/libdspl-2.0
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.hpull/6/merge
Dsplib
rodzic
a8a0b4a0c1
commit
6b28f93910
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
||||
|
|
|
|||
|
|
@ -2,31 +2,31 @@
|
|||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
|
||||
|
||||
#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;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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);
|
||||
|
||||
|
|
|
|||
|
|
@ -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*
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue