Repository: scipy/scipy
Pull request: BUG: signal: fixed issue #10360
#### Reference issue
Fixes #10360

#### What does this implement/fix?
This fixes few bugs in firwin2 function:

1. Ensures that freq array is not modified and can be passed as tuple.
2. Adds additional requirements for freq and raises an exception when they are not met.
3. Multiplies eps by Nyquist  frequency.

I also added two tests as described here #10360. 



Problem statement:
firwin2 inadvertently modifies input and may result in undefined behavior for nyq>1
<!-- 
Thank you for taking the time to report a SciPy issue.
Please describe the issue in detail, and for bug reports
fill in the fields below. You can delete the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
-->

I tried to adapt a sample code from the documentation of `scipy.signal.firls` and found that preprocessing in `scipy.signal.firwin2` has unwanted effects. In addition, the implementation potentially misses the special case of discontinuous gain (allowed in the docs) if nyq≫1.

### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be marked as a code block automatically 
-->


```
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
fig, axs = plt.subplots(2)
fs = 1000.0  # Hz
desired = (0, 0, 1, 1, 0, 0)
for bi, bands in enumerate(((0, 100, 200, 300, 400, 500), (0, 200, 200, 300, 300, 500))):
    ## the following would be required for producing a figure:
    ## for bi, bands in enumerate(((0, 100, 200, 300, 400, 500), [0, 200, 200, 300, 300, 500])):
    fir_firwin2 = signal.firwin2(73, bands, desired, fs=fs)
    hs = list()
    ax = axs[bi]
    freq, response = signal.freqz(fir_firwin2)
    hs.append(ax.semilogy(0.5*fs*freq/np.pi, np.abs(response))[0])
    for band, gains in zip(zip(bands[::2], bands[1::2]),
                           zip(desired[::2], desired[1::2])):
        ax.semilogy(band, np.maximum(gains, 1e-7), 'k--', linewidth=2)
    if bi == 0:
        ax.legend(hs, ('firwin2', ), loc='lower center', frameon=False)
    else:
        ax.set_xlabel('Frequency (Hz)')
    ax.grid(True)
    ax.set(title='Band-pass %d-%d Hz' % tuple(bands[2:4]), ylabel='Magnitude')

fig.tight_layout()
plt.show()
```

### Error message:
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```
Traceback (most recent call last):
  File "examples/filter-design-error.py", line 8, in <module>
    fir_firwin2 = signal.firwin2(73, bands, desired, fs=fs)
  File "C:\ProgramData\Anaconda3\envs\py36\lib\site-packages\scipy\signal\fir_filter_design.py", line 604, in firwin2
    freq[k] = freq[k] - eps
TypeError: 'tuple' object does not support item assignment
```

Digging into the code showed that the problematic part modifies the _array-like_ `freq`:
```
    # Tweak any repeated values in freq so that interp works.
    eps = np.finfo(float).eps
    for k in range(len(freq)):
        if k < len(freq) - 1 and freq[k] == freq[k + 1]:
            freq[k] = freq[k] - eps
            freq[k + 1] = freq[k + 1] + eps

    # Linearly interpolate the desired response on a uniform mesh `x`.
    x = np.linspace(0.0, nyq, nfreqs)
    fx = np.interp(x, freq, gain)
```

This can be avoided by making a copy of freq:
**Solution 1&2 part A**
```
    d = np.diff(freq)
    if (d == 0).any():
        freq = np.array(freq, copy=True)
```

Thinking further it turned out that adding and subtracting `eps` is effective only if `freq[k]<=1`, i.e., `nyq==1`. Therefore either of the following

**Solution 1 part B**
```
    # Tweak any repeated values in freq so that interp works.
    eps = np.finfo(float).eps
    freq = freq / nyq
    for k in range(len(freq)):
        if k < len(freq) - 1 and freq[k] == freq[k + 1]:
            freq[k] = freq[k] - eps
            freq[k + 1] = freq[k + 1] + eps

    # Linearly interpolate the desired response on a uniform mesh `x`.
    x = np.linspace(0.0, 1.0, nfreqs)  # nyq is now at freq == 1.0
    fx = np.interp(x, freq, gain)
```

or

