[Preprocessing] add causal backward filtering to correct hardware induced phase shift

[JuanPimientoCaicedo]

Hi, guys. I worked in a new class to include a causal backward filtering to the Spikeinterface preprocessing pipeline.

I did this to be able to correct for the hardware induced phase shift in the lf and ap band during neuropixel recordings. This is the first time I do this, so I don't know if this is the correct way of proposing these changes.

Please let me know what you think!

[JuanPimientoCaicedo]

PR related to issue #2943.

[alejoe91]

maybe it would be better to have causal_mode:

• filtfilt
• forward
• backward
  ?

[alejoe91]

Please revert these changes. There should be a space before the :

[JoeZiminski]

Hey @JuanPimientoCaicedo I think the git history went funny, this can happen if you are unlucky and a big formatting change is merged after you branched off from SpikeInterface main. One suggestion, as no one else as contributed to this branch, is rebasing off an updated main branch and then force pushing. To do this you would:

1. make a backup of your local repository that is synced to the remote repository. This is useful as force-pushing will overwrite the remote respository (this can just be copying the local repo after ensuring your are in sync with remote using git pull).
2. on your fork's main branch, use the 'Sync' button to sync to SpikeInterface main branch.
3. checkout main branch on your local fork and pull the most recent main branch
4. checkout your filter_update branch and do (if using the terminal) git rebase -i main. You will then see a screen like the second picture down in this nice guide that will let you keep / drop commits. I believe if you drop commits after and including 030c534 by changing the pick to drop you will keep all of your changes and drop all of the merging / reverting changes that came after.
5. Perform rebase and carefully check you didn't loose any work
6. do git push -f to force-push to your branch. If you realise something is missing you can go onto your backup repo and force-push from there to revert to what you have now.

There are other ways around this problem you may want to try and do not involve force pushing / rebasing, but this is the one that comes first to my mind. Please let me know if you want to try this approach and have any questions!

[JuanPimientoCaicedo]

I will sync the fork and make the changes again, there were just a few lines, thanks @JoeZiminski.

[h-mayorquin]

Hey, I find the semantics of this a bit confusing.

Do I understand correctly that you goal is to rever the phase shift introduced by harward causal filters? If that's correct I have some suggestions that I think can improve the semantics and the maintenance.

You are trying to rever the effects of a causal filter but the process that is used for that is not causal, hence, I don't think CausalFilterRecording is a good name or causal is a good keyword to route this behavior in the corresponding recording segment. I would call it by what you are doing which is BackwardsFilter and then explain in the docs that if you run this with the same parameters of the Forward Filter of the hardward you can get rid of the phase shift.

I have some more comments but I have to go. I will add more later.

[h-mayorquin]

This and the same comment for line 154 is unecessarily fliping the channels as well. I think you should only flip the time axis.

[JuanPimientoCaicedo]

Hi, yes. You are right about the channels flipping as well. that is missing: the correct line should look like this, right?:

filtered_traces = np.flip(scipy.signal.lfilter(b, a, np.flip(traces_chunk, axis = 0), axis=0))

Yes, BackwardsFilter sounds better and it is more self explanatory.

[h-mayorquin]

Yes, I am wondering if just using the reverse notation traces[::-1, :] is better. I think that using this instead of flip creates a view instead of a new memory allocation.

The scipy filter functions are really memory hungry. I will test this.

[JoeZiminski]

