VSIPL16 — Convolution with decimation!=1 accesses outside input

(Copied from an internal Mentor Graphics issue.)

Introduction
------------

C-VSIPL defines the following equation to compute the values y_n of a
minimum support 1-D convolution:

y_n = Sum k from 0 to M-1 of ( h_k * x_(n*D + (M-1) - k) )

Where:

N is the input vector length,
M is the kernel length,
D is the deimation,
x_j is the input vector (of size N)
h_k is the kernel vector (of size M)

C-VSIPL defines the length of the output vector (i.e. values n for
which y_n is defined) as:

floor( (N-1) / D ) - floor( (M-1) / D) + 1

In some cases, this equation requires the implementation to access
elements of x outside the x proper.

Example Illustrating the Problem
--------------------------------

To see this, first consider the case (where the equation is correct):

N = 5, M = 4, D = 1.

The expected output length is:

floor((5-1)/1) - floor((4-1)/1) + 1
= 4 - 3 + 1
= 2

Consider a diagram showing how each output is computed:

x_0 x_1 x_2 x_3 x_4

y_0 = h_3 h_2 h_1 h_0

y_1 = h_3 h_2 h_1 h_0

I.e.
y_0 = h_0*x_3 + h_1*x_2 + h_2*x_1 + h_3*x_0
y_1 = h_0*x_4 + h_1*x_3 + h_2*x_2 + h_3*x_1

In computing y_0 and y_1, only values

are required.

Next, consider the case:

N = 5, M = 4, D = 2.

The expected output length is:

floor((5-1)/2) - floor((4-1)/2) + 1
= floor(4/2) - floor(3/2) + 1
= 2 - 1 + 1
= 2

To compute y_n, we need to multiply h_k by x_(n D + (M-1) - k).
For n=1, k=0, the index to x is

n D + (M-1) - k
= 1 2 + (4-1) - 0
= 2 + 3 - 0
= 5

However, x_5 is not a valid element of x.

Pictorally:

x_0 x_1 x_2 x_3 x_4

y_0 = h_3 h_2 h_1 h_0

y_1 = h_3 h_2 h_1 h_0

Finally, consider two more cases:

N = 5, M = 4, D = 3 -> output length = 1

N = 5, M = 4, D = 4 -> output length = 2

(Why should D=3 be any different from the other cases?)

This does not make sense for several reasons:

• First, when D=2, a value (x_5) outside the support of x is required
  to compute the result. This might run counter to user's expectations
  of how the minimal output region of support should be defined.

• Second, changing the decimation

A New Equation for the Output Length
------------------------------------

What should the output length be?

Starting from the assumption that the formula for y_n is correct, we
need to determine the range of indices

for which all indices to x
are valid.

In particular, we want to find the largest n such that the following
holds for all k in 0 <= k <= M-1:

0 <= n*D + (M-1) - k < N

The lower bound is true for all n >= 0 since n*D >= 0 and (M-1)-k >= 0.

The upper bound:

n*D + (M-1) - k < N

The largest LHS occurs when k = 0:

n*D + (M-1) < N

n*D < N - M + 1

n < ceil((N-M+1) / D)

The largest n that we can compute y_n for is

n = ceil((N-M+1) / D)-1

Hence the output length is:

ceil((N-M+1) / D)

In the case where D=1, this reduces to:

N-M+1

(which is what the current equation reduces too when D=1).

Comparing the new equation with the old one for the previous examples:

N=5, M=4, D=1 old=2 new=2
N=5, M=4, D=2 old=2 new=1
N=5, M=4, D=3 old=1 new=1
N=5, M=4, D=4 old=2 new=1

Equation for Full Support Output Length
---------------------------------------

Similarly, the lengths for full support case are also incorrect.

For full support, the assumption is that each output vector should
have at least one value from x used to compute it:

n*D - k < N

(minimal index occurs when k = M-1)

n*D - (M-1) < N

n*D < N + (M-1)

n < ceil((N+M-1)/D)

However, here the error is less of an issue because the equation to
compute y_n is based on < x_k > rather than x_k. <x_k> is defined to
be 0 if k < 0 of k >= N