Preface
According to traditional textbook DSP, digital IIR filter structures achieve phase shifting by using cascaded single sample delay blocks. This is in contrast to ‘real’ filters in the analog domain, where phase shifting is due electronical properties of capacitors and inductors  these electronic circuits contain no signal delaying elements and can be then called zero delay.
It then seems like achieving zero delay filtering in the digital domain is impossible, since we are stuck with delay elements by definition. This is not the case. To see how this is possible, we first have to understand what digital IIR filters are mathematically.
Almost all the engineering tools in DSP are tools for simplifying the handling of Ordinary Differential Equations with the use of transfer functions. Handling and solving general ODEs requires a significant amount of mathematical baggage and this is the reason why DSP has not adopted the usage of these equations directly, but have instead developed their own jargon and engineering tradition based on the discretization of ODEs.
‘Infinite Impulse Response’ filter is then nothing more than a solver for a differential equation. This is also true for an analog circuit, which is just a physical implementation of a differential equation solver.
In any case, a practical ODE solver (in contrast to the exact closedform symbolic solution) is based on integration. Traditional DSP techniques involve transforming the ODE to a difference equation which is then integrated numerically.
This is usually done by means which is equivalent to the explicit Euler method, the most basic approximation method available for solving ODEs by explicit integration. The resulting structure is simple enough to be interpreted as a cascade of delay elements  hence all the silly diagrams of delay element networks in DSP textbooks.
Defining zero delay
What is then a ‘zero delay’ digital filter? Nothing more than a numerical integration method. Now, instead of working with explicit discrete delayed values, you use various ways to get a better implicit approximation of the derivative you are going to integrate the state of the equation with.
This might seem like cheating, to call something zero delay when you’re obviously still working with discrete delayed values and I agree to the notion to some degree. Zero delay seems like a confusing and misleading term in this context.
From a differential equation to an algorithm
Let’s define a state variable filter as an ordinary second order differential equation:
where $u(t)$ is the input function.
Knowing the initial conditions of the state functions and the input function, a closedform solution for $x_1$ and $x_2$ could be constructed, but we are interested in transforming this equation to an integration algorithm.
We are now going to look into explicit integration methods.
Euler method
The traditional way to do this is done by applying the Euler method to discretize the differential equation to an explicit difference equation.
In the method, the derivatives in the equation are evaluated as $y’(t)=f(t, y)$ and the algorithm for integrating the state is
where $h=t_{i+1}t_{i}$ is the time step between $t_{i+1}$ and $t_{i}$.
Applying this to the differential equation gets us
It is now easy to form a practical pseudocode integration algorithm out of this:
1 2 3 

To get the well known Chamberlin SVF you need to apply a slight modification and use the semiimplicit Euler method:
and in pseudocode, replacing the state variables with more descriptive names:
1 2 3 

Heun’s method
To improve the Euler method, you can form a corrector step after the first Euler integration step. This is called the Heun’s method:
Being an explicit method, although comparable to implicit trapezoidal integration, it is straightforward to apply to pseudocode:
1 2 3 4 5 6 7 8 9 10 11 

We have looked at the two simplest explicit integration methods. In the next part we will go into the implicit methods and build a ‘zero delay’ digital state variable filter.
Happy hacking!