Allpass Filter Gain Issue

Question

Background

I am having issues implementing an allpass filter to model wave dispersion in a stiff string. In order to simulate wave propagation in a string, I am using a digital waveguide. I implemented the matlab code referenced in this paper in c++, and my code produces the same coefficients. This leads me to believe that my error is in the implementation of the allpass filter. I used this as a reference when implementing the allpass biquads. Relevant c++ code is shown below:

//args are: x[n], x[n - 1], x[n - 2], y[n - 1], y[n - 2]
    float dispersion_filter(float x, float x1, float x2, float y1, float y2){
            double y = 0.0f;
            for (int i = 0; i < disp_filter_coeffs.size(); i++) {
                float a1 = disp_filter_coeffs[i][0];
                float a2 = disp_filter_coeffs[i][1];
                float b1 = a1/a2;
                float b2 = 1.0/a2;
                float g = a2;
                y = 0.45 * (x + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2);
                if (i < disp_filter_coeffs.size() - 1){
                x2 = x1;
                x1 = x;
                y2 = y1;
                y1 = y;
                x = y;
                    }
            }
            return y;
    }

And the function that generates the coefficients:

    void calculate_disp_filter_coeffs(double f0, double B, double bw, double beta){
        double tau0 = rate / f0;
        double pd0 = tau0 / std::sqrt(1.0 + B * std::pow(bw/f0, 2));
        double mu0 = pd0 / (1.0 + B * std::pow(bw/f0, 2));
        double phi0 = 2.0 * M_PI * (bw / rate) * pd0 - mu0 * 2.0 * M_PI * bw/rate;
        int n_ap = std::floor(phi0/(2.0 * M_PI)); //number of filters
        disp_filter_coeffs.resize(n_ap);

        double last_theta;
        for (int i = 0; i < n_ap; i++){
            double phi0 = M_PI * ((i * 2) + 1);
            double eta0 = i / (1.2 * n_ap) * bw;
            double pd = tau0 / std::sqrt(1 + B * std::pow(eta0/f0, 2));
            double tau = pd / (1 + B * std::pow(eta0/f0, 2));
            double phi1 = 2.0 * M_PI * (eta0 / rate) * pd;
            double theta = rate / (2.0 * M_PI) * (phi0 - phi1 + (2.0 * M_PI / rate) * eta0 * tau) / (tau - mu0);
            double delta = (i == 0) ? theta : (theta - last_theta) / 2.0;
            last_theta = theta;
            double cc = std::cos(theta * (2.0 * M_PI / rate));
            double eta1 = (1.0 - beta * std::cos(delta * 2.0 * M_PI / rate)) / (1.0 - beta);
            double alpha = std::sqrt(std::pow(eta1, 2) - 1.0) - eta1;

            double a1 = 2 * alpha * cc;
            double a2 = std::pow(alpha, 2);
            disp_filter_coeffs[i][0] = a1;
            disp_filter_coeffs[i][1] = a2;
            }
    }

And how the filter is run:

float tmp1 = dispersion_filter(r, disp_x1, disp_x2, disp_y1, disp_y2);
    disp_x2 = disp_x1;
    disp_x1 = r;
    disp_y2 = disp_y1;
    disp_y1 = tmp1;

My Issue

There are a couple issues with this implementation. Firstly, the filter is not at 0db gain at all frequencies. Below is the plot of freqz(b, a, 2048) of one allpass biquad section with the coefficients of a = [1, -1.74296, 0.76196], and b = [1, -2.28746, 1.312399]:

This shows that the filter has a gain of about +2.4db. This issue causes the system to be very unstable and generate increasing large values.

Secondly, besides the issue of instability, there is also the issue of the delay being incorrect. As shown in fig. 1, the phase delay is not correct. The delay should be greater for lower frequencies, and smaller for higher frequencies. The freqz plot shows the opposite of this.

The full(and unfinished) StiffString class, for context: Pastebin

Answer 1

That all seems rather tortured and needlessly complicated.

If you want to design an allpass filter (2nd order) at a specific resonant frequency and Q, you can simply use the Audio EQ Cookbook.

Here is a simple Matlab example

%% allpass design
fs = 48000; % sample rate
f0 = 1000; % corner frequency
Q = 2;  % desired Q

%% design
w0 = 2*pi*f0/fs; % normalized frequency
a0 = [1+sin(w0)/2/Q -2*cos(w0) 1-sin(w0)/2/Q]; % denominator
a = a0/a0(1); % normalzed denominator
b = [a(3) a(2) a(1)];  % simply flip the coefficients to get the numerator

%% plot and pretty up the graph a bit
freqz(b,a,1024,fs);
h = get(gcf,'Children');
set(h,'XScale','log');
set(h(2),'YTick',[-360:90:0]);
set(h(1),'YLim',[-1 1]);

This shows that the filter has a gain of about +2.4db

It sure does. Your filter scaling is wrong. Your numerator coefficients should be the reverse of the denominator coefficients.

As shown in fig. 1, the phase delay is not correct.

Fig. 1 doesn't show the phase delay, it shows the phase, which is indeed correct. The phase delay would be the phase divided by the negative frequency, so the shape of the graph would flip. However, for most applications, the relevant quantity if interest is the group delay (not the phase delay). That's defined as the negative derivative of the phase vs frequency and that is indeed be higher at low frequencies.