Групповая задержка фильтра и ее расчёт

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"


/*filter order */
#define ORD 8

/* Frequency response vector size */
#define T  201
#define N  2000


/*******************************************************************************
 * Нахожу значение y(x0) в табличной функции y(x)
 ******************************************************************************/
double find_value(double* x, double *y, int n, double x0)
{
    double val   = y[0];
    double dxmin = fabs(x[0] - x0);
    int i;
    double dx;
    for(i = 0; i < n; i++)
    {
        dx = fabs(x[i] - x0);
        if(dx < dxmin)
        {
            dxmin = dx;
            val = y[i];
        }
    }
    return val;
}


/*******************************************************************************
 * Main program
 ******************************************************************************/
int main(int argc, char* argv[])
{
    void* hdspl;              /* DSPL handle         */
    void* hplot[4];           /* GNUPLOT handle      */
    hdspl = dspl_load();      /* Load DSPL functions */

    double b[ORD+1];  /* Коэффициенты числителя H(z) эллиптического фильтра */
    double a[ORD+1];  /* Коэффициенты знаменателя H(z) эллиптического фильтра */
    double w[N];      /* нормированная часта */
    double tau_p[N];  /* Фазовая задержка фильтра */
    double tau_g[N];  /* Групповая задержка фильтра */
    double mag[N];    /* АЧХ фильтра */
    double phi[N];    /* ФЧХ фильтра */
    double t[T];      /* время */
    double x[T];      /* входной сигнал */
    double y[T];      /* выходной сигнал */
    double zx[T];     /* треугольная огибающая входного сигнала */
    double tzy[T];    /* время смещенное на величину групповой задержки */
    double zy[T];     /* треугольная огибающая на выходе фильтра */
    
    double g, tg, tp;
    
    int i;
    
    /* Рассчитываем эллиптический БИХ-фильтр */
    iir(0.5, 60.0, ORD, 0.15, 0.0, DSPL_FILTER_LPF | DSPL_FILTER_ELLIP, b, a);
    
    /* вывод коэффициентов фильтра на печать */
    for(i = 0; i <= ORD; i++)
       printf("%.9f    %.9f\n", b[i], a[i]);
    
    /* Заполняю массив нормированной частоты от 0 до pi/4 */
    linspace(0, M_PI*0.25, N, DSPL_PERIODIC, w);
    
    /* рассчет фазовой задержки */
    phase_delay(b,a,ORD, DSPL_FLAG_DIGITAL, w, N, tau_p); 
    
    /* рассчет групповой задержки */
    group_delay(b,a,ORD, DSPL_FLAG_DIGITAL, w, N, tau_g); 
    
    /* рассчет АЧХ и ФЧХ задержки */
    filter_freq_resp(b, a, ORD, w, N, DSPL_FLAG_DIGITAL, mag, phi, NULL);

    /* 
    находим значения АЧХ, фазовой и групповой задержки для частоты 
    входного сигнала pi/8 и выводим значения на печать.
    */
    g  = find_value(w, mag,   N, 0.125*M_PI);
    tg = find_value(w, tau_g, N, 0.125*M_PI);
    tp = find_value(w, tau_p, N, 0.125*M_PI);
    printf(" gain  = %.4f\n tau_g = %.3f\n tau_p = %.3f\n", g, tg, tp);

    /* Сохраняем АЧХ, ФЧХ, фазовую и групповую задержку в файлы для 
       построения графиков */
    writetxt(w, tau_g, N, "dat/tau_g.txt");
    writetxt(w, tau_p, N, "dat/tau_p.txt");
    writetxt(w, mag,   N, "dat/mag.txt");
    writetxt(w, phi,   N, "dat/phi.txt");
    
    /* Формируем сигнал и треугольные огибающие на входе и выходе */
    window(zx, T, DSPL_WIN_BARTLETT | DSPL_SYMMETRIC, 0);
    window(zy, T, DSPL_WIN_BARTLETT | DSPL_SYMMETRIC, 0);
    for(i = 0; i < T; i++)
    {
        t[i] = (double)i;  /* время */
        x[i] = sin(M_PI * 0.125 * t[i]) * zx[i]; /* сигнал на входе */
        tzy[i] = t[i] + tg; /* смещаю огибающую на групповую задержку */
        /* умножаю на найденное значение АЧХ на частоте pi/8 */
        zy[i] = zx[i] * g;  
    }
    /* пропускаю сигнал через фильтр */
    filter_iir(b,a,ORD,x,T,y);
    
    /* сохраняю входной и выходной сигналы, а также треугольные огибающие */
    writetxt(t, x, T, "dat/x.txt");
    writetxt(t, y, T, "dat/y.txt");
    writetxt(t, zx, T, "dat/zx.txt");
    writetxt(tzy, zy, T, "dat/zy.txt");


    /*------------------------------------------------------------------------*/
    /* Построение графиков GNUPLOT */
    gnuplot_create(argc, argv, 800, 600, "img/group_delay.png", hplot);
    gnuplot_cmd(hplot[0], "unset key");
    gnuplot_cmd(hplot[0], "set grid");
    gnuplot_cmd(hplot[0], "unset xlabel");
    gnuplot_cmd(hplot[0], "plot 'dat/tau_g.txt' with lines, 'dat/tau_p.txt' with lines");
    gnuplot_close(hplot[0]);
    
    gnuplot_create(argc, argv, 800, 600, "img/group_delay.png", hplot+1);
    gnuplot_cmd(hplot[1], "unset key");
    gnuplot_cmd(hplot[1], "set grid");
    gnuplot_cmd(hplot[1], "unset xlabel");
    gnuplot_cmd(hplot[1], "plot 'dat/mag.txt' with lines");
    gnuplot_close(hplot[1]);
    
    gnuplot_create(argc, argv, 800, 600, "img/group_delay.png", hplot+2);
    gnuplot_cmd(hplot[2], "unset key");
    gnuplot_cmd(hplot[2], "set grid");
    gnuplot_cmd(hplot[2], "unset xlabel");
    gnuplot_cmd(hplot[2], "plot 'dat/x.txt' with lines, 'dat/zx.txt' with lines");
    gnuplot_close(hplot[2]);
    
    gnuplot_create(argc, argv, 800, 600, "img/group_delay.png", hplot+3);
    gnuplot_cmd(hplot[3], "unset key");
    gnuplot_cmd(hplot[3], "set grid");
    gnuplot_cmd(hplot[3], "unset xlabel");
    gnuplot_cmd(hplot[3], "plot 'dat/y.txt' with lines, 'dat/zy.txt' with lines");
    gnuplot_close(hplot[3]);
    /* free dspl handle */
    dspl_free(hdspl);

    return 0;
}