Repository: scipy/scipy
Pull request: BUG: Fix length calculation in upfirdn
Closes #9844.

Implements the suggestion implemented in that issue by @grlee77 and @sirtom67, except I kept `ceil` instead of floor as now the length matches both the naive `upfirdn` calculation with `np.convolve` and slicing, as well as MATLAB outputs (calculated and put in the test).

Also refactors a few tests to use `pytest.mark.parametrize` instead of `for` in the tests (was useful for figuring out what was actually wrong) and fixes a few flake8 errors.

Would be good to get this into 1.5.0 if possible. @grlee77 @sirtom67 would either of you have time to look in the next few days? Sorry for the short notice...

Problem statement:
scipy.signal.upfirdn gives different length matrix versus MATLAB upfirdn
This is something I am not sure is a bug, a mistake in my code, or just a difference in design or implementation.

Now everything before my call to upfirdn looks to be the same in both languages, floating point precision notwithstanding. So I get to this call:

MATLAB:
>y = upfirdn(x,h,p,q);

Python:
>y = upfirdn(h, x, p, q)

x, h, p, q, are the same across both files.
MATLAB upfirdn gives me a matrix of 1x647.
SciPy upfirdn gives me a matrix of 1x651.

Here are the first and last ten values from MATLAB:

> 0	-0.000199467579547424	-0.000728046369939899	-0.00110566784707238	-0.000996359004971441	-0.000228970643181423	0.000615489415076301	0.00126315546064305	0.00142874951087505	0.000908693634358364
>.....
>4.03845240072376e-06	1.87473690616758e-05	2.31142819067192e-05	1.71887339187862e-05	5.37933346930270e-06	-9.91074437326740e-06	-1.88844901213745e-05	-1.71623305506084e-05	-9.10662941865632e-06	-5.37239384457409e-07

Now here are the first ten values from Python:
>0.00000000e+00 -1.97601366e-04 -7.21234789e-04 -1.09532325e-03
  -9.87037098e-04 -2.26828400e-04  6.09910612e-04  1.25175809e-03
   1.41591478e-03  9.00565745e-04

Now I think the differences there are due to floating point precision. 
Here are the last 14 of values from Python:
>4.03214296e-06  1.87159549e-05  2.30728961e-05
   1.71559510e-05  5.36843263e-06 -9.88322175e-06 -1.88320470e-05
  -1.71146699e-05 -9.08133986e-06 -5.35747444e-07  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00

So you can see that the first value in the last 14 values in Python... it matches up roughly with the last 10 values in MATLAB. Except I have four extra zeros being added to the end of the list. So basically my SciPy upfirdn is adding padding at the end.

Now is this just me misunderstanding the differences in implementation or is there something more at play here?

Could someone explain to me the differences between the MATLAB implementation and the SciPy one??

Also - bonus question... should I be worried by the differences in the output values regardless of this padding issue?

Thank you!

Hints from issue discussion:
So right now I'm just deleting those trailing zeros. It seems like in all the cases I've tried, I always get four of them. I don't know where they're coming from, but I'm using np.delete to remove them. Best I've got right now. I'd love to understand why there are differences between the Scipy and MATLAB methods.
I would plot the two results, it's possible `matlab` and `scipy.signal` are making different assumptions about what the output mode should be, e.g. output should be the same length as the input if `p=q`, or it should be the full length `len(x) + len(h) - 1`, etc. Then you can check the MATLAB source to see how it computes the output size. In SciPy you can check to see how it's computed here

