(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
{ x_j | 0 <= j < 5 }
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
{n}
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
conversions are defensible, and even