libdspl-2.0
Библиотека алгоритмов цифровой обработки сигналов
Расчет БИХ-фильтров

Функции

int butter_ap (double rp, int ord, double *b, double *a)
 Расчет передаточной характеристики $ H(s) $ аналогового нормированного ФНЧ Баттерворта. Подробнее...
 
int cheby1_ap (double rp, int ord, double *b, double *a)
 Расчет передаточной характеристики $ H(s) $ аналогового нормированного ФНЧ Чебышёва первого рода. Подробнее...
 
int cheby2_ap (double rs, int ord, double *b, double *a)
 Расчет передаточной характеристики $ H(s) $ аналогового нормированного ФНЧ Чебышёва второго рода. Подробнее...
 
int low2high (double *b, double *a, int ord, double w0, double w1, double *beta, double *alpha)
 Частотное преобразование ФНЧ-ФВЧ Подробнее...
 
int low2low (double *b, double *a, int ord, double w0, double w1, double *beta, double *alpha)
 Частотное преобразование ФНЧ-ФНЧ Подробнее...
 
int ratcompos (double *b, double *a, int n, double *c, double *d, int p, double *beta, double *alpha)
 Рациональная композиця Подробнее...
 
int bilinear (double *bs, double *as, int ord, double *bz, double *az)
 Билинейное преобразование передаточной характеристики аналогового фильтра $H(s)$, в передаточную характеристику цифрового фильтра $H(z)$. Подробнее...
 

Подробное описание

Функции расчета цифровых БИХ-фильтров.

Функции

◆ bilinear()

int bilinear ( double *  bs,
double *  as,
int  ord,
double *  bz,
double *  az 
)

Билинейное преобразование передаточной характеристики аналогового фильтра $H(s)$, в передаточную характеристику цифрового фильтра $H(z)$.


Функция рассчитывает коэффициенты передаточной характеристики $H(z)$ цифрового фильтра путем дробно-рациональной подстановки вида

\[ s \leftarrow \frac{1 - z^{-1}}{1 - z^{-1}}. \]

Порядок цифрового фильтра при этом остается равным порядку аналогового фильтра, а ось частот $\Omega$ аналогового фильтра связана c осью частот $\omega$ цифрового фильтра соотношением:

\[ \Omega = \tan(\omega / 2). \]

Аргументы
[in]bsУказатель на вектор коэффициентов числителя передаточной функции $H(s)$ исходного аналогового фильтра.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[in]asУказатель на вектор коэффициентов знаменателя передаточной функции $H(s)$ исходного аналогового фильтра.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[in]ordПорядок фильтра.
Количество коэффициентов числителя и знаменателя передаточных функций $H(s)$ и $H(z)$ аналогового и цифрового фильтров равно ord+1.

[out]bzУказатель на вектор коэффициентов числителя передаточной функции $H(z)$ полученного цифрового фильтра.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[out]azУказатель на вектор коэффициентов знаменателя передаточной функции $H(z)$ полученного цифрового фильтра.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

Возвращает
RES_OK Фильтр рассчитан успешно.

В противном случае код ошибки.

Пример использования функции bilinear

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define N 1000
#define ORD 4
int main(int argc, char* argv[])
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
double w[N], h[N];
complex_t hz[N];
double bs[ORD+1], as[ORD+1];
double bz[ORD+1], az[ORD+1];
int err, k;
// расчет аналогового фильтра Чебышева первого рода
err = cheby1_ap(1.0, ORD, bs, as);
if(err != RES_OK)
{
printf("cheby1_ap error code %d\n", err);
return err;
}
// Билинейное преобразование
err = bilinear(bs, as, ORD, bz, az);
if(err != RES_OK)
{
printf("bilinear error code %d\n", err);
return err;
}
// Печать коэффициентов
for(k = 0; k < ORD+1; k++)
printf("bz[%d] = %7.3f az[%d] = %7.3f\n", k, bz[k], k, az[k]);
// Расчет АЧХ полученного цифрового фильтра
linspace(0, M_PI, N, DSPL_PERIODIC, w);
freqz(bz, az, ORD, w, N, hz);
for(k = 0; k < N; k++)
{
h[k] = 10.0 * log10 (ABSSQR(hz[k])); // АЧХ в дБ
w[k] /= M_PI; // нормировка частоты от 0 до 1
}
writetxt(w,h,N,"dat/bilinear.txt");
/* run GNUPLOT script */
err = gnuplot_script(argc, argv, "gnuplot/bilinear_test.plt");
dspl_free(handle); // free dspl handle
return err;
}

