Convolution with decimation!=1 accesses outside input

Key: VSIPL1616

Legacy Issue Number: 18211

Status: open

Source: dpdx.net ( Brooks Moses)

Summary:
(Copied from an internal Mentor Graphics issue.)
Introduction
CVSIPL defines the following equation to compute the values y_n of a
minimum support 1D convolution:y_n = Sum k from 0 to M1 of ( h_k * x_(n*D + (M1)  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)CVSIPL defines the length of the output vector (i.e. values n for
which y_n is defined) as:floor( (N1) / D )  floor( (M1) / 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((51)/1)  floor((41)/1) + 1
= 4  3 + 1
= 2Consider 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_1In computing y_0 and y_1, only values
{ x_j  0 <= j < 5 }are required.
Next, consider the case:
N = 5, M = 4, D = 2.
The expected output length is:
floor((51)/2)  floor((41)/2) + 1
= floor(4/2)  floor(3/2) + 1
= 2  1 + 1
= 2To compute y_n, we need to multiply h_k by x_(n D + (M1)  k).
For n=1, k=0, the index to x isn D + (M1)  k
= 1 2 + (41)  0
= 2 + 3  0
= 5However, 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
{n}
need to determine the range of indicesfor 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 <= M1:0 <= n*D + (M1)  k < N
The lower bound is true for all n >= 0 since n*D >= 0 and (M1)k >= 0.
The upper bound:
n*D + (M1)  k < N
The largest LHS occurs when k = 0:
n*D + (M1) < N
n*D < N  M + 1
n < ceil((NM+1) / D)
The largest n that we can compute y_n for is
n = ceil((NM+1) / D)1
Hence the output length is:
ceil((NM+1) / D)
In the case where D=1, this reduces to:
NM+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=1Equation 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 = M1)
n*D  (M1) < N
n*D < N + (M1)
n < ceil((N+M1)/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  First, when D=2, a value (x_5) outside the support of x is required

Reported: VSIPL++ 1.2 — Tue, 23 Oct 2012 04:00 GMT

Updated: Fri, 6 Mar 2015 20:57 GMT