This question is possibly relevant to [this](https://dsp.stackexchange.com/questions/15763/resonant-peak-frequencies-for-phaser-effect) 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](https://en.wikipedia.org/wiki/Phaser_%28effect%29))

[![reference phaser design][1]][1]


  [1]: https://i.sstatic.net/0dI8y.png

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](https://dsp.stackexchange.com/questions/15763/resonant-peak-frequencies-for-phaser-effect#comment27282_15767) 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