Group delay of a transfer function

Question

I'm trying to find a straightforward way to find the group delay of a transfer function.

I've created a transfer function, and through the BodePlot function, I'm able to plot its phase characteristic. The group delay is defined as the negative of the derivative of the phase as a function of the frequency ($ g_{d} = - \frac{d\phi \left(\omega\right)}{d\omega} $), indicating that the group delay of this transfer function has a value of approximately 5 in the interval $ 0 \le \omega \le 1$.

I can obtain the phase through the Arg or the ArcTan functions, but as these functions output values in the range $ \left[ -\pi, \pi \right]$ and $ \left[ -\frac{\pi}{2}, \frac{\pi}{2} \right]$, respectively, the phase is wrapped between these limiting values. It seems that the Derivative function does not work with Arg, Re and Im functions (Derivative of real functions including Re and Im).

I think I can plot the group delay by computing the numerical derivative of the transfer function's phase. However, I need an unwrapped phase function, as that given by the phase plot, or I will obtain discontinuities in my group delay plot. Please, take a look at the code below.

ClearAll["Global`*"] (* clear all symbols *)
p[s_] := 0.3576 - 0.1394*s^2 + 0.1721*s^4 - 0.1404*s^6
e[s_] := 0.3620 + 1.6452*s + 3.7021*s^2 + 4.9510*s^3 + 4.5715*s^4 +   2.5235*s^5 + s^6
s21[s_] := p[s]/e[s] (* transfer function *)
tfm = TransferFunctionModel[s21[s], s]
BodePlot[tfm, {0.03, 1.5}, PlotLayout -> "Phase", ImageSize -> Large, GridLines -> Automatic, ScalingFunctions -> {"Linear", "Radian"}] (* phase plot of the transfer function *)
Plot[Arg[s21[I*bigOmega]], {bigOmega, -1.5, 1.5}, GridLines -> Automatic, PlotRange -> All, Exclusions -> None] (* Trying to plot the phase with the Arg function *)
Plot[ArcTan[Im[s21[I*bigOmega]]/Re[s21[I*bigOmega]]], {bigOmega, -1.5, 1.5}, GridLines -> Automatic, PlotRange -> All, Exclusions -> None] (* Trying to plot the phase with the ArcTan function *)
Plot[Derivative[Arg[s21[I*bigOmega]]], {bigOmega, 0.03, 0.5}, GridLines -> Automatic, PlotRange -> All, Exclusions -> None] (* do not work *)
Plot[Derivative[ArcTan[Im[s21[I*bigOmega]]/Re[s21[I*bigOmega]]]], {bigOmega, 0.03, 0.5}, GridLines -> Automatic, PlotRange -> All, Exclusions -> None] (* do not work *)
Needs["NumericalCalculus`"];
Plot[-ND[Arg[s21[I*bigOmega]], bigOmega, omega], {omega, -1.5, 1.5}, GridLines -> Automatic, PlotRange -> 10, Exclusions -> None] (* Find the group delay by computing the numerical derivative of the phase *)

In Implementing continuous phase/Arg function, the code below is suggested to plot the unwrapped phase. However, my knowledge of Mathematica is insufficient to understand the following snippet and adapt it to my needs.

SetAttributes[argPlot, HoldAll];
Options[argPlot] = Options[Plot];

argPlot[exp_, {x_, x0_, x1_}, opt : OptionsPattern[argPlot]] := 
 Module[{pts, pl}, pl = Plot[Arg[exp], {x, x0, x1}, PlotRange -> All, PlotPoints -> OptionValue[PlotPoints]];
  pts = SortBy[Cases[pl, Line[pts_] :> pts, Infinity], #[[1, 1]] &];
  pts = Reap[Fold[Module[{ptsn}, ptsn = #2;
          ptsn[[All, 2]] -= Round[ptsn[[1, 2]] - #1, 2 Pi];
          Sow[ptsn];
          ptsn[[-1, 2]]] &, 0, pts];][[2, 1]];
  ListLinePlot[Flatten[pts, 1], opt]]

argPlot[s21[I*omega], {omega, 0.03, 1.5}, GridLines -> Automatic]

So, does anyone knows how to obtain an unwrapped phase function so as it is possible to compute the numerical derivative and obtain the proper group delay? or, is there an easier way to obtain the group delay?

Answer 1

To wrap the phase range, use the option PhaseRange.

opts = Sequence[PlotLayout -> "Phase", PhaseRange -> {-\[Pi], \[Pi]}, 
ImageSize -> Large, GridLines -> Automatic, ScalingFunctions -> {"Linear", "Radian"}, PlotRange -> All];

BodePlot[tfm, {0, 2}, Evaluate@opts]

To compute the group delay, first compute the sinusoidal transfer function.

stf = tfm[I w][[1, 1]];

Assume that its real and imaginary parts are $x[w]$ and $y[w]$. Then the group delay can be computed as

gd = -D[ArcTan[x[w], y[w]], w]

Now compute the above values for the sinusoidal transfer function and substitute it to get group delay as a function of $w$.

values = Thread[{x[w], y[w]} -> (Together@ComplexExpand@#@stf & /@ {Re, Im})];
values = Join[values, D[values, w]];

Compute its maximum.

FindMaximum[gd /. values, w]

{5.57155, {w -> 1.05231}}

Plot it.

Plot[Evaluate[gd /. values], {w, 0, 2}, GridLines -> Automatic]

Answer 2

I have managed to find a relative simple solution to plot the group delay of a transfer functrion through the data taken from the BodePlot.

The following snippet plot the desired group delay.

pointsList = BodePlot[tfm, {0, 1.5}, ImageSize -> Large, GridLines -> Automatic, 
   ScalingFunctions -> {{"Linear", Automatic}, {"Linear", 
   "Radian"}}][[1, 2, 1, 1, 1, 1, 3, 2, 1]];
b = Interpolation[pointsList, Method -> "Spline"];
c = b';
Plot[Evaluate[-c[x]], {x, 0, 1.5}, GridLines -> Automatic]

Another option is

pointsList = BodePlot[tfm, {0, 1.5}, ImageSize -> Large, GridLines -> Automatic, 
   ScalingFunctions -> {{"Linear", Automatic}, {"Linear", "Radian"}},
   PlotPoints -> 25][[1, 2, 1, 1, 1, 1, 3, 2, 1]];
derA = Differences[pointsList] /. {x_, y_} -> y/x;
lst = Transpose[{Drop[Transpose[pointsList][[1]], 1], -derA}];
ListLinePlot[lst, GridLines -> Automatic, PlotRange -> {2, 6}]

In this second snippet, the ListLinePlot resolution can be improved by increasing the PlotPoints in the BodePlot function. The first snippet generally presents a good resolution, as the phase function is interpolated.

I'm using Mathematica version 11.0.0.0 for Windows. Maybe, in another platform or software version, one should adjust the object position (the [[1, 2, 1, 1, 1, 1, 3, 2, 1]]) to recover just the phase plot point coordinates, as the code does. Otherwise, one could improve de code through the use of pattern matching and the Module function.

Obviously, the BodePlot also provides me the unwrapped phase of the transfer function.