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

Question

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?

Answer 1

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$.

Answer 2

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

Youtube tutorial

Tutorial's Point tutorial