https://github.com/scipy/scipy/blob/v1.2.1/scipy/signal/_upfirdn_apply.pyx#L62-L74
@DrewTNBD did you end up figuring out where the differences come from? Is the difference acceptable?
According to the [online Matlab documentation](https://www.mathworks.com/help/signal/ref/upfirdn.html), the output size returned by their `upfirdn` is:
`ceil(((len(x) - 1) * p + len(h)) / q)` where `p` is the upsampling factor and `q` the downsampling factor.

I tried a range of up/down factors and signal/filter lengths and this often gives the same length as SciPy, but sometimes the Matlab length is slightly shorter than the SciPy one. 


The length used in SciPy's case is equivalent to:
`ceil((len(x) * p + len(h) - 1) / q)`
while Matlab's is:
`ceil(((len(x) - 1) * p + len(h)) / q)`

So the difference is where the `- 1` is located in the above formula. 

Note:
The SciPy convention here is what I would have naively expected:  FIR filtering of a signal upsampled by `p` would give a signal of length `len(x) * p + len(h) - 1` which is then downsampled by `q`. If the final samples really are always zero in the case where SciPy is longer than Matlab, then it may be we are keeping some extra samples that aren't strictly necessary.

Also, when `p=1` Matlab and SciPy always agree on the length.
> The length used in SciPy's case is equivalent to:
> ceil((len(x) * p + len(h) - 1) / q)
> while Matlab's is:
> ceil(((len(x) - 1) * p + len(h)) / q)

It looks like we probably should have the ` - 1` in there. @sirtom67 do you agree?
Sorry for the late response. Thanks for tagging me Eric. I'll give it some thought. I think we need to be careful about doing it just because MATLAB does it. But it might be the right thing to do. 
OK I rederived things and since Scipy's upfirdn is always starting at sample 0, I agree with the MATLAB expression for the output length.  HOWEVER I think it should be floor, not ceil.
Specifically, I think it should be 
  floor(((len(x) - 1) * p + len(h)) / q)

Actually, seeing something strange for length(x)=1:
In [28]: upfirdn([1,1],[1],1,1)
Out[28]: array([ 1.,  0.])

In [29]: upfirdn([1,1],[0,1],1,1)
Out[29]: array([ 0.,  1.,  1.])

So it looks like there are two bugs to be fixed:

1. Behavior at short lengths (`upfirdn([1, 1], [1], 1, 1)` should be `[1, 1]` but we get `[0, 1, 1]`)
2. Length calculation is off and should use `floor(((len(x) - 1) * p + len(h)) / q)`

One question -- why prefer `floor` to `ceil` here? It would be nice to be consistent with MATLAB (`ceil`) but if their formula is less accurate then it makes sense to deviate (`floor`).

@sirtom67 do you want to look into the code? If not I (or someone else) can find time at some point.

Commit message: BUG: Better protection

Patch to review:
diff --git a/scipy/signal/_upfirdn.py b/scipy/signal/_upfirdn.py
index 578fb1772415..c449a40ed6ba 100644
--- a/scipy/signal/_upfirdn.py
+++ b/scipy/signal/_upfirdn.py
@@ -70,8 +70,9 @@ def _check_mode(mode):
 
 
 class _UpFIRDn(object):
+    """Helper for resampling."""
+
     def __init__(self, h, x_dtype, up, down):
-        """Helper for resampling"""
         h = np.asarray(h)
         if h.ndim != 1 or h.size == 0:
             raise ValueError('h must be 1-D with non-zero length')
@@ -84,10 +85,11 @@ def __init__(self, h, x_dtype, up, down):
         # This both transposes, and "flips" each phase for filtering
         self._h_trans_flip = _pad_h(h, self._up)
         self._h_trans_flip = np.ascontiguousarray(self._h_trans_flip)
+        self._h_len_orig = len(h)
 
     def apply_filter(self, x, axis=-1, mode='constant', cval=0):
-        """Apply the prepared filter to the specified axis of a N-D signal x"""
-        output_len = _output_len(len(self._h_trans_flip), x.shape[axis],
+        """Apply the prepared filter to the specified axis of N-D signal x."""
+        output_len = _output_len(self._h_len_orig, x.shape[axis],
                                  self._up, self._down)
         # Explicit use of np.int64 for output_shape dtype avoids OverflowError
         # when allocating large array on platforms where np.int_ is 32 bits
@@ -103,7 +105,7 @@ def apply_filter(self, x, axis=-1, mode='constant', cval=0):
 
 
 def upfirdn(h, x, up=1, down=1, axis=-1, mode='constant', cval=0):
-    """Upsample, FIR filter, and downsample
+    """Upsample, FIR filter, and downsample.
 
     Parameters
     ----------
@@ -147,9 +149,6 @@ def upfirdn(h, x, up=1, down=1, axis=-1, mode='constant', cval=0):
     The algorithm is an implementation of the block diagram shown on page 129
     of the Vaidyanathan text [1]_ (Figure 4.3-8d).
 
-    .. [1] P. P. Vaidyanathan, Multirate Systems and Filter Banks,
-       Prentice Hall, 1993.
-
     The direct approach of upsampling by factor of P with zero insertion,
     FIR filtering of length ``N``, and downsampling by factor of Q is
     O(N*Q) per output sample. The polyphase implementation used here is
@@ -157,6 +156,11 @@ def upfirdn(h, x, up=1, down=1, axis=-1, mode='constant', cval=0):
 
     .. versionadded:: 0.18
 
+    References
+    ----------
+    .. [1] P. P. Vaidyanathan, Multirate Systems and Filter Banks,
+           Prentice Hall, 1993.
+
     Examples
     --------
     Simple operations:
@@ -204,7 +208,6 @@ def upfirdn(h, x, up=1, down=1, axis=-1, mode='constant', cval=0):
            [ 4.,  5.],
            [ 6.,  7.],
            [ 6.,  7.]])
