1
$\begingroup$

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]: freqz plot

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

$\endgroup$

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Thanks for the plug, Hilmar. I thought maybe one of W3C versions is a better place to go than the original .txt file. I hope my edit is okay with you. $\endgroup$ Commented May 16, 2023 at 18:46
  • $\begingroup$ Certainly ! I bookmarked this for future reference :-) $\endgroup$
    – Hilmar
    Commented May 16, 2023 at 19:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.