Repository: scipy/scipy
Pull request: ENH: Use faster solvers for firls
`Q` is symmetric, so we should use that to speed up solving. Also, we can try using `linalg.solve` first, and fall back to a `pinvh` if necessary.

I'm not 100% sure that it's okay to use `solve` in cases where it warns about rank deficiency. It seems to produce okay results in my tests, but if someone thinks it's better to catch the emitted warning and use `pinvh` in that case, too, I can change it to do so.

Closes #9835.

Problem statement:
scipy.signal.firls seems to be inefficient versus MATLAB firls
Hi everyone,
I've been porting some code over from MATLAB to Python, and I have run into an issue with scipy.signal.firls. It doesn't seem to be as fast/efficient as the MATLAB equivalent. I don't know if I'm using it wrong, but I wanted to report this and maybe get some guidance if I am indeed using it wrong.

Here is how I am executing the code in MATLAB. I've hard-coded the variables here for illustration purposes:
`firls(23220, [0, 0.000861326442721792, 0.000861326442721792, 1], [1 1 0 0]);`

This takes less than a couple of seconds to process.

Now in Python I do the following:
`scipy.signal.firls(23221, [0, 0.000861326442721792, 0.000861326442721792, 1], [1, 1, 0, 0])`

First thing you'll see is that argument 1 needs to have 1 added to it to make it an odd number. This I am guessing is related to the differences between Python and Matlab indexing functionality.

But when I run this, it takes over 5 minutes to process. I am guessing it's because of the large number of taps I am looking to process. 

Now if I sit and let it process it does give me the same output as MATLAB, roughly speaking. Which makes me think that for my purposes I'm not doing anything wrong. But can you spot anything here?

FWIW - I'm porting across a MATLAB library, which required me to also port across the MATLAB Resample class, and within that sits this firls line that is very slow in Python.

Any thoughts on how I can get past this?

Hints from issue discussion:
The MATLAB documentation for their firls implementation states:

> b = firls(n,f,a) returns row vector b containing the n+1 coefficients of the order n FIR filter whose frequency-amplitude characteristics approximately match those given by vectors f and a. The output filter coefficients, or “taps,” in b obey the symmetry relation.

So is 'n' in their implementation equivalent to scipy.signal.firls(numtaps) parameter?
Can you see if #9838 solves your problem? And yes MATLAB does `n+1` and requires `n` to be odd, we instead have `ntaps` odd, and you get a `ntaps` length filter back.
@larsoner 

UBER thanks for responding so quicky. That solved MUCH quicker, and on the face of it, it looks like I am getting the correct values. However the complex numbers returned are slightly different in some cases.

For example, the return array begins with in MATLAB:

> -3.35759847762220e-20 
Which resolves to 0.

But in your solution it begins with:

> 7.632048298997159e-20 
Which resolves to -0.

I don't know if that is a problem necessarily? 

Here are the first 10 numbers in case you wanted to check them on your end. I wont spam you with the entire matrix!!

MATLAB:
> -3.35759847762220e-20	-7.41946241368262e-08	-1.48401488327497e-07	-2.22620052386448e-07	-2.96849775949765e-07	-3.71090118479317e-07	-4.45340539266368e-07	-5.19600497435580e-07	-5.93869451949603e-07	-6.68146861611719e-07

Python:

> 7.63204830e-20 -7.41946241e-08 -1.48401488e-07 -2.22620052e-07
 -2.96849776e-07 -3.71090118e-07 -4.45340539e-07 -5.19600497e-07
 -5.93869452e-07 -6.68146862e-07

The numbers are slightly different because of precision differences between MATLAB and Python, but I'm not so worried about that.

I've patched my local Scipy.signal.fir_filter_design.py file. Should I leave it like that? Will your fix make it to the next official release?

Thanks again!
I could not get `octave` and SciPy outputs to match. But I did check that the properties of the returned filters matched those of the design criteria. Can you check to see if they satisfy yours using `freqz`?
(and yes assuming the PR is merged before 1.3.0 branching date, it will be in that release)
The formatting of the complex numbers looks a bit different, but if I whack them into Complex Number Calculator, then what I see in MATLAB's 'freqz' result is the same as what I'm seeing in Python. I *think*.

Again, here are some reference numbers. The first 10 results in each.

Matlab:
>1160.99999999999 + 0.00000000000000i
0.00931142233356466 + 0.0151114824035457i
0.0123139146781979 - 0.0244635468636103i
0.0190376493996473 - 0.00163942434520292i
0.00327291561885030 + 0.00441301585403953i
0.00111064687613271 - 0.00277725458056106i
0.00460084848837975 - 0.000798322943340352i
0.00141683250770980 + 0.00160231810952528i
0.000256602799033961 - 0.000845906063621180i
0.00196579592028704 - 0.000518124162136642i

Python:
>1.16100000e+03+0.00000000e+00j  9.31142233e-03+1.51114824e-02j
  1.23139147e-02-2.44635469e-02j  1.90376494e-02-1.63942435e-03j
  3.27291562e-03+4.41301585e-03j  1.11064688e-03-2.77725458e-03j
  4.60084849e-03-7.98322943e-04j  1.41683251e-03+1.60231811e-03j
  2.56602799e-04-8.45906064e-04j  1.96579592e-03-5.18124162e-04j

You can see the formatting is a bit different, but the resolved numbers are the same.

Commit message: ENH: Use lstsq

