Convolution by multiplying list of numbers in memory, so an inverse convolution algorithm?

307 views Asked by At

While "convolutions are multiplications in the frequency domain", they also seem also to be multiplications much more literally. For example if I adjoin numbers adjacently in memory into a sort of list:

5:  (byte) 00000101
22: (byte) 00010110
7:  (byte) 00000111
8:  (byte) 00001000

Becomes:

(Int32) 00000101000101100000011100001000

Then multiply this resulting list/number by another such list, the result is a convolution of one list with the other. The following C# program demonstrates this:

using System;

public class Program
{
    public static void Main(string[] args)
    {
        var A = new byte[]{1, 2, 3, 4, 5, 6, 7, 8};
        var a = combine_to_Int64(A);
        var B = new byte[]{1, 2, 3, 4, 5, 6, 7, 8};
        var b = combine_to_Int64(B);
        var c = a * b;
        var C = separate_to_8_byte(c);
        var d = deconvolution(c, b); // not correct
        var D = separate_to_8_byte(d);
        Console.WriteLine(
            "The convolution of (" + to_string(A) + ") with ("
            + to_string(B) + ") is: (" + to_string(C) + ")");
        Console.WriteLine(
            "The deconvolution of (" + to_string(C) + ") with ("
            + to_string(B) + ") is: (" + to_string(D) + ")");
    }

    public static Int64 convolution(Int64 a, Int64 b)
    {
        return a * b;
    }

    private static Int64 combine_to_Int64(byte[] n)
    {
        if (n.Length == 4)
            return n[0] + 65536 * n[1] + 4294967296 * n[2] +
                281474976710656 * n[3];
        else if (n.Length == 8)
            return n[0] + 256 * n[1] + 65536 * n[2] + 16777216 * n[3] +
                4294967296 * n[4] + 1099511627776 * n[5] +
                281474976710656 * n[6] + 72057594037927936 * n[7];
        else
            throw new ArgumentException();
    }

    private static Int16[] separate_to_4_Int16(Int64 a)
    {
        return new []{(Int16) a, (Int16) (a >> 0x10),
            (Int16) (a >> 0x20), (Int16) (a >> 0x30)};
    }

    private static byte[] separate_to_8_byte(Int64 a)
    {
        return new []{(byte) a, (byte) (a >> 8), (byte) (a >> 0x10),
            (byte) (a >> 0x18), (byte) (a >> 0x20), (byte) (a >> 0x28),
            (byte) (a >> 0x30), (byte) (a >> 0x38)};
    }

    public static string to_string(byte[] values)
    {
        string s = "";
        foreach (var val in values)
            s += (s == "" ? "" : " ") + val;
        return s;
    }

    public static string to_string(Int16[] values)
    {
        string s = "";
        foreach (var val in values)
            s += (s == "" ? "" : " ") + val;
        return s;
    }

    // ---

    public static Int64 deconvolution(Int64 a, Int64 b)
    {
        var large = System.Numerics.BigInteger.Parse(
            "1000000000000000000000000000000");
        return (Int64)(((
            (new System.Numerics.BigInteger(a) * large)
            / new System.Numerics.BigInteger(b))
            * 72057594037927936) / large);
    }
}

Try this code at: rextester.com/YPKFA14408

The output is:

The convolution of (1 2 3 4 5 6 7 8) with (1 2 3 4 5 6 7 8) is: (1 4 10 20 35 56 84 120)

The deconvolution of (1 4 10 20 35 56 84 120) with (1 2 3 4 5 6 7 8) is: (219 251 194 172 10 94 253 14)

A demonstration that the first line of the output is correct can be found here: https://oeis.org/A000292

"Online Encyclopedia of Integer Sequences" … "Tetrahedral (or triangular pyramidal) numbers" … "the convolution of the natural numbers with themselves. - Felix Goldberg"


Q: What is wrong with the deconvolution function in the code provided, or with my understanding of what the result should be?

1

There are 1 answers

7
Ben Voigt On BEST ANSWER

The difference is carries.

Treat the inputs as polynomials...

00000101 becomes x^2 + x^0
00000111 becomes x^2 + x^1 + x^0

Polynomial multiplication and convolution are exactly the same operation (let x be z-1, the z-transform primitive for a delay... now multiplication is taking place in the transform domain after all)

But arithmetic multiplication is not the same as polynomial multiplication, because of carries. The result of

(x^1 + x^0) * (x^1 + x^0)

is

x^2 + 2 x^1 + x^0

and indeed the convolution

conv([1 1], [1 1]) = [1 2 1]

but although

0011 * 0011

is similarly

0100 + 2 * 0010 + 0001

in arithmetic, 2 * 0010 becomes 0100, and the final result of

1001

is totally different from convolution (polynomial multiplication), because of the carrying that took place.

The arithmetic result only holds when x = 2... but in convolution x must be the time delay operator.

The same effect occurs if you are working with groups of 8 bits each representing the amplitude of one time sample. Arithmetic within a group is ok, but as soon as you have a carry across a group boundary, the equivalence between convolution and multiplication is lost.