For me, this should be generalised to a causal filter implementation rather than its implementation tailored to a specific goal. Using a causal filter to reverse the phase shift from the NP hardware would be just a specific use-case of the more generally implemented causal filter e.g causal_filter(recording, direction="backwards")). There are reasons that people might want to use a causal rather than non-causal IIR filter implementation (e.g. to replicate other processing pipelines e.g. CatGT biquad filters) so I think it makes sense to include as a generic causal filter. To me (imaging I don't know about the reasoning to reverse the head-stage filtering) it is a little confusing that the only causal filter implementation is to backwards filter the signal.

BTW @JuanPimientoCaicedo apologies for the delay reply to the message in the thread, I think in general what you wrote is correct but note Jennifer comment on the Neuropixels slack thread, the actual filter implementation in the headstage may be more complex than first order RC filters and may not be well modelled by first order butterworth. ATM I am looking a little more into these switched capacitor filters, the transfer function of some typical implementations can be found here (slide 37) so was interested in looking at how different their phase response is from the standard RC filter and responding properly ASAP. But like Jennifer said I think the safest way is to get the exact implementation details from imec, maybe we can meet with her and anyone else interested to help get a response from imec.

[JuanPimientoCaicedo]

Hi, @JoeZiminski. No worries!

You are right, I had the idea from a previous slack comment that they were single pole RC filters. It looks like instead of assuming some configuration, it is important to use the specific parameters as @oliche pointed out. Looking forward to having that information at hand.

Jennifer said that she was going to ask the imec engineers almost 2 weeks ago, however, she has not posted anything new since then... I wonder if they replied back.

[h-mayorquin]

@JoeZiminski I am thinking

For me, this should be generalised to a causal filter implementation rather than its implementation tailored to a specific goal. Using a causal filter to reverse the phase shift from the NP hardware would be just a specific use-case of the more generally implemented causal filter e.g causal_filter(recording, direction="backwards")). There are reasons that people might want to use a causal rather than non-causal IIR filter implementation (e.g. to replicate other processing pipelines e.g. CatGT biquad filters) so I think it makes sense to include as a generic causal filter. To me (imaging I don't know about the reasoning to reverse the head-stage filtering) it is a little confusing that the only causal filter implementation is to backwards filter the signal.

This makes me reconsider that backwards might be a bad name but I am afraid of getting too much into details that will bog down the contribution of this PR.

Right now we have two kind of preprocessing filters in the API:

• Purposes specific filters like bandpass and notch
• Infraestructure / Machinery filters like FilterRecording.

Note that FilterRecording is our NonCausalFilter as it is implemented with scipy.filtfilt or the sos counterpart.

I feel that at the implementation, testing and interface of this would be easier if we think on this contribution as a purpose specific filter: That is, we are adding TheFilterThatAllowsYouToRevertThePhaseEffectsOfCasualFilters and we can focus on making this sensible.

I follow the heuristic that in most cases it is better to create an example and then generalize after but maybe I am wrong here.

I really think that this should be the call of @JuanPimientoCaicedo on what he wants to do because both scopes seem very different. He can either want to create a specific purpose filter with interface tailored for this case (the internals we can adjust later) or maybe he wants to create general machinery for 1) causal filters, 2) applying filters with reverse time, than then the user can use for this purpose of something else. Note that the documentation and testing of both cases would be very different with the latter being harder I feel.

[JoeZiminski]

In general I would agree with you that first making a single example is better and then generalise later to avoid premature generalisation. However in this case I think it is worth making it general for the following reasons:

1. It is very close to being general, the only change is to expose the direction and it becomes a general-purpose causal filter. I don't think this will add much overhead in terms of testing or documentation.
2. It breaks with the existing convention that all filters are simple exposure of the underlying scipy function. Although we have filter recording / bandpass / notch, these are all still simple wrappers around the scipy functions. Instead of now wrapping the sosfilt we are making additional restrictions / assumptions on the applied use case.
3. It is inevitable that someone is going to want to perform forward causal filtering in SI at some point. Then we will either have to change this implementation anyway or create a duplicate forward version.

[JuanPimientoCaicedo]

I like the generalized causal filter implementation with the causal_filter(recording, direction="forward/backward")). I do understand that backwards filtering is just one of the possible uses of the function and as long as you are able to define the direction of the filter, I am happy with it.

Now, @h-mayorquin points out something important, It might be better to develop that function outside of the current FilterRecording functionality. I believe that is your call.

Jennifer pointed out she does not have the filter details yet. I am wondering if it is possible to obtain the filter responses in any other way, in order to try to advance on the specific details of this implementation.

[h-mayorquin]

I like the generalized causal filter implementation with the causal_filter(recording, direction="forward/backward"))

@JoeZiminski @JuanPimientoCaicedo this makes sense to me. I just did not want to burden Juan with something more general.

Creating a new CausalFilter class based on lfilt and sosfilt instead of (filtfilt and sosfiltfilt) seems like a good addition to the library to me.

Jennifer pointed out she does not have the filter details yet. I am wondering if it is possible to obtain the filter responses in any other way, in order to try to advance on the specific details of this implementation.

How do you think the details of the implementation will affect the PR here? If we are going for the more general version, the API should expose generic parameters that define the filter like this:
https://github.com/h-mayorquin/spikeinterface/blob/199bf60d13f144a70e29fee126eb946d52cfc2c3/src/spikeinterface/preprocessing/filter.py#L35-L51

