использовал приведенный данные для расчета. Расчет делал своей библиотекой. Код программы:
Код: Выделить всё
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
/* Порядок фильтра ФНЧ */
#define ORDLP 5
/* Порядок фильтра полосового */
#define ORDBP 10
/* размер векторов частотной характеристики фильтра */
#define N 1000
int main(int argc, char* argv[])
{
void* hdspl; /* DSPL handle */
void* hplot; /* GNUPLOT handle */
/* Коэффициенты числителя ФНЧ H(s) */
double blp[ORDLP+1] = {0.229827, 0.000000, 0.220033,
0.000000, 0.046968, 0.000000};
/* Коэффициенты знаменателя ФНЧ H(s) */
double alp[ORDLP+1] = {0.229827, 0.788084, 1.129145,
1.847068, 0.923344, 1.000000};
/* Коэффициенты числителя полосового H(s) */
double bbp[ORDBP+1] = {0.000000, 14811342.626230, 0.000000,
2023340.076052, 0.000000, 71598.588280,
0.000000, 360.148756, 0.000000, 0.469268, 0.000000};
/* Коэффициенты знаменателя полосового H(s) */
double abp[ORDBP+1] = {2365742167.178879, 291177069.320213,
235456493.051836, 21865990.781777,
7907224.974367, 502677.680863, 105494.693148,
3892.083925, 559.152440, 9.225358, 1.000000};
double bz[ORDBP+1];
double az[ORDBP+1];
double w[N]; /* Циклическая частота рад/c */
double mag_lp[N]; /* АЧХ фильтра ФНЧ */
double mag_bp[N]; /* АЧХ фильтра полосового */
double mag_z[N]; /* АЧХ фильтра цифрового */
int k;
/* Load DSPL function */
hdspl = dspl_load();
if(!hdspl)
{
printf("libdspl loading error!\n");
return -1;
}
/* вектор частоты в логарифмической шакале от 0.01 до 100 рад/c */
logspace(-2.0, 2.0, N , DSPL_SYMMETRIC, w);
/* частотные характеристика фильтра */
filter_freq_resp(blp, alp, ORDLP, w, N,
DSPL_FLAG_LOGMAG| DSPL_FLAG_ANALOG, mag_lp, NULL, NULL);
filter_freq_resp(bbp, abp, ORDBP, w, N,
DSPL_FLAG_LOGMAG| DSPL_FLAG_ANALOG, mag_bp, NULL, NULL);
/* Сохранить характеристики для построения графиков */
writetxt(w, mag_lp, N, "dat/mag_lp.txt");
writetxt(w, mag_bp, N, "dat/mag_bp.txt");
/* билинейное преобразование */
bilinear(bbp, abp, ORDBP, bz, az);
/* печать коэффициентов */
for(k = 0; k < ORDBP+1; k++)
printf("b[%2d] = %- 22.7f a[%2d] = %- 22.7f\n", k, bz[k], k, az[k]);
/* частота цифрового фильтра от 0 до pi */
linspace(0, M_PI, N, DSPL_PERIODIC, w);
filter_freq_resp(bz, az, ORDBP, w, N, DSPL_FLAG_LOGMAG, mag_z, NULL, NULL);
/* Сохранить характеристики для построения графиков */
writetxt(w, mag_z, N, "dat/mag_z.txt");
/* plotting by GNUPLOT */
gnuplot_create(argc, argv, 920, 260, "img/forum_low2bp.png", &hplot);
gnuplot_cmd(hplot, "set logscale x");
gnuplot_cmd(hplot, "unset key");
gnuplot_cmd(hplot, "set grid");
gnuplot_cmd(hplot, "set xlabel 'частота, рад/с'");
gnuplot_cmd(hplot, "set multiplot layout 1,3 rowsfirst");
gnuplot_cmd(hplot, "set ylabel 'АЧХ ФНЧ, дБ'");
gnuplot_cmd(hplot, "set yrange [-100:5]");
gnuplot_cmd(hplot, "plot 'dat/mag_lp.txt' with lines");
gnuplot_cmd(hplot, "set ylabel 'АЧХ полосового, дБ'");
gnuplot_cmd(hplot, "plot 'dat/mag_bp.txt' with lines");
gnuplot_cmd(hplot, "unset logscale x");
gnuplot_cmd(hplot, "set xrange [0:pi]");
gnuplot_cmd(hplot, "set ylabel 'АЧХ цифрового, дБ'");
gnuplot_cmd(hplot, "set xlabel 'нормированная частота 0...pi'");
gnuplot_cmd(hplot, "plot 'dat/mag_z.txt' with lines");
gnuplot_cmd(hplot, "unset multiplot");
gnuplot_close(hplot);
/* free dspl handle */
dspl_free(hdspl);
return 0;
}
В результате выполнения рассчитывает билинейное преобразование, выводятся к-ты цифрового полосового фильтра:
Все вроде нормально. Если надо можете отладить свой код по моей программе. Отладку моей программы можно сделать в среде CodeBlocks.
PS. ФНЧ график совпадает. Для полосового фильтра графики ваш и мой различаются. В частности шкала частот разная на графиках, хотя качественно фильтры похожи.