VSIPL 1.6 RTF Avatar
  1. OMG Issue

VSIPL16 — Convolution with decimation!=1 accesses outside input

  • Key: VSIPL16-16
  • Legacy Issue Number: 18211
  • Status: open  
  • Source: dpdx.net ( Brooks Moses)
  • Summary:

    (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

  • Reported: VSIPL++ 1.2 — Tue, 23 Oct 2012 04:00 GMT
  • Updated: Fri, 6 Mar 2015 20:57 GMT