**Solution 2 part B**
```
    # Tweak any repeated values in freq so that interp works.
    eps = np.finfo(float).eps * nyq
    for k in range(len(freq)):
        if k < len(freq) - 1 and freq[k] == freq[k + 1]:
            freq[k] = freq[k] - eps
            freq[k + 1] = freq[k + 1] + eps

    # Linearly interpolate the desired response on a uniform mesh `x`.
    x = np.linspace(0.0, nyq, nfreqs)
    fx = np.interp(x, freq, gain)
```

Please confirm the above issue and suggest which one to choose between Solution 1 and 2. Then I will be happy to submit a pull request.


### Scipy/Numpy/Python version information:
<!-- You can simply run the following and paste the result in a code block
```
import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
```
-->
1.2.1 1.16.2 sys.version_info(major=3, minor=6, micro=8, releaselevel='final', serial=0)

Hints from issue discussion:
This does indeed look like a bug. I would multiply `eps` by `nyq` as in 2B. And if you can add two tests that fail in `master`:

1. `tuple` input works
2. `ndarray` input is not modified by the function

and pass on your PR, that would be great!
Whoops, my bugs.  A pull request would be great.

An alternative to `eps = np.finfo(float).eps * nyq` is to make `eps` depend on `freqs[k]` by making it a multiple of `numpy.spacing(freqs[k])`.  Something like this:

    # Tweak any repeated values in freq so that interp works.
    for k in range(len(freq)):
        if k < len(freq) - 1 and freq[k] == freq[k + 1]:
            eps = M * np.spacing(freq[k])
            freq[k] = freq[k] - eps
            freq[k + 1] = freq[k + 1] + eps

where M is some fixed constant, say 10. 

There are more bugs lurking in that code (what was I thinking?):
* What if 0 is repeated?  Then the first frequency is changed to `-eps`.  And what if the Nyquist frequency is repeated?  Then the last frequency is changed to `nyq + eps`.  These adjustments violate the assumption that the frequencies are between 0 and `nyq`.
* What if `freq` contains the sequence `f, f, f + tiny`, where `f` is some constant and `tiny` is less than or equal to the `eps` that is used to adjust the repeated frequency?  After the adjustment, we have `f - eps, f + eps, f + tiny`, and now those second two are out of order, or repeated if `eps == tiny`.

The repeated 0 or Nyquist frequencies could be handled by adjusting just one of the repeats, i.e. `[0, 0, ...]` -> `[0, eps, ...]`, and `[..., nyq, nyq]` -> `[..., nyq-eps, nyq]`.  But is there really any need to allow either of these frequencies to be repeated?  I don't see how a discontinuity at the end of the interval has any practical use.

Alternatively, and to handle the second lurking bug mentioned above,  a simple way to handle these cases is to check the bounds and strict monotonicity requirement again after making the adustments, and if the requirements are not satisfied, give up and raise an exception, preferably with an error message that is as informative as possible.  E.g. "the frequency 0 must not be repeated", "the Nyquist frequency must not be repeated", "a third frequency is too close to a repeated frequency" (or something like that).
@stippingerm would you be up for making a pull request to implement @WarrenWeckesser's last suggestion:

> Alternatively, and to handle the second lurking bug mentioned above, a simple way to handle these cases is to check the bounds and strict monotonicity requirement again after making the adustments, and if the requirements are not satisfied, give up and raise an exception, preferably with an error message that is as informative as possible. E.g. "the frequency 0 must not be repeated", "the Nyquist frequency must not be repeated", "a third frequency is too close to a repeated frequency" (or something like that).
Hey, if no one is working on this issue, I'd like to give it a shot. @larsoner could I make a PR with @WarrenWeckesser 's suggestion?
@grassknoted, I don't see a response from @stippingerm to @larsoner's question, so I'd say go for it.

Commit message: DOC: add new requirements for input

Patch to review:
diff --git a/scipy/signal/fir_filter_design.py b/scipy/signal/fir_filter_design.py
index d84e4aedd091..99c69015b925 100644
--- a/scipy/signal/fir_filter_design.py
+++ b/scipy/signal/fir_filter_design.py
@@ -505,7 +505,8 @@ def firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=None,
         Nyquist.  The Nyquist frequency is half `fs`.
         The values in `freq` must be nondecreasing.  A value can be repeated
         once to implement a discontinuity.  The first value in `freq` must