-
     """
     x = np.asarray(x)
     ufd = _UpFIRDn(h, x.dtype, up, down)
diff --git a/scipy/signal/_upfirdn_apply.pyx b/scipy/signal/_upfirdn_apply.pyx
index 83860c8c8262..415ebb554c5c 100644
--- a/scipy/signal/_upfirdn_apply.pyx
+++ b/scipy/signal/_upfirdn_apply.pyx
@@ -63,14 +63,15 @@ def _output_len(np.intp_t len_h,
                 np.intp_t up,
                 np.intp_t down):
     """The output length that results from a given input"""
-    cdef np.intp_t nt
-    cdef np.intp_t in_len_copy
-    in_len_copy = in_len + (len_h + (-len_h % up)) // up - 1
-    nt = in_len_copy * up
-    cdef np.intp_t need = nt // down
-    if nt % down > 0:
-        need += 1
-    return need
+    return _c_output_len(len_h, in_len, up, down)
+
+
+cdef np.intp_t _c_output_len(np.intp_t len_h,
+                             np.intp_t in_len,
+                             np.intp_t up,
+                             np.intp_t down) nogil:
+    # ceil(((in_len - 1) * up + len_h) / down), but using integer arithmetic
+    return (((in_len - 1) * up + len_h) - 1) // down + 1
 
 
 # Signal extension modes
@@ -288,6 +289,7 @@ def _apply(np.ndarray data, DTYPE_t [::1] h_trans_flip, np.ndarray out,
     cdef DTYPE_t *filter_ptr
     cdef DTYPE_t *out_ptr
     cdef int retval
+    cdef np.intp_t len_out = out.shape[axis]
 
     data_info.ndim = data.ndim
     data_info.strides = <np.intp_t *> data.strides
@@ -305,7 +307,7 @@ def _apply(np.ndarray data, DTYPE_t [::1] h_trans_flip, np.ndarray out,
         retval = _apply_axis_inner(data_ptr, data_info,
                                    filter_ptr, len_h,
                                    out_ptr, output_info,
-                                   up, down, axis, <MODE>mode, cval)
+                                   up, down, axis, <MODE>mode, cval, len_out)
     if retval == 1:
         raise ValueError("failure in _apply_axis_inner: data and output arrays"
                          " must have the same number of dimensions.")
@@ -325,7 +327,8 @@ cdef int _apply_axis_inner(DTYPE_t* data, ArrayInfo data_info,
                            DTYPE_t* h_trans_flip, np.intp_t len_h,
                            DTYPE_t* output, ArrayInfo output_info,
                            np.intp_t up, np.intp_t down,
-                           np.intp_t axis, MODE mode, DTYPE_t cval) nogil:
+                           np.intp_t axis, MODE mode, DTYPE_t cval,
+                           np.intp_t len_out) nogil:
     cdef np.intp_t i
     cdef np.intp_t num_loops = 1
     cdef bint make_temp_data, make_temp_output
@@ -403,7 +406,8 @@ cdef int _apply_axis_inner(DTYPE_t* data, ArrayInfo data_info,
 
         # call 1D upfirdn
         _apply_impl(data_row, data_info.shape[axis],
-                    h_trans_flip, len_h, output_row, up, down, mode, cval)
+                    h_trans_flip, len_h, output_row, up, down, mode, cval,
+                    len_out)
 
         # Copy from temporary output if necessary
         if make_temp_output:
@@ -423,8 +427,8 @@ cdef int _apply_axis_inner(DTYPE_t* data, ArrayInfo data_info,
 cdef void _apply_impl(DTYPE_t *x, np.intp_t len_x, DTYPE_t *h_trans_flip,
                       np.intp_t len_h, DTYPE_t *out,
                       np.intp_t up, np.intp_t down, MODE mode,
-                      DTYPE_t cval) nogil:
-    cdef np.intp_t h_per_phase = len_h / up
+                      DTYPE_t cval, np.intp_t len_out) nogil:
+    cdef np.intp_t h_per_phase = len_h // up
     cdef np.intp_t padded_len = len_x + h_per_phase - 1
     cdef np.intp_t x_idx = 0
     cdef np.intp_t y_idx = 0
@@ -435,6 +439,8 @@ cdef void _apply_impl(DTYPE_t *x, np.intp_t len_x, DTYPE_t *h_trans_flip,
     cdef bint zpad
 
     zpad = (mode == MODE_CONSTANT and cval == 0)
+    if len_out == 0:
+        return
 
     while x_idx < len_x:
         h_idx = t * h_per_phase
@@ -453,8 +459,10 @@ cdef void _apply_impl(DTYPE_t *x, np.intp_t len_x, DTYPE_t *h_trans_flip,
             h_idx += 1
         # store and increment
         y_idx += 1
+        if y_idx >= len_out:
+            return
         t += down
-        x_idx += t / up  # integer div
+        x_idx += t // up
         # which phase of the filter to use
         t = t % up
 
@@ -472,6 +480,8 @@ cdef void _apply_impl(DTYPE_t *x, np.intp_t len_x, DTYPE_t *h_trans_flip,
             out[y_idx] += xval * h_trans_flip[h_idx]
             h_idx += 1
         y_idx += 1
+        if y_idx >= len_out:
+            return
         t += down
         x_idx += t / up  # integer div
         t = t % up
diff --git a/scipy/signal/tests/test_upfirdn.py b/scipy/signal/tests/test_upfirdn.py
index 218a0062ddf5..e7fc352d4170 100644
--- a/scipy/signal/tests/test_upfirdn.py
+++ b/scipy/signal/tests/test_upfirdn.py
@@ -39,13 +39,13 @@
 from pytest import raises as assert_raises
 import pytest
 
-from scipy.signal import upfirdn, firwin, lfilter
+from scipy.signal import upfirdn, firwin
 from scipy.signal._upfirdn import _output_len, _upfirdn_modes
 from scipy.signal._upfirdn_apply import _pad_test
 
 
 def upfirdn_naive(x, h, up=1, down=1):
-    """Naive upfirdn processing in Python
+    """Naive upfirdn processing in Python.
 
     Note: arg order (x, h) differs to facilitate apply_along_axis use.
     """
