Deriving Frequency Response for 2-pole Zero-Delay Feedback State Variable Filter

I have an existing zero-delay feedback (ZDF) 2-pole state variable filter implementation (along the lines of the theory presented in VA Filter Design by V. Zavalishin), and I wish to determine the frequency response so that I can plot it. My usual approach with other filters (e.g. biquads) has been to determine the difference equation by inspection, apply the Z-transform, and then substitute $z=e^{j\omega}$ and plot phase and magnitude against $\omega$ from $0$ to $\pi$.

However I am having difficulty determining the difference equation of this filter.

Here is a diagrammatic representation of the implementation I have: Here's the implementation of an iteration, corresponding to a single input sample, in pseudo-code. Fs is the sample rate, F is the filter frequency of interest, R is the resonance, x is the input sample, and lp, bp, hp represent the low, band and high-pass outputs respectively:

# z1 and z2 contain state from previous iteration

k = tan(pi / Fs * F)
b = R + k
a = 1 / (1 + b * k)

hp = (x - (z2 + z1 * b)) * a
x1 = hp * k  # temporary
bp = x1 + z1
x2 = bp * k  # temporary
lp = x2 + z2

z1 = x1 + bp
z2 = x2 + lp

For analysis, I assign $v[n]$ to be the output of z1, and $w[n]$ to be the output of z2. From this I am able to determine by inspection the following equations:

$$hp[n] = x[n] - (b.v[n] + w[n])$$ $$v[n] = k.hp[n-1] + bp[n-1]$$ $$w[n] = k.bp[n-1] + lp[n-1]$$ $$bp[n] = k.hp[n] + v[n]$$ $$lp[n] = k.bp[n] + w[n]$$

I'm not sure how to handle such a system of equations algebraically. For example, if I wish to analyse the low-pass response, then I will need to determine a difference equation for $lp[n]$ in terms of $x[n]$. However I have not had success manipulating these equations by rearrangement and substitution even with the help of the Z transform. Is it possible to do this or are there too many unknown variables?

If an algebraic solution is impossible, is there an alternative method to determine the frequency response of this filter?

Solving for the input-output relations of such a structure may be a bit tedious, but not difficult. For the sake of generality I define the value of the first multiplier $k$ (closer to the input) as $c$, and the value of the second multiplier $k$ as $d$. In your case $c=d=k$. Furthermore I define an auxiliary signal $w[n]$ (with $\mathcal{Z}$-transform $W(z)$) at the output of the adder at the bottom of your figure, which is fed back to the input. The output signals are $y_h[n]$, $y_b[n]$, and $y_l[n]$ for the high-pass, band-pass and low-pass signals, respectively.

In the $\mathcal{Z}$-transform domain you get the following $4$ equations for the $4$ unknowns ($W(z)$, $Y_h(z)$, $Y_b(z)$, $Y_l(z)$):

\begin{align}W(z)&=bz^{-1}[cY_h(z)+Y_b(z)]+z^{-1}[dY_b(z)+Y_l(z)]\\ Y_h(z)&=a[X(z)-W(z)]\\ Y_b(z)&=c(1+z^{-1})Y_h(z)+Y_b(z)z^{-1}\\ Y_l(z)&=d(1+z^{-1})Y_b(z)+Y_l(z)z^{-1}\end{align}

This system of equations can be solved for the (transforms of the) three output signals in terms of the input $X(z)$. The solutions have the form

\begin{align}\frac{Y_b(z)}{X(z)}&=g_b\frac{(1+z^{-1})(1-z^{-1})}{P(z^{-1})}\\\frac{Y_h(z)}{X(z)}&=g_h\frac{(1-z^{-1})^2}{P(z^{-1})}\\\frac{Y_l(z)}{X(z)}&=g_l\frac{(1+z^{-1})^2}{P(z^{-1})}\end{align}

where the polynomial $P(z^{-1})$ has the form

$$P(z^{-1})=1+Az^{-1}+Bz^{-2}$$

and where all constants ($A$, $B$, $g_b$, $g_h$, $g_l$) of course depend on the multiplier values $a$, $b$, $c$, and $d$.

• Thank you. My algebra was leading my into trouble as I ended up with [n-1], [n-2], [n-3], ... terms. I did not think to convert to the Z transform earlier. However, the bit that I would ask for an additional hint with please - how do you "solve" that system of equations? Does the standard linear algebra technique using matrices for solving four equations in four unknowns apply? I just need prompting in the right direction here as it's been a long time since I've tackled this sort of thing. It is reassuring to know that it is possible though - I just wasn't sure if I was wasting my time. May 20 '17 at 23:58
• @meowsqueak: It's just like solving any other system of linear equations. In this case you could: plug Eq1 into Eq2, use Eq3 to express $Y_h$ in terms of $Y_b$ and plug that into Eq2. Then you can express $Y_b$ in terms of $X$. Now you can use Eqs 3 and 4 to express $Y_h$ and $Y_l$ in terms of $Y_b$. May 21 '17 at 8:59
• I think $W(z)$ should be $W(z)=bz^{-1}[cY_h(z)+Y_b(z)] + z^{-1}(dY_b(z)+Y_l(z))$ otherwise only the feedback from the first integrator stage is used. May 22 '17 at 3:40
• @meowsqueak: You're right, I forgot to type the second part of the equation. Corrected. May 22 '17 at 7:00

I do not see any reason why this could not be derived algebraically. But this would be a good application for Mason's gain formula; the references below may be helpful for that approach:

Rick Lyon's Tutorial

Mason's Gain Formula - Wikipedia