Данная программа производит расчет передаточной характеристики аналогового фильтра Чебышева первого рода, с частотой среза равной 1 рад/с, и производит билинейное преобразование в цифровой, с частотой среза равной 0.5.

Результат работы программы:

bz[0] =   0.246    az[0] =   4.425
bz[1] =   0.983    az[1] =  -3.318
bz[2] =   1.474    az[2] =   4.746
bz[3] =   0.983    az[3] =  -2.477
bz[4] =   0.246    az[4] =   1.034
err = 0

Кроме этого производится расчет АЧХ полученного цифрового фильтра и строится график АЧХ пакетом GNUPLOT

bilinear.png

Скрипт GNUPLOT для построения графиков из текстовых файлов:

if(!exists("plotterm")) plotterm = 'wxt'
if(plotterm eq "pngcairo") set output 'img/bilinear.png'
unset key
set grid
set xlabel "normalized frequency"
set terminal plotterm size 520, 380 enhanced font 'Verdana,8'
set ylabel "Magnitude, dB"
set yrange [-80:5]
plot 'dat/bilinear.txt' with lines
Автор
Бахурин Сергей www.dsplib.org

См. определение в файле filter_iir.c строка 80

◆ butter_ap()

int butter_ap ( double  Rp,
int  ord,
double *  b,
double *  a 
)

Расчет передаточной характеристики $ H(s) $ аналогового нормированного ФНЧ Баттерворта.


Функция рассчитывает коэффициенты передаточной характеристики $H(s)$ аналогового нормированного ФНЧ Баттерворта порядка ord с частотой среза 1 рад/с по уровню $ -R_p $ дБ.

Аргументы
[in]RpНеравномерность АЧХ в полосе пропускания (дБ).
Параметр задает уровень искажений в полосе от 0 до 1 рад/с.
Значение должно быть положительным.

[in]ordПорядок фильтра.
Количество коэффициентов числителя и знаменателя передаточной функции $H(s)$ равно ord+1.

[out]bУказатель на вектор коэффициентов числителя передаточной функции $H(s)$.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[out]aУказатель на вектор коэффициентов знаменателя передаточной функции $H(s)$.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

Возвращает
RES_OK Фильтр рассчитан успешно.

В противном случае код ошибки.
Автор
Бахурин Сергей www.dsplib.org

См. определение в файле filter_ap.c строка 33

◆ cheby1_ap()

int cheby1_ap ( double  Rp,
int  ord,
double *  b,
double *  a 
)

Расчет передаточной характеристики $ H(s) $ аналогового нормированного ФНЧ Чебышёва первого рода.


Функция рассчитывает коэффициенты передаточной характеристики $ H(s)$ аналогового нормированного ФНЧ Чебышёва первого рода порядка ord с частотой среза 1 рад/с по уровню $-R_p$ дБ.

Особенностью фильтра Чебышёва первого рода являются равноволновые пульсации АЧХ в полосе пропускания.

Аргументы
[in]RpНеравномерность АЧХ в полосе пропускания (дБ).
Параметр задает уровень искажений в полосе от 0 до 1 рад/с.
Значение должно быть положительным.

[in]ordПорядок фильтра.
Количество коэффициентов числителя и знаменателя передаточной функции $ H(s)$ равно ord+1.

[out]bУказатель на вектор коэффициентов числителя передаточной функции $H(s)$.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[out]aУказатель на вектор коэффициентов знаменателя передаточной функции $H(s)$.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

Возвращает
RES_OK Фильтр рассчитан успешно.

В противном случае код ошибки.

Пример использования функции cheby1_ap:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
// Порядок фильтра
#define ORD 4
// размер векторов частотной характериситки фильтра
#define N 1000
int main(int argc, char* argv[])
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
double a[ORD+1], b[ORD+1]; // коэффицинеты H(s)
double Rp = 3.0; // неравномерность в полосе пропускания 3дБ
// Частота (w), АЧХ (mag), ФЧХ (phi) и ГВЗ (tau)
double w[N], mag[N], phi[N], tau[N];
int k;
// рассчитываем нормированный ФНЧ Чебышева 1 рода
int res = cheby1_ap(Rp, ORD, b, a);
if(res != RES_OK)
printf("error code = 0x%8x\n", res);
// печать коэффициентов фильтра
for(k = 0; k < ORD+1; k++)
printf("b[%2d] = %9.3f a[%2d] = %9.3f\n", k, b[k], k, a[k]);
// вектор частоты в логарифмической шакале от 0.01 до 100 рад/c
logspace(-2.0, 2.0, N , DSPL_SYMMETRIC, w);
//частотные характеристика фильтра
filter_freq_resp(b, a, ORD, w, N,
DSPL_FLAG_LOGMAG|DSPL_FLAG_UNWRAP | DSPL_FLAG_ANALOG,
mag, phi, tau);
// Сохранить характеристики для построения графиков
writetxt(w, mag, N, "dat/cheby1_ap_test_mag.txt");
writetxt(w, phi, N, "dat/cheby1_ap_test_phi.txt");
writetxt(w, tau, N, "dat/cheby1_ap_test_tau.txt");
/* run GNUPLOT script */
res = gnuplot_script(argc, argv, "gnuplot/cheby1_ap_test.plt");
dspl_free(handle); // free dspl handle
return res;
}