@@ -91,7 +91,11 @@ def __call__(self):
     def scrub(self, x, axis=-1):
         yr = np.apply_along_axis(upfirdn_naive, axis, x,
                                  self.h, self.up, self.down)
+        want_len = _output_len(len(self.h), x.shape[axis], self.up, self.down)
+        assert yr.shape[axis] == want_len
         y = upfirdn(self.h, x, self.up, self.down, axis=axis)
+        assert y.shape[axis] == want_len
+        assert y.shape == yr.shape
         dtypes = (self.h.dtype, x.dtype)
         if all(d == np.complex64 for d in dtypes):
             assert_equal(y.dtype, np.complex64)
@@ -106,6 +110,9 @@ def scrub(self, x, axis=-1):
         assert_allclose(yr, y)
 
 
+_UPFIRDN_TYPES = (int, np.float32, np.complex64, float, complex)
+
+
 class TestUpfirdn(object):
 
     def test_valid_input(self):
@@ -113,42 +120,77 @@ def test_valid_input(self):
         assert_raises(ValueError, upfirdn, [], [1], 1, 1)  # h.ndim != 1
         assert_raises(ValueError, upfirdn, [[1]], [1], 1, 1)
 
-    def test_vs_lfilter(self):
-        # Check that up=1.0 gives same answer as lfilter + slicing
+    @pytest.mark.parametrize('len_h', [1, 2, 3, 4, 5])
+    @pytest.mark.parametrize('len_x', [1, 2, 3, 4, 5])
+    def test_singleton(self, len_h, len_x):
+        # gh-9844: lengths producing expected outputs
+        h = np.zeros(len_h)
+        h[len_h // 2] = 1.  # make h a delta
+        x = np.ones(len_x)
+        y = upfirdn(h, x, 1, 1)
+        want = np.pad(x, (len_h // 2, (len_h - 1) // 2), 'constant')
+        assert_allclose(y, want)
+
+    def test_shift_x(self):
+        # gh-9844: shifted x can change values?
+        y = upfirdn([1, 1], [1.], 1, 1)
+        assert_allclose(y, [1, 1])  # was [0, 1] in the issue
+        y = upfirdn([1, 1], [0., 1.], 1, 1)
+        assert_allclose(y, [0, 1, 1])
+
+    # A bunch of lengths/factors chosen because they exposed differences
+    # between the "old way" and new way of computing length, and then
+    # got `expected` from MATLAB
+    @pytest.mark.parametrize('len_h, len_x, up, down, expected', [
+        (2, 2, 5, 2, [1, 0, 0, 0]),
+        (2, 3, 6, 3, [1, 0, 1, 0, 1]),
+        (2, 4, 4, 3, [1, 0, 0, 0, 1]),
+        (3, 2, 6, 2, [1, 0, 0, 1, 0]),
+        (4, 11, 3, 5, [1, 0, 0, 1, 0, 0, 1]),
+    ])
+    def test_length_factors(self, len_h, len_x, up, down, expected):
+        # gh-9844: weird factors
+        h = np.zeros(len_h)
+        h[0] = 1.
+        x = np.ones(len_x)
+        y = upfirdn(h, x, up, down)
+        assert_allclose(y, expected)
+
+    @pytest.mark.parametrize('down, want_len', [  # lengths from MATLAB
+        (2, 5015),
+        (11, 912),
+        (79, 127),
+    ])
+    def test_vs_convolve(self, down, want_len):
+        # Check that up=1.0 gives same answer as convolve + slicing
         random_state = np.random.RandomState(17)
         try_types = (int, np.float32, np.complex64, float, complex)
         size = 10000
-        down_factors = [2, 11, 79]
 
         for dtype in try_types:
             x = random_state.randn(size).astype(dtype)
             if dtype in (np.complex64, np.complex128):
                 x += 1j * random_state.randn(size)
 
-            for down in down_factors:
-                h = firwin(31, 1. / down, window='hamming')
-                yl = lfilter(h, 1.0, x)[::down]
-                y = upfirdn(h, x, up=1, down=down)
-                assert_allclose(yl, y[:yl.size], atol=1e-7, rtol=1e-7)
-
-    def test_vs_naive(self):
-        tests = []
-        try_types = (int, np.float32, np.complex64, float, complex)
-
-        # Simple combinations of factors
-        for x_dtype, h in product(try_types, (1., 1j)):
-            tests.append(UpFIRDnCase(1, 1, h, x_dtype))
-            tests.append(UpFIRDnCase(2, 2, h, x_dtype))
-            tests.append(UpFIRDnCase(3, 2, h, x_dtype))
-            tests.append(UpFIRDnCase(2, 3, h, x_dtype))
+            h = firwin(31, 1. / down, window='hamming')
+            yl = upfirdn_naive(x, h, 1, down)
+            y = upfirdn(h, x, up=1, down=down)
+            assert y.shape == (want_len,)
+            assert yl.shape[0] == y.shape[0]
+            assert_allclose(yl, y, atol=1e-7, rtol=1e-7)
 
-        # mixture of big, small, and both directions (net up and net down)
-        # use all combinations of data and filter dtypes
-        factors = (100, 10)  # up/down factors
-        cases = product(factors, factors, try_types, try_types)
-        for case in cases:
-            tests += self._random_factors(*case)
+    @pytest.mark.parametrize('x_dtype', _UPFIRDN_TYPES)
+    @pytest.mark.parametrize('h', (1., 1j))
+    @pytest.mark.parametrize('up, down', [(1, 1), (2, 2), (3, 2), (2, 3)])
+    def test_vs_naive_delta(self, x_dtype, h, up, down):
+        UpFIRDnCase(up, down, h, x_dtype)()
 
+    @pytest.mark.parametrize('x_dtype', _UPFIRDN_TYPES)
+    @pytest.mark.parametrize('h_dtype', _UPFIRDN_TYPES)
+    @pytest.mark.parametrize('p_max, q_max',
+                             list(product((10, 100), (10, 100))))
+    def test_vs_naive(self, x_dtype, h_dtype, p_max, q_max):
+        tests = self._random_factors(p_max, q_max, h_dtype, x_dtype)
         for test in tests:
             test()
 


Code review comments:
diff_hunk: @@ -63,14 +63,15 @@ def _output_len(np.intp_t len_h,
                 np.intp_t up,
                 np.intp_t down):
     """The output length that results from a given input"""
-    cdef np.intp_t nt
-    cdef np.intp_t in_len_copy
-    in_len_copy = in_len + (len_h + (-len_h % up)) // up - 1
-    nt = in_len_copy * up
-    cdef np.intp_t need = nt // down
-    if nt % down > 0:
-        need += 1
-    return need
+    return _c_output_len(len_h, in_len, up, down)
+
+
+cdef np.intp_t _c_output_len(np.intp_t len_h,
path: scipy/signal/_upfirdn_apply.pyx
text: @user1:
not sure there is any benefit to splitting this `_c_output_len` out from `_output_len`? It doesn't look like any C/Cython code calls `_c_output_len` directly.

That said, I am fine with it as is.

@author:
Ahh yes that's cruft from when I used it in C code, pushing a commit to get rid of that

diff_hunk: @@ -472,6 +480,8 @@ cdef void _apply_impl(DTYPE_t *x, np.intp_t len_x, DTYPE_t *h_trans_flip,
             out[y_idx] += xval * h_trans_flip[h_idx]
             h_idx += 1
         y_idx += 1
+        if y_idx >= len_out:
+            return
         t += down
         x_idx += t / up  # integer div
path: scipy/signal/_upfirdn_apply.pyx
text: @user1:
for consistency with above

```suggestion
        x_idx += t // up
```