2
$\begingroup$

I have an FIR filter that is a low pass at 11025Hz for 44100Hz, 461 taps. The reference implementation, i.e. naive convolution, works as expected. Now I've written a vectorized implementation to take advantage of SIMD/AVX acceleration but there is something wrong with it.

From a hearing standpoint, both versions are almost indistinguishable from one another. I systematically check a filter's spectrum to see if I did things right.

And for that vectorized convolution, I found the following difference:

This is a test against a white noise signal.

The green spectrum is the FFT of the reference filter that works as expected.

The blue spectrum is the FFT of vectorized convolution (wrong) mixed with its phase inverted with the regular convolution (good) to see the difference.

Basically, it appears that there is a high pass filtering occurring as well.

enter image description here

Now when I compare spectrums with a song, they are 99.9% close as you can see:

enter image description here

But when looking at peaks, the differences are drastic:

This is the unfiltered signal:

enter image description here

This is the difference between both filters:

enter image description here

In that difference, you can hear the original signal, filtered with what appears to be a high pass filter.

Ultimately, this isn't a big deal as both filters sound the same, but I would like to understand why is this happening.

This is the vectorized convolution I wrote:

(although it's C#, it gets SIMD-accelerated by Unity Burst)

public void VectorizedOuterInner(float[] source, float[] target, int length)
{
    VectorizedCheck(source, target, length);

    var h = H;
    var z = Z;

    for (var s = 0; s < length; s += 4)
    {
        var dt = Tables[Offset];

        var d0 = dt[0];
        var d1 = dt[1];
        var d2 = dt[2];
        var d3 = dt[3];

        var s0 = s + 0;
        var s1 = s + 1;
        var s2 = s + 2;
        var s3 = s + 3;

        z[d0] = z[d0 + Length] = source[s0];
        z[d1] = z[d1 + Length] = source[s1];
        z[d2] = z[d2 + Length] = source[s2];
        z[d3] = z[d3 + Length] = source[s3];

        ref var y0 = ref target[s0];
        ref var y1 = ref target[s1];
        ref var y2 = ref target[s2];
        ref var y3 = ref target[s3];

        y0 = y1 = y2 = y3 = 0.0f;

        var t = 0;

        for (; t < Length - 4; t += 4)
        {
            var t0 = t + 0;
            var t1 = t + 1;
            var t2 = t + 2;
            var t3 = t + 3;
            var t4 = t + 4;
            var t5 = t + 5;
            var t6 = t + 6;

            var h0 = h[t0];
            var h1 = h[t1];
            var h2 = h[t2];
            var h3 = h[t3];

            var k0 = dt[t0];
            var k1 = dt[t1];
            var k2 = dt[t2];
            var k3 = dt[t3];
            var k4 = dt[t4];
            var k5 = dt[t5];
            var k6 = dt[t6];

            var z0 = z[k0];
            var z1 = z[k1];
            var z2 = z[k2];
            var z3 = z[k3];
            var z4 = z[k4];
            var z5 = z[k5];
            var z6 = z[k6];

            y0 += h0 * z0 + h1 * z1 + h2 * z2 + h3 * z3;
            y1 += h0 * z1 + h1 * z2 + h2 * z3 + h3 * z4;
            y2 += h0 * z2 + h1 * z3 + h2 * z4 + h3 * z5;
            y3 += h0 * z3 + h1 * z4 + h2 * z5 + h3 * z6;
        }

        for (; t < Length; t += 1)
        {
            var h0 = h[t];
            var z0 = z[t + d0];
            var z1 = z[t + d1];
            var z2 = z[t + d2];
            var z3 = z[t + d3];

            y0 += h0 * z0;
            y1 += h0 * z1;
            y2 += h0 * z2;
            y3 += h0 * z3;
        }

        UpdateOffset(4);
    }
}

I pre-compute tables of indices, else it runs slower than naive convolution:

public static int[][] GetTables(int coefficients, int vectorization)
{
    var tables = new int[coefficients][];

    var length = coefficients + (vectorization * 2 - 2);

    for (var i = 0; i < coefficients; i++)
    {
        tables[i] = new int[length];
    }

    for (var i = 0; i < coefficients; i++) // offset
    {
        for (var j = 0; j < length; j++) // tap
        {
            var k = i - j; // delay

            while (k < 0)
            {
                k += coefficients;
            }

            tables[i][j] = k;
        }
    }

    return tables;
}

This is for updating the delay line index:

private void CheckOffset()
{
    Assert.IsTrue(Offset >= 0, Offset.ToString());
}

private void UpdateOffset(int count)
{
    Offset -= count;

    if (Offset < 0)
    {
        Offset += Length;
    }

    CheckOffset();
}

I've tried the following things but with zero improvement:

  • use 64-bit precision instead of 32-bit precision, no changes
  • check/fiddle with the delay/indices, output is immediately wrong, i.e. looks just right

If someone can shed some light on what's going on, that'd be much appreciated!

$\endgroup$

1 Answer 1

2
$\begingroup$

I did not read your code.

I would try to play with delay (like you did) and try feeding both filters a simple digital dirac. Then you can compare the impulse response directly. Perhaps different behaviour is more evident?

Also, try simplifying the convolution kernel to bare minimum (4 taps?) and execute the code.

Notr that depending on compilet flags, the compiler could choose to re-arrange floating point ops so as to change numerical output. That should give only minute differences for something like a FIR filter though.

$\endgroup$
2
  • $\begingroup$ I've followed your first suggestions of trying an impulse and a smaller filter. Been able to understand what does happen, delays are 'escaping' the next iterations and thus, result is zero. However, I haven't yet been able to figure out how to address that. 🤣 $\endgroup$ Jul 24, 2023 at 12:32
  • $\begingroup$ Basically, I was looking at the opposite direction when picking samples, it works now!🥳 $\endgroup$ Jul 27, 2023 at 10:59

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.