Результат работы программы:

b[ 0] =     0.125     a[ 0] =     0.177
b[ 1] =     0.000     a[ 1] =     0.405
b[ 2] =     0.000     a[ 2] =     1.169
b[ 3] =     0.000     a[ 3] =     0.582
b[ 4] =     0.000     a[ 4] =     1.000



В каталоге dat будут созданы три файла:

cheby1_ap_test_mag.txt    АЧХ фильтра   
cheby1_ap_test_phi.txt    ФЧХ фильтра
cheby1_ap_test_tau.txt    ГВЗ фильтра

Кроме того программа GNUPLOT произведет построение следующих графиков по сохраненным в файлах данным:

cheby1_ap_test.png

Скрипт GNUPLOT для построения графиков из текстовых файлов:

if(!exists("plotterm")) plotterm = 'wxt'
if(plotterm eq "pngcairo") set output 'img/cheby1_ap_test.png'
set logscale x
unset key
set grid
set xlabel "frequency, rad/s"
set terminal plotterm size 920, 260 enhanced font 'Verdana,8'
set multiplot layout 1,3 rowsfirst
set ylabel "Magnitude, dB"
set yrange [-100:5]
plot 'dat/cheby1_ap_test_mag.txt' with lines
set ylabel "Phase response, rad"
unset yrange
plot 'dat/cheby1_ap_test_phi.txt' with lines
set ylabel "Groupdelay, sec"
unset yrange
plot 'dat/cheby1_ap_test_tau.txt' with lines
unset multiplot
Автор
Бахурин Сергей www.dsplib.org

См. определение в файле filter_ap.c строка 129

◆ cheby2_ap()

int cheby2_ap ( double  Rs,
int  ord,
double *  b,
double *  a 
)

Расчет передаточной характеристики $ H(s) $ аналогового нормированного ФНЧ Чебышёва второго рода.


Функция рассчитывает коэффициенты передаточной характеристики $H(s)$ аналогового нормированного ФНЧ Чебышёва второго рода порядка ord с частотой заграждения 1 рад/с по уровню $-R_s$ дБ.
Особенностью фильтра Чебышёва второго рода являются:
1) равноволновые пульсации АЧХ в полосе заграждения.
2) уровень АЧХ $H(j\cdot 1) = -R_s$ дБ.

Аргументы
[in]RsУровень подавления в полосе пропускания (дБ).
Значение должно быть положительным.

[in]ordПорядок фильтра.
Количество коэффициентов числителя и знаменателя передаточной функции $H(s)$ равно ord+1.

[out]bУказатель на вектор коэффициентов числителя передаточной функции $H(s)$.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[out]aУказатель на вектор коэффициентов знаменателя передаточной функции $H(s)$.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

Возвращает
RES_OK Фильтр рассчитан успешно.
В противном случае код ошибки.
Автор
Бахурин Сергей www.dsplib.org

См. определение в файле filter_ap.c строка 239

◆ low2high()

int low2high ( double *  b,
double *  a,
int  ord,
double  w0,
double  w1,
double *  beta,
double *  alpha 
)

Частотное преобразование ФНЧ-ФВЧ


Функция производит перобразование передаточной функции $ H(s) $ аналогового ФНЧ с частотой среза w0 рад/c в передаточную функцию $ F(s) $ аналоговго ФВЧ с частотой среза w1 рад/c.

Неравномерность АЧХ в полосе пропускания, уровень подавления в полосе заграждения и порядок фильтра остаются неизменными.

Аргументы
[in]bУказатель на вектор коэффициентов числителя передаточной функции $H(s)$ исходного аналогового ФНЧ.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[in]aУказатель на вектор коэффициентов знаменателя передаточной функции $H(s)$ исходного аналогового ФНЧ.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[in]ordПорядок исходного фильтра и фильтра после переобразования.

