#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
/*filter order */
#define ORD 8
/* signal length */
#define T 201
/* Frequency response vector size */
#define N 2000
/*******************************************************************************
* Find y(x0) in table 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) coefficients */
double a[ORD+1]; /* H(z) coefficients */
double w[N]; /* normalized frequency */
double tau_p[N]; /* phase delay */
double tau_g[N]; /* group delay */
double mag[N]; /* magnitude */
double phi[N]; /* phase response */
double t[T]; /* time */
double x[T]; /* input signal */
double y[T]; /* output signal */
double zx[T]; /* triangle envelope in the input */
double tzy[T]; /* Envelope time shif vector */
double zy[T]; /* triangle envelope in the output */
double g, tg, tp;
int i;
/* Elliptic IIR filter order 8 */
iir(0.5, 60.0, ORD, 0.15, 0.0, DSPL_FILTER_LPF | DSPL_FILTER_ELLIP, b, a);
/* print filter coefficients */
for(i = 0; i <= ORD; i++)
printf("%.9f %.9f\n", b[i], a[i]);
/* fill frequency array from 0 to pi/4 */
linspace(0, M_PI*0.25, N, DSPL_PERIODIC, w);
/* phase delay calculation */
phase_delay(b,a,ORD, DSPL_FLAG_DIGITAL, w, N, tau_p);
/* group delay calculation */
group_delay(b,a,ORD, DSPL_FLAG_DIGITAL, w, N, tau_g);
/* magnitude and frequency response */
filter_freq_resp(b, a, ORD, w, N, DSPL_FLAG_DIGITAL, mag, phi, NULL);
/*
find magnitude, phase and group delay corresponds to pi/8 frequency .
*/
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);
/* save filter parameters for plotting */
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");
/* input signal generation */
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; /* time */
x[i] = sin(M_PI * 0.125 * t[i]) * zx[i]; /* signal */
tzy[i] = t[i] + tg; /* time shift to tthe group delay */
zy[i] = zx[i] * g; /* scale output envelope */
}
/* filter output */
filter_iir(b,a,ORD,x,T,y);
/* save signals and envelopes */
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");
/*------------------------------------------------------------------------*/
/* plottyng by 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;
}
Group delay and it calculation
Content
Found a mistake?
Select it with the mouse
and press
Introduction
Analog and digital filters are characterized by
transfer functions
and 
that can be written as rational functions of complex variables 
and 
:
(1)
it is assumed that 
.
From filter transfer function we can get a complex frequency response. The analog filter complex frequency response 
we can get if substitute 
, and for digital filter 
.
Frequency responses and 
for analog and digital filters are:
(2)
Complex frequency response can be described in terms of module and phase, then we obtain magnitude and phase response characteristics.
Magnitude and frequency response of analog filter:
(3)
Magnitude and frequency response of analog filter:
(4)
Frequency response characterizes the selective properties of the filter, and phase response — phase shifts that will acquire different signal frequencies when passing through this filter.
Relationship of phase response and time delay of a harmonic signal. Filter phase delay
For a harmonic sinusoidal signal, the phase shift corresponds to the time shift, as shown in the figure 1.
Figure 1. Phase shift and time delay
Delayed signal can be written as:
(5)
here delay 
corresponds to phase shift 
on the 
frequency:
(6)
Then, when the signal passes through a filter with a phase response 
, the time delay of the harmonic signal at the output will depend on the frequency:
(7)
Parameter 
is called the phase delay of the filter and determines the harmonic signal delay.
Relationship of phase response and signal envelope delay. Filter group delay
If the input signal is not just a harmonic oscillation with frequency , but has some envelope in the band from to 
, then different frequencies from a given band will have different delays.
We can formally consider the increments:
(8)
(9)
Subtracting the first from the second we get:
(10)
(11)
In the passage to the limit at 
, we obtain the derivative of the phase response with respect to frequency, and the delay 
characterizes the group delay:
(12)
The meaning of the group delay 
is that the signal envelope is delayed by this amount. The harmonic carrier is delayed by the amount of the phase delay 
.
Let's demonstrate this with an example.
Consider the frequency response 
, phase response 
, phase delay 
and group time 8th order digital elliptic lowpass filter with coefficients:
(13)
Figure 2 shows the graphs of the frequency response 
, phase response , phase delay and group delay 
of this filter.
Figure 2. Frequency response, phase response, phase delay and group delay of the 8th order digital elliptic lowpass filter
Let also we have 
samples of the harmonic signal at the normalized frequency 
with triangular amplitude envelope 
:
(14)
as it showed on the figure 3a.
From the figure 2 we can see that filter phase delay at the frequency equals 
samples, but the group delay 
almost two times more.
Let us pass the signal 
through this elliptical filter and analyze the signal at the output of the filter, shown in Figure
b with a solid line. You can see that the output is delayed relative to the input.
Figure 3. Passing a triangular envelope signal through a specified filter
Consider the delay of a harmonic signal at the frequency 
, as shown in Figure 3 in when zoomed in along the time axis. You can see that the harmonic oscillation is delayed by the amount of the phase delay 
samples. This can be seen from both swing highs and zero crossing points.
Consider now the delay of the triangular envelope at the filter output, as shown with an upscaling in Figure 3 d. You can see that the triangular envelope is delayed by the group delay amount 
samples.
Thus, we have demonstrated that the signal passing through the filter is characterized by two parameters: the phase delay 
, which characterizes the delay of the harmonic signal as a result of the phase shift, and also the value of the group delay time (group delay) , which determines the delay of the envelope of the harmonic signal, and can differ significantly from the phase delay.
When transmitting information, the harmonic oscillation is not informational, therefore, in practice, the value of the group delay time is a more frequently used parameter, since it is it that determines the delay of the information envelope when the signal passes through the filter.
Digital IIR filters group delay calculation
The calculation of the group delay 
based on the numerical differentiation of the phase response according to (12) can be difficult due to the presence of phase response gaps, as shown in Figure 3. Differentiation results in outliers in the group delay that interfere with the analysis.
To improve the quality of the group delay calculation, you can use the algorithm given in this paragraph for the digital IIR filter.
Digital IIR filter Frequency response :
(15)
According to Euler's formula:
(16)
here
(17)
In the future, derivatives will also come in handy:
(18)
Let's continue the transformation of the complex frequency response (16):
(19)
Let us omit the dependence on 
to simplify the notation:
(20)
(21)
Simplified notation of :
(22)
and group delay equation is:
(23)
Calculate the derivative of the arctangent:
(24)
Let's rewrite and get:
(25)
Denote the variables:
(26)
also take account that
(27)
Finally we can write:
(28)
Note that when calculating 
and 
analytical derivatives are used, as a result. (28), will not experience distortions of the numerical differentiation of the phase response.
Analog filters group delay calculation
Complex frequency response of analog filter is:
(29)
Group delay equation (28) can be used for analog filters if we will find in form:
(30)
Rewrite (30) in follow form:
(31)
here
(32)
(33)
The functions 
, 
, 
and 
for analog filters are equals:
(34)
(35)
(36)
(37)
We also need to calculate derivatives according (28):
(38)
(39)
(40)
(41)
Analog filter group delay can be calculated according (28) if we substitute equations (37) for , , and and also its derivatives (41).
Group delay function in DSPL-2.0
The given algorithm for calculating the group delay of analog and digital filters is implemented in software in the DSPL-2.0. library. The phase_delay function implements the calculation of the phase delay of the filter specified by the vectors of the transfer function coefficients, and the group_delay function calculates the group delay of the filter.
The data for the test case showing the difference between the phase and group delay of the filter is also calculated using the DSPL-2.0.
The source code of the group_delay_test program for calculating data for building Figures2 and3 is shown in the following listing:
Conclusions
In this section, we examined such characteristics of digital and analog filters as phase and group delay.
It was shown that the value of the phase delay determines the time delay of the harmonic signal at the filter output, which is determined by the phase shift in accordance with the value of the phase response at a given frequency.
Unlike phase delay, the group delay determines the amount of time delay of the harmonic envelope.
Analytical expressions for calculating the group delay of an analog and digital filter without numerical differentiation of the phase response are given.
Reference
[1]
Approximation Methods for Electronic Filter Design.
New York, McGraw-Hill, 1974, 388 p.
[3]
Discrete-Time Signal Processing (3rd Edition).
Prentice-Hall, 2009.
Page update: 30.08.2020 (13:10:15)