I did use Differential Evolution to calculate the coefficients. But you can re-design the filter pair easily using the HIIR library by Laurent de Soras (its source code will automatically unzip to a subdirectory hiir). You can use this C++ HilbertDesign.cpp source and compile with g++ using the compile-command quoted on the first line:
// -*- compile-command: "g++ HilbertDesign.cpp -std=c++11 -msse -I. -g3 -O0 hiir/PolyphaseIir2Designer.cpp -o HilbertDesign" -*-
#include <stdio.h>
#include "hiir/PolyphaseIir2Designer.h"
const int numCoefs = 8; // Number of coefficients, must be even
double transition = 2*20.0/44100; // Sampling frequency is 44.1 kHz. Approx. 90 deg phase difference band is from 20 Hz to 22050 Hz - 20 Hz. The transition bandwidth is twice 20 Hz.
double coefs[numCoefs];
int main() {
hiir::PolyphaseIir2Designer::compute_coefs_spec_order_tbw (coefs, numCoefs, transition);
printf("Phase reference path c coefficients:\n");
for (int i = 1; i < numCoefs; i += 2) {
printf("%.20f,", coefs[i]);
}
printf("\n+90 deg path c coefficients:\n");
for (int i = 0; i < numCoefs; i += 2) {
printf("%.20f,", coefs[i]);
}
printf("\n");
return 0;
}
Running HilbertDesign outputs:
Phase reference path c coefficients:
0.47944111608296202665,0.87624358989504858020,0.97660296916871658368,0.99749940412203375040,
+90 deg path c coefficients:
0.16177741706363166219,0.73306690130335572242,0.94536301966806279840,0.99060051416704042460,
You'd use a coefficient $c$ in a section with transfer function:
$$H(z) = \frac{c - z^{-2}}{1 - cz^{-2}},\tag{1}$$
implementable by:
$$y[k] = c\,\left(x[k] + y[k-2]\right) - x[k-2],\tag{2}$$
with input $x$ and output $y$.
Each of the two paths is a cascade of a number of such sections. The phase reference path should be delayed by an additional one-sample delay. The phase difference between the two paths will be approximately 90 degrees over a band centered about sampling frequency / 4. This symmetry gives the somewhat sparse transfer function (Eq. 1) with a simple computational implementation of the filter (Eq. 2). The coefficients are optimal in equiripple sense.
I used HIIR v1.20. The c coefficients are equivalent to my a^2 coefficients. The square roots of the coefficients from HIIR are quite close to my original pole locations, a, here interleaved for the two paths which gives them in ascending order:
a sqrt(c) a-sqrt(c)
0.4021921162 0.402215635 -2.35187886997168E-005
0.6923878 0.6924168658 -2.9065827921726E-005
0.8561710882 0.8561932617 -2.21734129581819E-005
0.9360654323 0.9360788374 -1.34051398409385E-005
0.9722909546 0.972297804 -6.84943738937793E-006
0.9882295227 0.9882322446 -2.72186419525422E-006
0.9952884791 0.9952891611 -0.000000682
0.9987488453 0.9987489195 -7.4186057630321E-008
There are two ways to obtain a true Hilbert transformer by forward-backward (ping-pong) IIR filtering, here expressed using the frequency responses $H_{\text{ref}}(\omega)$ and $H_{\text{ref}+90^\circ}(\omega)$ of the all-pass branches, with the subscript denoting the approximate phase shift. The time reversal of a conjugate-symmetric frequency response $H(\omega)$ results in a frequency response $H(-\omega)$, which is the signature of a backwards pass on the left sides of:
$$H_{\text{ref}+90^\circ}(\omega)H_{\text{ref}}(-\omega) = H_{\text{ref}+90^\circ}(\omega)H_{-\text{ref}}(\omega) = H_{90^\circ}(\omega),$$
and
$$-H_{\text{ref}}(\omega)H_{\text{ref}+90^\circ}(-\omega) = -H_{\text{ref}}(\omega)H_{-\text{ref}-90^\circ}(\omega) = -H_{-90^\circ}(\omega) = H_{90^\circ}(\omega).$$
Furthermore, Martin Vicanek shows in a 2015—2022 article A New Reverse IIR Filtering Algorithm how to obtain a causal, truncated approximation of each the forward IIR filter impulse response and a causal delayed version of the backward IIR filter impulse response, using low-complexity FIR filters, to obtain an approximate delayed Hilbert transformer.