And then, once you have those details you can create a filter, pass the parameters as in the above with direction=backward and then the object should perform backward causal filter with the filters. Am I making sense?

[JuanPimientoCaicedo]

Sorry for the delay, I am going to work on this later this week. Yeah, it makes sense. So, to be clear, I am planning to create a new class called CausalFilter that works independently of the FilterRecording class. The main differences there are the implementation of sosfilt and lfilt, instead of the zero-phase versions and the forward/backward functionality. Is that what you are expecting?

[h-mayorquin]

That's what I am proposing to avoid having way too many branches on the segment of the current filter.

[JoeZiminski]

I like the idea of not overloading FilterRecording but I'm not sure how easy it would be to develop CausalFilterRecording independently without duplicating a lot of the code. At present I don't think the addition of causal makes the get_traces() too heavy, I would more fear errors creeping in from maintaining duplicate implementations of the filtering machinery, but maybe there is an efficient way of doing this I've overlooked 🤔

True though as you mentioned though once we find the information about the headstage filters we might need to expose more parameters anyway. Maybe the easiest way will be to take directly b, a or sos coefficients generated from the transfer function of the switched capacitor filters. But I wonder at this stage if it is easiest to implement the causal aspect here and maybe have a bigger refactor where we could address things like increasing complexity of FilterRecording?

[h-mayorquin]

Yeah, we come at this from very different places. We have had this disagreement before. Outside of code by PhD students or data scientists in my previous job I never seen a place that really needs more DRY. Most coders love abstractions and the concept of generalizing. On the other hand, I have seen too many early generalizations, conditional structures with number_of_branches > 5, wrong abstractions, gigantic one-liners, etc. I would be curios what do you think of the common criticisms to DRY as a principle but let's have this discussion outside of this thread.

My concrete argument here:
The proposed changes introduced two more branches to the segment: cause and direction . The direction parameter does not really make sense in the context of FilterRecording so you would need to handle that interaction both in the code and in the docstring. What do you gain by paying this cost? Avoid having to code managing the margins which is done most preprocessors anyway.

[h-mayorquin]

Anyway, I think that the call should be of @JuanPimientoCaicedo . This is just what I believe is best based on the reasons above but I think both implementations work.

[JoeZiminski]

Yes I think it is a philosophical difference in approaches, it will be good to discuss more some time 😄 For me the complex structures that can occur when trying to apply DRY are a problem but my preferred solution is to better modularise them but not at the expense of DRY. One area that I found DRY key was making it easy to add new features to codebases without introducing bugs. Anyways, a discussion for another day!

I agree @JuanPimientoCaicedo please go ahead in whatever way feels natural. I was thinking it would require complete rewriting of FilterRecording. But from your comment I see did you mean to subclass FilterRecording and only override get_traces(). I didn't think of that option! Agree it is weird to pass "backwards"/"forwards" to FilterRecording.

[JuanPimientoCaicedo]

Hi, guys. I just started working in the new class and noticed that the get_traces() method is from the FilterRecordingSegment class. And now I am wondering which could be the best way to add this feature without affecting the current functionality. should I create two subclasses like this?:

class CausalFilterRecording(FilterRecording):
    """
    """
class ForwardBackwardFilter(FilterRecordingSegment):
    """
    """
    def get_traces(self, start_frame, end_frame, channel_indices):

Is there a better way to do it?

[h-mayorquin]

Yes, if you create a new class you will need to create both a recording and a recording segment. I don't think you need to inherit from the other classess directly. You can inherit from the abstract class themselves.

[samuelgarcia]

We are reading this implementation and discuss with lots of delays.
Sorry guys.
Alessio and I have the feeling that injection more option (direction) to FilterRecording would more efficient that making a new class which is very simlar for only one parameters more.
Then the direction would "forward-backward" (default) | "forward" | "backward"

We could propose to modify the FilterRecording instead but we could have a function 'causal_filter` that use this class and explicitly expose the direction parameter.
What do you think ?

[h-mayorquin]

That's what @JoeZiminski was proposing, @samuelgarcia .

I don't like it because the filtfilt family functions are forward-backward only and I think you will end up with a more complicated functions with more if branches.

The final proposal though was that @JuanPimientoCaicedo implements it in the way that he finds it easier and we can refactor aftewards.

[JuanPimientoCaicedo]

Hi, @samuelgarcia.

I can create another PR with that implementation and that way you can compare both approaches. You can choose the one you prefer. How does that sound?

[alejoe91]

Closing in favor of #3133