This question is possibly relevant to this one, but it's not duplicate as far as I'm aware.

I'm trying to make phaser with feedback control as shown in the picture below: (Got it from here)

First problem for me was that feedback loop here doesn't have any delay, so at time $n$ I should pass allpass cascade output $y_n$ to its input at this exact time which seemed impossible to me. However, thanks to this comment I was able to figure out how to transform IIR filter cascade to incorporate this feedback into it in terms of its coefficients:

Recursive equation for $k$-th order IIR filter: $$y_n = b_0x_n + b_1x_{n-1} + ... + b_kx_{n-k} - a_1y_{n-1} - a_2y_{n-2} - ... - a_ky_{n-k}$$ As we pass system output through explicit feedback loop our input sequence $x_n$ is defined as follows: $$x_n = d_n + fy_n$$ where $d_n$ is system "dry" input and $f$ is feedback coefficient such as $|f| < 1$

If we plug this into the first equation we'll get $$y_n = b_0(d_n + fy_n) + b_1(d_{n-1} + fy_{n-1}) + ... + b_k(d_{n-k} + fy_{n-k}) - a_1y_{n-1} - a_2y_{n-2} - ... - a_ky_{n-k}=\\ =b_0d_n + b_0fy_n+...+b_kd_{n-k}+b_kfy_{n-k}-a_1y_{n-1}-...-a_ky_{n-k}=\\ = b_0fy_n+b_0d_n+...+b_kd_{n-k}-(a_1-b_1f)y_{n-1}-...-(a_k-b_kf)y_{n-k}$$

Moving all $y_n$ terms to the left and grouping them we'll finally get our desired system coefficients: $$(1-b_0f)y_n = b_0d_n+b_1d_{n-1}+...+b_kd_{n-k}-(a_1-b_1f)y_{n-1}-(a_2-b_2)y_{n-2} -...-(a_k-b_k)y_{n-k}$$

Now, the problem is: we all know we cannot directly use filters with order $k > 2$ without stumbling upon various stability problems. Say, I want to implement phaser with four notches, meaning I'll have allpass cascade of overall order $k = 8$ In that case I would probably want to find poles/zeros of my system analytically (instead of coefficients) and split it to second order sections afterwards or something like that, but I cannot figure out how to find poles/zeros instead of coefficients or is it even possible to do so?

All I've got is that using the fact that for allpass of order $k$ and having poles $[p_1, p_2, ..., p_k]$ its zeros will be always equal to $[1/p_1, 1/p_2, ..., 1/p_k]$ and knowing that all poles of initial system (without feedback) are equal I can write down my problem as follows:

$(x - p)^8 + a(x-1/p)^8= 0$

$p,a \in ℝ$

$|p| < 1$

$|a| < 1$

Also, below you can find matlab code which calculates filter coefficients for the whole phaser, you can see unstable case on [bf, af] = notch_reso2(100, 48000, 0.5, 0.5);

function [ b, a ] = notch_reso2(  fc, fs, feedback, depth  )
order = 8;
w = pi * fc / fs;
a1 = (1 - cot(w)) / (1 + cot(w));
k = a1 ^ order;
p = -a1;
z = 1/p;

zs = repmat(z, 1, order);
ps = repmat(p, 1, order);

b = poly(zs) .* k;
a = poly(ps);

f = feedback;

af = a - b .* f;
bf = b ./ af(1);
af = af ./ af(1);

bf = bf .* depth + af * (1 - depth);
b = bf;
a = af;
end

This question is possibly relevant to this one, but it's not duplicate as far as I'm aware.

I'm trying to make phaser with feedback control as shown in the picture below: (Got it from here)

First problem for me was that feedback loop here doesn't have any delay, so at time $n$ I should pass allpass cascade output $y_n$ to its input at this exact time which seemed impossible to me. However, thanks to this comment I was able to figure out how to transform IIR filter cascade to incorporate this feedback into it in terms of its coefficients:

Recursive equation for $k$-th order IIR filter: $$y_n = b_0x_n + b_1x_{n-1} + ... + b_kx_{n-k} - a_1y_{n-1} - a_2y_{n-2} - ... - a_ky_{n-k}$$ As we pass system output through explicit feedback loop our input sequence $x_n$ is defined as follows: $$x_n = d_n + fy_n$$ where $d_n$ is system "dry" input and $f$ is feedback coefficient such as $|f| < 1$

If we plug this into the first equation we'll get $$y_n = b_0(d_n + fy_n) + b_1(d_{n-1} + fy_{n-1}) + ... + b_k(d_{n-k} + fy_{n-k}) - a_1y_{n-1} - a_2y_{n-2} - ... - a_ky_{n-k}=\\ =b_0d_n + b_0fy_n+...+b_kd_{n-k}+b_kfy_{n-k}-a_1y_{n-1}-...-a_ky_{n-k}=\\ = b_0fy_n+b_0d_n+...+b_kd_{n-k}-(a_1-b_1f)y_{n-1}-...-(a_k-b_kf)y_{n-k}$$

Moving all $y_n$ terms to the left and grouping them we'll finally get our desired system coefficients: $$(1-b_0f)y_n = b_0d_n+b_1d_{n-1}+...+b_kd_{n-k}-(a_1-b_1f)y_{n-1}-(a_2-b_2)y_{n-2} -...-(a_k-b_k)y_{n-k}$$

Now, the problem is: we all know we cannot directly use filters with order $k > 2$ without stumbling upon various stability problems. Say, I want to implement phaser with four notches, meaning I'll have allpass cascade of overall order $k = 8$ In that case I would probably want to find poles/zeros of my system analytically (instead of coefficients) and split it to second order sections afterwards or something like that, but I cannot figure out how to find poles/zeros instead of coefficients or is it even possible to do so?

All I've got is that using the fact that for allpass of order $k$ and having poles $[p_1, p_2, ..., p_k]$ its zeros will be always equal to $[1/p_1, 1/p_2, ..., 1/p_k]$ and knowing that all poles of initial system (without feedback) are equal I can write down my problem as follows:

$(x - p)^8 + a(x-1/p)^8= 0$

$p,a \in ℝ$

$|p| < 1$

$|a| < 1$

Also, below you can find matlab code which calculates filter coefficients for the whole phaser, you can see unstable case on [bf, af] = notch_reso2(100, 48000, 0.5, 0.5);

function [ b, a ] = notch_reso2(  fc, fs, feedback, depth  )
order = 8;
w = pi * fc / fs;
a1 = (1 - cot(w)) / (1 + cot(w));
k = a1 ^ order;
p = -a1;
z = 1/p;

zs = repmat(z, 1, order);
ps = repmat(p, 1, order);

b = poly(zs) .* k;
a = poly(ps);

f = feedback;

af = a - b .* f;
bf = b ./ af(1);
af = af ./ af(1);

bf = bf .* depth + af * (1 - depth);
b = bf;
a = af;
end