Patch to review:
diff --git a/scipy/signal/fir_filter_design.py b/scipy/signal/fir_filter_design.py
index 23c7e0b37d48..4e8fb8d74199 100644
--- a/scipy/signal/fir_filter_design.py
+++ b/scipy/signal/fir_filter_design.py
@@ -9,7 +9,8 @@
 import numpy as np
 from numpy.fft import irfft, fft, ifft
 from scipy.special import sinc
-from scipy.linalg import toeplitz, hankel, pinv
+from scipy.linalg import (toeplitz, hankel, solve, LinAlgError, LinAlgWarning,
+                          lstsq)
 from scipy._lib.six import string_types
 
 from . import sigtools
@@ -965,6 +966,8 @@ def firls(numtaps, bands, desired, weight=None, nyq=None, fs=None):
     bands = np.asarray(bands).flatten() / nyq
     if len(bands) % 2 != 0:
         raise ValueError("bands must contain frequency pairs.")
+    if (bands < 0).any() or (bands > 1).any():
+        raise ValueError("bands must be between 0 and 1 relative to Nyquist")
     bands.shape = (-1, 2)
 
     # check remaining params
@@ -1032,8 +1035,20 @@ def firls(numtaps, bands, desired, weight=None, nyq=None, fs=None):
     b[1:] += m * np.cos(n[1:] * np.pi * bands) / (np.pi * n[1:]) ** 2
     b = np.dot(np.diff(b, axis=2)[:, :, 0], weight)
 
-    # Now we can solve the equation (use pinv because Q can be rank deficient)
-    a = np.dot(pinv(Q), b)
+    # Now we can solve the equation
+    try:  # try the fast way
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            a = solve(Q, b, sym_pos=True, check_finite=False)
+        for ww in w:
+            if ww.category == LinAlgWarning and \
+                    str(ww.message).startswith('Ill-conditioned matrix'):
+                raise LinAlgError(str(ww.message))
+    except LinAlgError:  # in case Q is rank deficient
+        # This is faster than pinvh, even though we don't explicitly use
+        # the symmetry here. gelsy was faster than gelsd and gelss in
+        # some non-exhaustive tests.
+        a = lstsq(Q, b, lapack_driver='gelsy')[0]
 
     # make coefficients symmetric (linear phase)
     coeffs = np.hstack((a[:0:-1], 2 * a[0], a[1:]))
diff --git a/scipy/signal/tests/test_fir_filter_design.py b/scipy/signal/tests/test_fir_filter_design.py
index 3714fe642b92..ff5a517d776f 100644
--- a/scipy/signal/tests/test_fir_filter_design.py
+++ b/scipy/signal/tests/test_fir_filter_design.py
@@ -5,9 +5,10 @@
                            assert_equal, assert_,
                            assert_allclose, assert_warns)
 from pytest import raises as assert_raises
+import pytest
 
+from scipy.linalg import LinAlgWarning
 from scipy.special import sinc
-
 from scipy.signal import kaiser_beta, kaiser_atten, kaiserord, \
         firwin, firwin2, freqz, remez, firls, minimum_phase
 
@@ -536,6 +537,27 @@ def test_compare(self):
         taps = firls(7, (0, 1, 2, 3, 4, 5), [1, 0, 0, 1, 1, 0], nyq=10)
         assert_allclose(taps, known_taps)
 
+        with pytest.raises(ValueError, match='between 0 and 1'):
+            firls(7, [0, 1], [0, 1], nyq=0.5)
+
+    def test_rank_deficient(self):
+        # solve() runs but warns (only sometimes, so here we don't use match)
+        x = firls(21, [0, 0.1, 0.9, 1], [1, 1, 0, 0])
+        w, h = freqz(x, fs=2.)
+        assert_allclose(np.abs(h[:2]), 1., atol=1e-5)
+        assert_allclose(np.abs(h[-2:]), 0., atol=1e-6)
+        # switch to pinvh (tolerances could be higher with longer
+        # filters, but using shorter ones is faster computationally and
+        # the idea is the same)
+        x = firls(101, [0, 0.01, 0.99, 1], [1, 1, 0, 0])
+        w, h = freqz(x, fs=2.)
+        mask = w < 0.01
+        assert mask.sum() > 3
+        assert_allclose(np.abs(h[mask]), 1., atol=1e-4)
+        mask = w > 0.99
+        assert mask.sum() > 3
+        assert_allclose(np.abs(h[mask]), 0., atol=1e-4)
+
 
 class TestMinimumPhase(object):
 


Code review comments:
diff_hunk: @@ -1032,8 +1035,20 @@ def firls(numtaps, bands, desired, weight=None, nyq=None, fs=None):
     b[1:] += m * np.cos(n[1:] * np.pi * bands) / (np.pi * n[1:]) ** 2
     b = np.dot(np.diff(b, axis=2)[:, :, 0], weight)
 
-    # Now we can solve the equation (use pinv because Q can be rank deficient)
-    a = np.dot(pinv(Q), b)
+    # Now we can solve the equation
+    try:  # try the fast way
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            a = solve(Q, b, sym_pos=True, check_finite=False)
+        for ww in w:
+            if ww.category == LinAlgWarning and \
path: scipy/signal/fir_filter_design.py
text: @user1:
I think we try to discourage the backslash line continuations in favor of using parentheses.

Less important is that there are no test cases where a different kind of warning is captured (partial line coverage), but I think we don't have to be that picky!