[in]w0Частота среза исходного ФНЧ.

[in]w1Требуемая частота среза ФВЧ после преобразования.

[in,out]betaУказатель на вектор коэффициентов числителя передаточной функции $F(s)$ ФВЧ после преобразования.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[in,out]alphaУказатель на вектор коэффициентов знаменателя передаточной функции $F(s)$ аналогового ФВЧ после преобразования.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

Возвращает
RES_OK Преоборазование расчитано успешно.

В противном случае код ошибки.
Автор
Бахурин Сергей www.dsplib.org

См. определение в файле filter_ft.c строка 151

◆ low2low()

int low2low ( double *  b,
double *  a,
int  ord,
double  w0,
double  w1,
double *  beta,
double *  alpha 
)

Частотное преобразование ФНЧ-ФНЧ


Функция производит перобразование передаточной функции $ H(s) $ аналогового ФНЧ с частотой среза w0 рад/c в передаточную функцию $ F(s) $ аналоговго ФНЧ с частотой среза w1 рад/c.

Неравномерность АЧХ в полосе пропускания, уровень подавления в полосе заграждения и порядок фильтра остаются неизменными.

Аргументы
[in]bУказатель на вектор коэффициентов числителя передаточной функции $H(s)$ исходного аналогового ФНЧ.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[in]aУказатель на вектор коэффициентов знаменателя передаточной функции $H(s)$ исходного аналогового ФНЧ.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[in]ordПорядок исходного фильтра и фильтра после переобразования.

[in]w0Частота среза исходного ФНЧ.

[in]w1Требуемая частота среза ФНЧ после преобразования.

[in,out]betaУказатель на вектор коэффициентов числителя передаточной функции $F(s)$ ФНЧ после преобразования.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

[in,out]alphaУказатель на вектор коэффициентов знаменателя передаточной функции $F(s)$ аналогового ФНЧ после преобразования.
Размер вектора [ord+1 x 1].
Память должна быть выделена.

Возвращает
RES_OK Преоборазование расчитано успешно.

В противном случае код ошибки.
Автор
Бахурин Сергей www.dsplib.org

См. определение в файле filter_ft.c строка 179

◆ ratcompos()

int ratcompos ( double *  b,
double *  a,
int  n,
double *  c,
double *  d,
int  p,
double *  beta,
double *  alpha 
)

Рациональная композиця


Функция рассчитывает композицию вида $Y(s) = (H \circ F)(s) = H(F(s))$, где

\[ H(s) = \frac{\sum\limits_{m = 0}^{n} b_m s^m}{\sum\limits_{k = 0}^{n} a_k s^k}, \quad F(s) = \frac{\sum\limits_{m = 0}^{p} d_m s^m}{\sum\limits_{k = 0}^{p} c_k s^k}, \quad Y(s) = \frac{\sum\limits_{m = 0}^{n p} \beta_m s^m}{\sum\limits_{k = 0}^{n p} \alpha_k s^k} \]

Функция рациональной композиции необходима для произведения частотных преобразований передаточных характеристик аналоговых и цифровых фильтров, а также для билинейного преобразования передаточных характеристик аналоговых фильтров в соответствующие передаточные характеристики цифровых фильтров.

Аргументы
[in]bУказатель на вектор коэффициентов числителя функции $H(s)$.
Размер вектора [n+1 x 1].
Память должна быть выделена.

[in]aУказатель на вектор коэффициентов знаменателя функции $H(s)$.
Размер вектора [n+1 x 1].
Память должна быть выделена.

[in]nПорядок полиномов рациональной функции $H(s)$.

[in]cУказатель на вектор коэффициентов числителя функции $F(s)$.
Размер вектора [p+1 x 1].
Память должна быть выделена.

[in]dУказатель на вектор коэффициентов знаменателя функции $F(s)$.
Размер вектора [p+1 x 1].
Память должна быть выделена.

[in]pПорядок полиномов рациональной функции $F(s)$.

[in,out]betaУказатель на вектор коэффициентов числителя функции $Y(s) = (H \circ F)(s)$.
Размер вектора [n*p+1 x 1].
Память должна быть выделена.

[in,out]alphaУказатель на вектор коэффициентов знаменателя функции $Y(s) = (H \circ F)(s)$.
Размер вектора [n*p+1 x 1].
Память должна быть выделена.

Возвращает
RES_OK Рациональная композиция рассчитана успешно.

В противном случае код ошибки.
Автор
Бахурин Сергей www.dsplib.org

См. определение в файле filter_ft.c строка 204

Используется в bilinear(), low2high() и low2low().