{{ message }}
ENH: signal.filtfilt supports complex-coefficient filters #22014
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Added unit test for padded (default) and Gustafsson's method. Fixes scipygh-17877.
|
I mistakenly thought scipy.signal doesn't support zpk->sos for complex-coefficient filters, though. The diff below tests sosfiltfilt using a custom, naive zpk->sos for the sosfiltfilt tests. That's good from the point-of-view that sosfiltfilt gets the right update, but the capability would be hard to use: users would need to compute their own second-order-sections. diff --git a/scipy/signal/_signaltools.py b/scipy/signal/_signaltools.py
index de4bbbb82..4cab95dca 100644
--- a/scipy/signal/_signaltools.py
+++ b/scipy/signal/_signaltools.py
@@ -4835,7 +4835,7 @@ def sosfiltfilt(sos, x, axis=-1, padtype='odd', padlen=None):
x_0 = axis_slice(ext, stop=1, axis=axis)
(y, zf) = sosfilt(sos, ext, axis=axis, zi=zi * x_0)
y_0 = axis_slice(y, start=-1, axis=axis)
- (y, zf) = sosfilt(sos, axis_reverse(y, axis=axis), axis=axis, zi=zi * y_0)
+ (y, zf) = sosfilt(sos.conj(), axis_reverse(y, axis=axis), axis=axis, zi=zi.conj() * y_0)
y = axis_reverse(y, axis=axis)
if edge > 0:
y = axis_slice(y, start=edge, stop=-edge, axis=axis)
diff --git a/scipy/signal/tests/test_signaltools.py b/scipy/signal/tests/test_signaltools.py
index 8013d824d..80aae9131 100644
--- a/scipy/signal/tests/test_signaltools.py
+++ b/scipy/signal/tests/test_signaltools.py
@@ -2307,6 +2307,27 @@ class TestLFilterZI:
assert_equal(np.real(signal.lfilter_zi(b, a)).dtype, dtype)
+def pairwise(flat):
+ """Make a list pairwise, with 0-tail padding"""
+ pairs = list(zip(flat[::2],flat[1::2]))
+ if len(flat) % 2 == 1:
+ pairs.append((flat[-1],0))
+ return pairs
+
+
+def zpk2sos_cplx(z, p, k):
+ """Naive zpk -> sos transformation for complex-coefficient filters"""
+ pairz = pairwise(z)
+ pairp = pairwise(p)
+ order = max(len(pairz), len(pairp))
+ pairz.extend([(0,0)] * (order - len(pairz)))
+ pairp.extend([(0,0)] * (order - len(pairp)))
+ soslist = []
+ for z, p in zip(pairz, pairp, strict=True):
+ soslist.append(np.hstack(signal.zpk2tf(z, p, 1)))
+ soslist[0][:3] *= k
+ return np.array(soslist)
+
class TestFiltFilt:
filtfilt_kind = 'tf'
@@ -2316,7 +2337,10 @@ class TestFiltFilt:
b, a = zpk2tf(*zpk)
return filtfilt(b, a, x, axis, padtype, padlen, method, irlen)
elif self.filtfilt_kind == 'sos':
- sos = zpk2sos(*zpk)
+ try:
+ sos = zpk2sos(*zpk)
+ except ValueError:
+ sos = zpk2sos_cplx(*zpk)
return sosfiltfilt(sos, x, axis, padtype, padlen)
def test_basic(self):
@@ -2405,8 +2429,6 @@ class TestFiltFilt:
def test_complex_coeff(self):
# gh-17877: complex coefficient filters
- if self.filtfilt_kind != 'tf':
- pytest.skip('complex coefficents only for tf')
# centre frequency for filter [/time]
fcentre = 50 |
Minor conflict in test_signaltools.py, kept both sides.
|
The merge conflict itself wasn't complicated. I updated the new test to use |
|
the failures are in I can easily enough change back to |
ev-br
reviewed
Dec 14, 2024
There was a problem hiding this comment.
We probably want an opinion of @larsoner and/or @ilayn since complex-valued story is scipy.signal is a bit tricky.
I personally have no opinion and would be happy to go with what you three agree on :-).
Most of my comments below are purely janitorial and assume you wanted to make the new test you added array api compatible. If not, it's fine to decorate the test with @skip_xp_backends(np_only=True) + a code comment to remind us that it's just lack of time not an inherent reason it cannot run on alternative array backends.
Also note that this part of scipy.signal is in a bit of a flux w.r.t. array api compatibility, and domain expert help in converting it would be great, actually. This is a bit baffling at first but is no rocket science and we will help :-).
filtfilt code itself is a bit weird: it's all numpy + xp.asarray at the very end because lfilter is numpy-only compiled code. So there's little reason to generalize filtfilt internals, they can keep being what they are.
| @@ -4331,8 +4331,16 @@ def _filtfilt_gust(b, a, x, axis=-1, irlen=None): | |||
| ---------- | |||
| b : scalar or 1-D ndarray | |||
| Numerator coefficients of the filter. | |||
|
|
|||
| ..versionchanged:: 1.15.0 | |||
| Complex-valued `b` is correctly handled. | |||
There was a problem hiding this comment.
TBH, "is correctly handled" does not tell a reader much. Maybe this needs a longer explanation in the Notes section which you refer the reader from here.
Also, is this change backwards compatibile? If not, is it a bug fix or a behavior change --- in the latter case you migh need to argue about a deprecation process.
There was a problem hiding this comment.
TBH, "is correctly handled" does not tell a reader much. Maybe this needs a longer explanation in the Notes section which you refer the reader from here.
Thanks, that is a good idea.
Also, is this change backwards compatibile? If not, is it a bug fix or a behavior change --- in the latter case you migh need to argue about a deprecation process
A good point.
It's backwardly compatible for real-coefficient filters, but not for complex-coefficient filters. I can't believe anyone is relying on the current behaviour for complex-coefficient filters, but it's remotely conceivable. The effect of the unchanged code is to apply the filter and its conjugate, but only once each.
There was a problem hiding this comment.
Sounds like it's clear a bugfix then, right?
There was a problem hiding this comment.
Sounds like it's clear a bugfix then, right?
Yes. Shall I change the PR title?
Here's what I've done for the "version changed" wording:
--- a/scipy/signal/_signaltools.py
+++ b/scipy/signal/_signaltools.py
@@ -4523,15 +4523,15 @@ def filtfilt(b, a, x, axis=-1, padtype='odd', padlen=None, method='pad',
b : (N,) array_like
The numerator coefficient vector of the filter.
- .. versionchanged:: 1.15.0
- Complex `b` is correctly handled.
+ .. versionchanged:: 1.16.0
+ Change in behaviour for complex `b`, see ``Notes``.
a : (N,) array_like
The denominator coefficient vector of the filter. If ``a[0]``
is not 1, then both `a` and `b` are normalized by ``a[0]``.
- .. versionchanged:: 1.15.0
- Complex `a` is correctly handled.
+ .. versionchanged:: 1.16.0
+ Change in behaviour for complex `a`, see ``Notes``.
x : array_like
The array of data to be filtered.
@@ -4586,6 +4586,13 @@ def filtfilt(b, a, x, axis=-1, padtype='odd', padlen=None, method='pad',
The option to use Gustaffson's method was added in scipy version 0.16.0.
+ .. versionchanged:: 1.16.0
+
+ `b` and `a` are complex-conjugated when applying the backward
+ filter. This conjugation reflects the filter's frequency
+ spectrum through zero, corresponding to the effect on the
+ signal's frequency spectrum when it is reversed in time.
+
References
----------
.. [1] F. Gustaffson, "Determining the initial states in forward-backwardThere was a problem hiding this comment.
Looks reasonable to me!
One last test idea, would it be worth testing some simple complex filter (2-tap IIR?) on a short signal (5 samples?) with MATLAB or Octave with simple edge padding (constant?) and hard-coding the result in tests to check against?
There was a problem hiding this comment.
One last test idea, would it be worth testing some simple complex filter (2-tap IIR?) on a short signal (5 samples?) with MATLAB or Octave with simple edge padding (constant?) and hard-coding the result in tests to check against?
Good idea. I'll see if Octave supports complex-coefficient filters.
| # sample rate [/time] | ||
| fs = 1e3 | ||
|
|
||
| z, p, k = signal.butter(2, 2*xp.pi*fwidth/2, output='zpk', fs=fs) |
There was a problem hiding this comment.
as long as butter is numpy-only , you need to wrap z and p into xp.asarray(). Please also add a comment # XXX while butter is numpy-only so that we can search-and-remove when butter is converted.
To expand that comment:
TL;DR:
EDIT:
Locally, use the If you want to experiment in ipython, add the env variable: |
|
I lost track of this so if you folks are happy with it, I'm OK with it, @DietBru wanna take a look? |
TL;DR: Since you defined In general, numpy is an odd one out which allows python scalars in and returns arrays. The next revision of the spec will most likely allow python PODs in binary functions (so that you will be able to |
|
Thanks @ev-br for all the help.
This was a good tip. We still need to resolve the backward-compatibility issue, and I will meanwhile reword the explanation of the change. |
|
All in all, this PR looks quite plausible to me. I have to admit that I do not fully understand Gustaffson's method and its implementation. Toying with a simple filter gives: b, a = np.array([1j]), np.array([1, -0.5])
x = np.array([0, 0, 0, 0, 1, 0, 0, 0, 0], dtype=complex)
fig1, (ax10, ax11) = plt.subplots(2, 1, sharex='all', tight_layout=True)
ax10.set(title=f"Response for {b=}, {a=}", ylabel="Real")
ax11.set(xlabel="Samples", ylabel="Imag")
for n_ in ('pad', 'gust'):
y = filtfilt(b, a, x, method=n_)
ax10.step(y.real, '.-', where='post', alpha=0.5, label=n_)
ax11.step(y.imag, '.-', where='post', alpha=0.5, label=n_)
for ax_ in (ax10, ax11):
ax_.grid()
ax_.legend()
plt.show()
|
Interesting test case, thanks; I'll take a closer look. |
|
@larsoner wrote:
I don't have access to Matlab. I read their documentation for Both Octave's Although I'm fairly sure my proposed change is correct, this overall lack of support could point to some convention I'm ignorant of, though it's also plausible this is just too narrow a use case to come up often. Octave: Julia: |
|
Summary: I think the Gustafsson result is correct. @DietBru wrote:
Gustafsson's method finds filter initial conditions such that one gets the same output whether the forward filter is run first, and then the reverse, or the reverse and then the forward. As I understand, these initial conditions are unique, assuming filter input longer than filter order (or thereabouts: it may be filter order + 1, or - 1). From the script below this forward-reverse condition is satisfied by the PR code, but the result is not the same as the # gh-22014 discussion re gustafsson vs pad result
import numpy as np
from scipy import signal
from scipy.signal._signaltools import _filtfilt_gust
b, a = np.array([-1j]), np.array([1, -0.5])
x = np.array([0, 0, 0, 0, 1, 0, 0, 0, 0])
ygust, x0, x1 = _filtfilt_gust(b, a, x)
# forward then reverse
yfr = signal.lfilter(b.conj(), a.conj(),
signal.lfilter(b, a, x, zi=x0)[0][::-1], zi=x1)[0][::-1]
# compare against filtfilt output; computed in slightly different way
print('y-forward-reverse vs ygust')
print(f' {max(abs(yfr-ygust))=}')
# reverse then forward
yrf = signal.lfilter(b, a,
signal.lfilter(b.conj(), a.conj(),
x[::-1], zi=x1)[0][::-1], zi=x0)[0]
# compare forward and reverse
print('y-forward-reverse vs y-reverse-forward')
print(f' {max(abs(yrf-yfr))=}')
print()
print('symmetry of ygust')
print(f' {ygust[:4].real - ygust[5:][::-1].real=}')
print(f' {ygust[:4].imag + ygust[5:][::-1].imag=}')
# symmetry of ypad
ypad = signal.filtfilt(b, a, x)
print()
print('symmetry of ypad')
print(f' {ypad[:4].real - ypad[5:][::-1].real=}')output |
Torch doesn't support negative strides. pytorch/pytorch#59786
Looking at your little script, I was able to reverse the symmetry by make the input b, a = np.array([-1j]), np.array([1, -0.5])
x = np.array([0, 0, 0, 0, 1j, 0, 0, 0, 0])It smells to me like a complex conjugation is missing or misplaced somewhere. I do not know if the problem is in Gustavsson's paper or in the implementation. I am sorry that currently I do not have the bandwidth to dive deeper into Gustavsson's method to write you something more meaningful. |
|
Thanks for your patience -- this is a obscure corner of a technical topic.
(For anyone else reading, real part of output is now odd & imaginary part is now even.) I think Why I think the current method is right (this is not a watertight argument): If I'm right, there's a is counter-intuitive consequence: the
Results, using this branch: |
|
I finally found the time to study Gustafsson paper a bit, but there are still quite a few question marks for me: I do not understand how to integrate the conjugation into the following properties: E.g., the second line Hence, how does the conjugation of Also Matlab's filtfilt does not produce the same result: >> b = 1i; a = [-1, 0.5];
>> x = [0, 0, 0, 0, 1, 0, 0, 0, 0];
>> y = filtfilt(b, a, x)
y =
-0.0833 -0.1667 -0.3333 -0.6667 -1.3334 -0.6667 -0.3335 -0.1670 -0.0840
PS: You can use Matlab online for 20h per month for free, if you register. |
|
There seems to be an ongoing discussion between @roryyorke and @DietBru , and it's unlikely the discussion resolves before branching 1.16. Am bumping the milestone. Would be great to arrive at a conclusion, ideally early in the 1.17 cycle, so that a potentially good change does not keep lingering. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Reference issue
Closes gh-17877.
What does this implement/fix?
Correctly handle complex-coefficient filters in
filtfiltfor both default and Gustafsson's method.Additional information
Arguably a fix, not an enhancement.
I'm not sure about the
.. versionchangeddirectives I added:Notesor similar?I haven't addressed the need for possibly different symmetry in imaginary part, nor making the Gustafsson's unit test more efficient (see gh-17877 for remarks on both of these).