Repository: scipy/scipy
Pull request: ENH: add optional argument, "include_nyquist", for freqz()
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://docs.scipy.org/doc/numpy/dev/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
To fix #8265

#### What does this implement/fix?
add optional argument, "nyquist", for freq()


Problem statement:
freqz() does not include evaluation at pi by default
The default, or when specifying worN as an integer, the frequencies are evaluated at

```
linspace(0, 2*pi if whole else pi, N, endpoint=False)
```

However, I believe it would make more sense to *include* the endpoint, as the complex value of the frequency at pi is valuable information that I believe should be visible by default. This would parallel the often-taken step of taking a N-point FFT, and, for even N, keeping the values at indexes 0..N/2 *inclusive*. The last value represents the frequency at the Nyquist frequency.

Hints from issue discussion:
I agree in principle but changing this would be backward incompatible, and not worth introducing a deprecation cycle to change because the number of affected scripts and downstream libs would probably be quite large.

We could add a new keyword argument like `nyquist=False`, or perhaps simpler, allow `whole='nyquist'`, meaning "False but include Nyquist" or "only go to Nyquist". I might try to make a PR to do this to see how it would look.
Hello, I would like to work on this issue. Could you please point me in the right direction as to where I should start looking to fix this? Thank you!

Commit message: add optional argment, nyquist, for freqz

Patch to review:
diff --git a/scipy/signal/filter_design.py b/scipy/signal/filter_design.py
index 3b537f18fd33..b66535bc7d6f 100644
--- a/scipy/signal/filter_design.py
+++ b/scipy/signal/filter_design.py
@@ -272,7 +272,7 @@ def freqs_zpk(z, p, k, worN=200):
     return w, h
 
 
-def freqz(b, a=1, worN=512, whole=False, plot=None, fs=2*pi):
+def freqz(b, a=1, worN=512, whole=False, plot=None, fs=2*pi, nyquist=False):
     """
     Compute the frequency response of a digital filter.
 
@@ -301,7 +301,7 @@ def freqz(b, a=1, worN=512, whole=False, plot=None, fs=2*pi):
         If a single integer, then compute at that many frequencies (default is
         N=512). This is a convenient alternative to::
 
-            np.linspace(0, fs if whole else fs/2, N, endpoint=False)
+            np.linspace(0, fs if whole else fs/2, N, endpoint=nyquist)
 
         Using a number that is fast for FFT computations can result in
         faster computations (see Notes).
@@ -321,6 +321,8 @@ def freqz(b, a=1, worN=512, whole=False, plot=None, fs=2*pi):
         radians/sample (so w is from 0 to pi).
 
         .. versionadded:: 1.2.0
+    nyquist : bool, optional
+        If `nyquist` is True, `freqz` includes last frequency (Nyquist frequency) 
 
     Returns
     -------
@@ -436,7 +438,8 @@ def freqz(b, a=1, worN=512, whole=False, plot=None, fs=2*pi):
         if N < 0:
             raise ValueError('worN must be nonnegative, got %s' % (N,))
         lastpoint = 2 * pi if whole else pi
-        w = np.linspace(0, lastpoint, N, endpoint=False)
+        # if nyquist is true, w should include end point
+        w = np.linspace(0, lastpoint, N, endpoint=nyquist)
         if (a.size == 1 and N >= b.shape[0] and
                 sp_fft.next_fast_len(N) == N and
                 (b.ndim == 1 or (b.shape[-1] == 1))):
diff --git a/scipy/signal/tests/test_filter_design.py b/scipy/signal/tests/test_filter_design.py
index c7882496f8ed..13aef9c07ff8 100644
--- a/scipy/signal/tests/test_filter_design.py
+++ b/scipy/signal/tests/test_filter_design.py
@@ -810,6 +810,20 @@ def test_w_or_N_types(self):
             assert_array_almost_equal(w_out, [8])
             assert_array_almost_equal(h, [1])
 
+    def test_nyquist(self):
+        w, h = freqz([1.0], worN=8, nyquist=True)
+        assert_array_almost_equal(w, np.pi * np.arange(8) / 7.)
+        assert_array_almost_equal(h, np.ones(8))
+        w, h = freqz([1.0], worN=9, nyquist=True)
+        assert_array_almost_equal(w, np.pi * np.arange(9) / 8.)
+        assert_array_almost_equal(h, np.ones(9))
+
+        for a in [1, np.ones(2)]:
+            w, h = freqz(np.ones(2), a, worN=0, nyquist=True)
+            assert_equal(w.shape, (0,))
+            assert_equal(h.shape, (0,))
+            assert_equal(h.dtype, np.dtype('complex128'))
+
 
 class TestSOSFreqz(object):
 


Code review comments:
diff_hunk: @@ -436,7 +438,8 @@ def freqz(b, a=1, worN=512, whole=False, plot=None, fs=2*pi):
         if N < 0:
             raise ValueError('worN must be nonnegative, got %s' % (N,))
         lastpoint = 2 * pi if whole else pi
-        w = np.linspace(0, lastpoint, N, endpoint=False)
+        # if nyquist is true, w should include end point
+        w = np.linspace(0, lastpoint, N, endpoint=nyquist)
path: scipy/signal/filter_design.py
text: @user1:
As per the idea above:
```suggestion
        w = np.linspace(0, lastpoint, N, endpoint=nyquist and not whole)
```

@author:
@user1  I fixed the code based on your reviews.

diff_hunk: @@ -321,6 +321,8 @@ def freqz(b, a=1, worN=512, whole=False, plot=None, fs=2*pi):
         radians/sample (so w is from 0 to pi).
 
         .. versionadded:: 1.2.0
+    nyquist : bool, optional
+        If `nyquist` is True, `freqz` includes last frequency (Nyquist frequency) 
path: scipy/signal/filter_design.py
text: @user1:
Could maybe use a `.. versionadded:: 1.5.0` directive, but wait for Eric to review first.

@user2:
Yes it will indeed need this directive.

Also `nyquist` should probably only be used/considered in the `linspace` call when `not whole`. Otherwise in `whole` mode you'll include 0 and 2*pi, which is probably not desired. So the docs would then look something like:
```suggestion
        If `whole` is False and `worN` is an integer, setting `nyquist` to True
        will include the last frequency (Nyquist frequency).
```

@user3:
I would also reserve nyquist keyword to something else. `endpoint` or something less serious would be better.

@user2:
`include_nyquist` perhaps? Because it's really about whether or not to include the Nyquist freq. Agreed `nyquist` by itself is not great because it's not immediately clear you're not specifying the Nyquist rate or Nyquist freq of the data. 

@author:
Thank you for review. I fixed the code based on your reviews.