-        be 0, and the last value must be ``fs/2``.
+        be 0, and the last value must be ``fs/2``. Values 0 and ``fs/2`` must
+        not be repeated.
     gain : array_like
         The filter gains at the frequency sampling points. Certain
         constraints to gain values, depending on the filter type, are applied,
@@ -607,6 +608,10 @@ def firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=None,
     d2 = d[:-1] + d[1:]
     if (d2 == 0).any():
         raise ValueError('A value in freq must not occur more than twice.')
+    if freq[1] == 0:
+        raise ValueError('Value 0 must not be repeated in freq')
+    if freq[-2] == nyq:
+        raise ValueError('Value fs/2 must not be repeated in freq')
 
     if antisymmetric:
         if numtaps % 2 == 0:
@@ -632,13 +637,23 @@ def firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=None,
     if nfreqs is None:
         nfreqs = 1 + 2 ** int(ceil(log(numtaps, 2)))
 
+    d = np.diff(freq)
+    if (d == 0).any():
+        freq = np.array(freq, copy=True)
+
     # Tweak any repeated values in freq so that interp works.
-    eps = np.finfo(float).eps
+    eps = np.finfo(float).eps * nyq
     for k in range(len(freq)):
         if k < len(freq) - 1 and freq[k] == freq[k + 1]:
             freq[k] = freq[k] - eps
             freq[k + 1] = freq[k + 1] + eps
 
+    # Check if freq is strictly increasing after tweak
+    d = np.diff(freq)
+    if (d <= 0).any():
+        raise ValueError("Freq cannot contain number that is too "
+                         "close to some repeated value")
+
     # Linearly interpolate the desired response on a uniform mesh `x`.
     x = np.linspace(0.0, nyq, nfreqs)
     fx = np.interp(x, freq, gain)
diff --git a/scipy/signal/tests/test_fir_filter_design.py b/scipy/signal/tests/test_fir_filter_design.py
index 97baa76791af..c57443aa8ea7 100644
--- a/scipy/signal/tests/test_fir_filter_design.py
+++ b/scipy/signal/tests/test_fir_filter_design.py
@@ -392,6 +392,18 @@ def test_fs_nyq(self):
         taps2 = firwin2(80, [0.0, 30.0, 60.0], [1.0, 1.0, 0.0], nyq=60.0)
         assert_array_almost_equal(taps1, taps2)
 
+    def test_tuple(self):
+        taps1 = firwin2(150, (0.0, 0.5, 0.5, 1.0), (1.0, 1.0, 0.0, 0.0))
+        taps2 = firwin2(150, [0.0, 0.5, 0.5, 1.0], [1.0, 1.0, 0.0, 0.0])
+        assert_array_almost_equal(taps1, taps2)
+
+    def test_input_modyfication(self):
+        freq1 = np.array([0.0, 0.5, 0.5, 1.0])
+        freq2 = np.array(freq1)
+        firwin2(80, freq1, [1.0, 1.0, 0.0, 0.0])
+        assert_equal(freq1, freq2)
+
+
 class TestRemez(object):
 
     def test_bad_args(self):


Code review comments:
diff_hunk: @@ -632,13 +637,23 @@ def firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=None,
     if nfreqs is None:
         nfreqs = 1 + 2 ** int(ceil(log(numtaps, 2)))
 
+    d = np.diff(freq)
+    if (d == 0).any():
+        freq = np.array(freq, copy=True)
+
     # Tweak any repeated values in freq so that interp works.
-    eps = np.finfo(float).eps
+    eps = np.finfo(float).eps * nyq
     for k in range(len(freq)):
         if k < len(freq) - 1 and freq[k] == freq[k + 1]:
             freq[k] = freq[k] - eps
             freq[k + 1] = freq[k + 1] + eps
 
+    # Check if freq is strictly increasing after tweak
+    d = np.diff(freq)
+    if (d <= 0).any():
+        raise ValueError("Freq cannot contain number that is too "
+                         "close to some repeated value")
path: scipy/signal/fir_filter_design.py
text: @user1:
"freq cannot contain numbers that are too close (within eps * nyquist: %s) to a repeated value" or so