diff --git a/dspl/src/conv.c b/dspl/src/conv.c index 5a687bb..804672a 100644 --- a/dspl/src/conv.c +++ b/dspl/src/conv.c @@ -993,15 +993,16 @@ dat/sf.txt - сигнал на выходе фильтра. int DSPL_API filter_iir(double* b, double* a, int ord, double* x, int n, double* y) { - double* buf = NULL; - double* an = NULL; + double *buf = NULL; + double *an = NULL; + double *bn = NULL; double u; int k; int m; int count; if(!b || !x || !y) - return ERROR_PTR; + return ERROR_PTR; if(ord < 1 || n < 1) return ERROR_SIZE; @@ -1016,10 +1017,19 @@ int DSPL_API filter_iir(double* b, double* a, int ord, memset(buf, 0, count*sizeof(double)); if(!a) + { memset(an, 0, count*sizeof(double)); + bn = b; + } else + { + bn = (double*) malloc(count*sizeof(double)); for(k = 0; k < count; k++) + { an[k] = a[k] / a[0]; + bn[k] = b[k] / a[0]; + } + } for(k = 0; k < n; k++) { @@ -1032,13 +1042,15 @@ int DSPL_API filter_iir(double* b, double* a, int ord, buf[0] = x[k] - u; y[k] = 0.0; for(m = 0; m < count; m++) - y[k] += buf[m] * b[m]; + y[k] += buf[m] * bn[m]; } if(buf) free(buf); if(an) free(an); + if(bn && (bn != b)) + free(bn); return RES_OK; } diff --git a/dspl/src/filter_an.c b/dspl/src/filter_an.c index ffaebe3..d2700f9 100644 --- a/dspl/src/filter_an.c +++ b/dspl/src/filter_an.c @@ -31,7 +31,7 @@ #ifdef DOXYGEN_ENGLISH /*! **************************************************************************** \ingroup FILTER_ANALYSIS_GROUP -int DSPL_API group_delay(double* b, double* a, int ord, int flag, +\fn int DSPL_API group_delay(double* b, double* a, int ord, int flag, double* w, int n, double* tau) \brief @@ -40,7 +40,7 @@ Group delay calculation for digital or analog filter corresponds to Group delay is describes as: \f[ -\tau_g(\omega) = - \dfrac{d\Phi(\omega)}{d\omega}, +\tau_g(\omega) = - \frac{d\Phi(\omega)}{d\omega}, \f] here \f$\Phi(\omega)\f$ -- filter phase response, \f$\omega\f$ is angular frequency for analog filter, or normalized frequency for digital filter. @@ -98,14 +98,14 @@ Else \ref ERROR_CODE_GROUP "code error". #ifdef DOXYGEN_RUSSIAN /*! **************************************************************************** \ingroup FILTER_ANALYSIS_GROUP -int DSPL_API group_delay(double* b, double* a, int ord, int flag, +\fn int DSPL_API group_delay(double* b, double* a, int ord, int flag, double* w, int n, double* tau) \brief Расчет группового времени запаздывания цифрового или аналогового фильтра. Групповое время запаздывания определяется как: \f[ -\tau_g(\omega) = - \dfrac{d\Phi(\omega)}{d\omega}, +\tau_g(\omega) = - \frac{d\Phi(\omega)}{d\omega}, \f] где \f$\Phi(\omega)\f$ -- ФЧХ фильтра, \f$\omega\f$ циктическая частот в случае аналогового фильтра, или нормированная частота цифрового фильтра. @@ -157,7 +157,7 @@ DSPL_FLAG_DIGITAL Коэффициенты относятся к цифрово Память должна быть выделена. \n \return -`RES_OK` Параметры фильтра рассчитаны успешно. \n +`RES_OK` групповая задержка фильтра рассчитана успешно. \n В противном случае \ref ERROR_CODE_GROUP "код ошибки". \n \author Бахурин Сергей www.dsplib.org @@ -1053,10 +1053,139 @@ exit_label: #ifdef DOXYGEN_ENGLISH +/*! **************************************************************************** +\ingroup FILTER_ANALYSIS_GROUP +\fn int DSPL_API phase_delay(double* b, double* a, int ord, int flag, + double* w, int n, double* tau) +\brief +Phase delay calculation for digital or analog filter corresponds to +\f$H(s)\f$, or \f$H(z)\f$ transfer function. + +Group delay is describes as: +\f[ +\tau_\{\varphi}(\omega) = - \frac{\Phi(\omega)}{\omega}, +\f] +here \f$\Phi(\omega)\f$ -- filter phase response, \f$\omega\f$ is angular +frequency for analog filter, or normalized frequency for digital filter. + +\param[in] b +Pointer to the \f$ H(s) \f$ or \f$H(z)\f$ transfer function +numerator coefficients vector. \n +Vector size is `[ord+1 x 1]`. \n \n + +\param[in] a +Pointer to the \f$ H(s) \f$ or \f$H(z)\f$ transfer function +denominator coefficients vector. \n +Vector size is `[ord+1 x 1]`. \n \n + +\param[in] ord +Filter order. \n +Transfer function \f$ H(s) \f$ or \f$H(z)\f$ numerator +and denominator coefficients number equals `ord+1`. \n \n + +\param[in] flag +Binary flags to set calculation rules: \n +\verbatim +DSPL_FLAG_ANALOG Coefficients corresponds to analog filter +DSPL_FLAG_DIGITAL Coefficients corresponds to digital filter +\endverbatim +\n \n + +\param[in] w +Pointer to the angular frequency \f$ \omega \f$ (rad/s), +which used for analog filter characteristics calculation +(flag sets as `DSPL_FLAG_ANALOG`). \n +For digital filter (flag sets as `DSPL_FLAG_DIGITAL`), + parameter `w` describes normalized frequency of +frequency response \f$ H \left(\mathrm{e}^{j\omega} \right) \f$. +Digital filter frequency response is \f$ 2\pi \f$-periodic function, +and vector `w` advisable to set from 0 to \f$ \pi \f$, +or from 0 to \f$ 2\pi \f$, or from \f$ -\pi \f$ to \f$ \pi \f$. +Vector size is `[n x 1]`. \n \n + +\param[in] n +Size of frequency vector `w`. \n \n + +\param[out] tau +Pointer to the phase delay vector. \n +Vector size is `[n x 1]`. \n +Memory must be allocated. \n \n + +\return +\return `RES_OK` if function is calculated successfully. \n +Else \ref ERROR_CODE_GROUP "code error". + +\author Sergey Bakhurin www.dsplib.org +***************************************************************************** */ #endif #ifdef DOXYGEN_RUSSIAN +/*! **************************************************************************** +\ingroup FILTER_ANALYSIS_GROUP +\fn int DSPL_API phase_delay(double* b, double* a, int ord, int flag, + double* w, int n, double* tau) +\brief +Расчет задержки цифрового или аналогового фильтра. +Фазовая задержка определяется как: +\f[ +\tau_g(\omega) = - \frac{\Phi(\omega)}{\omega}, +\f] +где \f$\Phi(\omega)\f$ -- ФЧХ фильтра, \f$\omega\f$ циктическая частот в случае +аналогового фильтра, или нормированная частота цифрового фильтра. + +\param[in] b +Указатель на вектор коэффициентов числителя передаточной функции +аналогового фильтра \f$ H(s) \f$ или цифрового фильтра \f$ H(z) \f$. \n +Размер вектора `[ord+1 x 1]`. \n \n + +\param[in] a +Указатель на вектор коэффициентов числителя передаточной функции +аналогового фильтра \f$ H(s) \f$ или цифрового фильтра \f$ H(z) \f$. \n +Размер вектора `[ord+1 x 1]`. \n +Параметр может быть `NULL`. В этом случае расчет производится для цифрового +КИХ-фильтра с коэффициентами, заданными вектором `b`. \n\n + +\param[in] ord +Порядок фильтра. Количество коэффициентов +числителя и знаменателя передаточной +функции \f$ H(s) \f$ или \f$ H(z) \f$ равно `ord+1`. \n \n + +\param[in] flag +Флаг который задает тип фильтра: \n +\verbatim +DSPL_FLAG_ANALOG Коэффициенты относятся к аналоговому фильтру +DSPL_FLAG_DIGITAL Коэффициенты относятся к цифровому фильтру +\endverbatim + +\param[in] w +Указатель на вектор значений циклической частоты \f$ \omega \f$ (рад/с), +для которого будет рассчитаны АЧХ, ФЧХ и ГВЗ аналогового фильтра, +если установлен флаг `DSPL_FLAG_ANALOG`. \n +В случае если флаг `DSPL_FLAG_ANALOG` не установлен, то вектор частоты `w` +используется как нормированная частота комплексного коэффициента передачи +\f$ H \left(\mathrm{e}^{j\omega} \right) \f$ цифрового фильтра. \n +В этом случае характеристика цифрового фильтра является +\f$ 2\pi \f$-периодической, и вектор частоты может содержать +произвольные значения, однако целесообразно задавать +его от 0 до \f$ \pi \f$, а такжет от 0 до \f$ 2\pi \f$, или +от \f$ -\pi \f$ до \f$ \pi \f$. \n +Размер вектора `[n x 1]`. \n \n + +\param[in] n +Размер вектора циклической частоты `w`. \n \n + +\param[out] tau +Указатель на вектор фазовой задержки. \n +Размер вектора `[n x 1]`. \n +Память должна быть выделена. \n + +\return +`RES_OK` фазовая задержка фильтра рассчитана успешно. \n +В противном случае \ref ERROR_CODE_GROUP "код ошибки". \n + +\author Бахурин Сергей www.dsplib.org +***************************************************************************** */ #endif int DSPL_API phase_delay(double* b, double* a, int ord, int flag, double* w, int n, double* tau) @@ -1069,13 +1198,13 @@ int DSPL_API phase_delay(double* b, double* a, int ord, int flag, else return ERROR_SIZE; - err = filter_freq_resp(b, a, ord, w, n, - flag | DSPL_FLAG_UNWRAP, NULL, phi, NULL); + err = filter_freq_resp(b, a, ord, w, n, flag | DSPL_FLAG_UNWRAP, NULL, phi, NULL); if(err!=RES_OK) goto exit_label; for(i = 0; i < n; i++) + { tau[i] = w[i] ? ( - phi[i] / w[i]) : ( - phi[i] / (w[i] + 1E-9) ); - + } exit_label: if(phi) free(phi); diff --git a/dspl/src/filter_iir.c b/dspl/src/filter_iir.c index f7b0dd9..e1c4dc5 100644 --- a/dspl/src/filter_iir.c +++ b/dspl/src/filter_iir.c @@ -410,7 +410,8 @@ int DSPL_API iir(double rp, double rs, int ord, double w0, double w1, double *at = NULL; double wa0, wa1, ws; int err, ord_ap = ord; - + int i; + if(((type & DSPL_FILTER_TYPE_MASK) == DSPL_FILTER_LPF) || ((type & DSPL_FILTER_TYPE_MASK) == DSPL_FILTER_HPF)) { @@ -477,6 +478,14 @@ int DSPL_API iir(double rp, double rs, int ord, double w0, double w1, err = bilinear(bt, at, ord, b, a); + + for(i = 1; i <= ord; i++) + { + a[i] /= a[0]; + b[i] /= a[0]; + } + b[0] /= a[0]; + a[0] = 1.0; error_proc: