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.
Now when I compare spectrums with a song, they are 99.9% close as you can see:
But when looking at peaks, the differences are drastic:
This is the unfiltered signal:
This is the difference between both filters:
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!
