Back to the main page.
Bug 780 - filter instability
| Status | CLOSED FIXED |
| Reported | 2011-06-28 14:24:00 +0200 |
| Modified | 2013-04-23 12:06:31 +0200 |
| Product: | FieldTrip |
| Component: | preproc |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Mac OS |
| Importance: | P1 enhancement |
| Assigned to: | Robert Oostenveld |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
Robert Oostenveld - 2011-06-28 14:24:54 +0200
On 24 Jun 2011, at 16:04, litvak.vladimir@gmail.com wrote: Hi Robert, I've been getting several complaints recently about filter instability that we discussed before (twice in the past week). So it is not so uncommon after all. You mentioned at the time that this could be circumvented by combining two filters. Would it be possible to implement that to put this to rest once and for all? Thanks, Vladimir
Robert Oostenveld - 2011-06-28 14:25:34 +0200
On 26 Jun 2011, at 10:53, Vladimir Litvak wrote: Hi Robert, I'm not at HBM. Maybe next year. Yes, there is an error now, but the problem is that it happens quite often and not for some weird settings. See: http://www.jiscmail.ac.uk/cgi-bin/wa.exe?A2=SPM;c8cd3970.1106 http://www.jiscmail.ac.uk/cgi-bin/wa.exe?A2=SPM;3b25fefb.1106 I remember that when we first discussed it you said that it's possible to avoid the problem by combining two filters with odd and even order. So what I was asking is if that'd be possible to implement that solution to avoid people complaining. Vladimir ============================= On Sun, Jun 26, 2011 at 3:45 PM, Robert Oostenveld <r.oostenveld@donders.ru.nl> wrote: Hi Vladimir Stefan was the expert on filtering and he made preproc/private/filter_with_correction. That is used in mbp> grep with_corr *.m ft_preproc_bandpassfilter.m:filt = filter_with_correction(B,A,dat,dir); ft_preproc_bandstopfilter.m:filt = filter_with_correction(B,A,dat,dir); ft_preproc_highpassfilter.m:filt = filter_with_correction(B,A,dat,dir); ft_preproc_lowpassfilter.m:filt = filter_with_correction(B,A,dat,dir); The filter_with_correction should give an error if the filter is instable. Could you provide a clear case of instability that is not caught by this low-level function?
Robert Oostenveld - 2011-06-28 14:28:51 +0200
In the first SPM mail it is not clear to me whether the filter function throws an explicit error, or whether it returns NaNs. It would be good to stop with the debugget at the begin of preproc_bandpassfilter and to know the exact inputs (includiong the data) to that function. The second mail is more explicit about the details and there the instability is correctly "caught" as error, including instructions on how to resolve it. That one allowed me to reproduce it as follows: >> dat = randn(1,1000); >> ft_preproc_highpassfilter(dat, 1000, 0.1) poles = 1.0032 + 0.0021i 1.0032 - 0.0021i 0.9996 + 0.0042i 0.9996 - 0.0042i 0.9960 + 0.0021i 0.9960 - 0.0021i ??? Error using ==> filter_with_correction at 42 Calculated filter coefficients have poles on or outside the unit circle and will not be stable. Try a higher cutoff frequency or a different type/order of filter. Error in ==> ft_preproc_highpassfilter at 79 filt = filter_with_correction(B,A,dat,dir); In the code I have now found the piece where I had implemented the repeated application of the filters, see ft_preproc_bandstopfilter. It also contains the comment from Stefan Klanke. The problem (or at least one problem that I understand) with the repeated lower-order filtering is that the filter order in the result is not correct. You specify the frequency at which the 3 dB point should be, i.e. the corner point in the transfer function. If you do two filters in succession, the corner point is not 3 dB down (i.e. the spectral content is attenuated with a factor of 2x), but twice that amount. The consequence is also that the exact specification of the passband is not controlled any more, i.e. the 3dB point shifts into the passband. I feel that Stefans point is correct, i.e. the repeated application does not solve the problem in the desired way. Perhaps better would be to automatically reduce the filter order, i.e. like >> ft_preproc_highpassfilter(dat, 10000, 0.01, 5) >> ft_preproc_highpassfilter(dat, 10000, 0.01, 4) >> ft_preproc_highpassfilter(dat, 10000, 0.01, 3) the 5 and 4 fail, but the 3 works. If you agree with the proposed solution, this would have to be implemented with a try-catch around the filter_with_correction in each of the lp/hp/bp/bs functions.
Vladimir Litvak - 2011-06-28 14:36:25 +0200
It's likely that the first mail pertains to older SPM version before the instability detection was implemented. What are the consequences of reducing filter order? Does it also change the attenuation? If the results are not dramatically affected then perhaps we should do it automatically. Or perhaps you should make it an option that at least in SPM will be enabled by default but at least the user will know that filter order can change. Vladimir (In reply to comment #2) > In the first SPM mail it is not clear to me whether the filter function throws > an explicit error, or whether it returns NaNs. It would be good to stop with > the debugget at the begin of preproc_bandpassfilter and to know the exact > inputs (includiong the data) to that function. > > The second mail is more explicit about the details and there the instability is > correctly "caught" as error, including instructions on how to resolve it. That > one allowed me to reproduce it as follows: > > >> dat = randn(1,1000); > >> ft_preproc_highpassfilter(dat, 1000, 0.1) > > poles = > > 1.0032 + 0.0021i > 1.0032 - 0.0021i > 0.9996 + 0.0042i > 0.9996 - 0.0042i > 0.9960 + 0.0021i > 0.9960 - 0.0021i > > ??? Error using ==> filter_with_correction at 42 > Calculated filter coefficients have poles on or outside the unit circle and > will not be stable. Try a higher cutoff frequency or > a different type/order of filter. > > Error in ==> ft_preproc_highpassfilter at 79 > filt = filter_with_correction(B,A,dat,dir); > > In the code I have now found the piece where I had implemented the repeated > application of the filters, see ft_preproc_bandstopfilter. It also contains the > comment from Stefan Klanke. > > The problem (or at least one problem that I understand) with the repeated > lower-order filtering is that the filter order in the result is not correct. > You specify the frequency at which the 3 dB point should be, i.e. the corner > point in the transfer function. If you do two filters in succession, the corner > point is not 3 dB down (i.e. the spectral content is attenuated with a factor > of 2x), but twice that amount. The consequence is also that the exact > specification of the passband is not controlled any more, i.e. the 3dB point > shifts into the passband. > > I feel that Stefans point is correct, i.e. the repeated application does not > solve the problem in the desired way. Perhaps better would be to automatically > reduce the filter order, i.e. like > > >> ft_preproc_highpassfilter(dat, 10000, 0.01, 5) > >> ft_preproc_highpassfilter(dat, 10000, 0.01, 4) > >> ft_preproc_highpassfilter(dat, 10000, 0.01, 3) > > the 5 and 4 fail, but the 3 works. > > If you agree with the proposed solution, this would have to be implemented with > a try-catch around the filter_with_correction in each of the lp/hp/bp/bs > functions.
Robert Oostenveld - 2011-06-28 16:06:54 +0200
YeP, reducing filter order reduces the attenuation. Designing a filter with a very steep cutoff is non trivial. If needed, someone should spend some time on extending the options for ft_preproc_xxxfilter to add another type of filter which is able to have a steeper cutoff and stronger attenuation. But I don't see the scientific reason for trying to make the attenuation so large. How should the interface look like for making this optional? Should we change the functions to key-value style, or just have a few keyvals at the end of the existing list, or just add an extra input argument like this function [filt] = ft_preproc_bandpassfilter(dat, Fs, Fbp, N, type, dir, instabilityfix) if nargin<7 instabilityfix = 'reduce'; % can be 'reduce' or 'split' end ... try filt = filter_with_correction(B,A,dat,dir); catch ME switch instabilityfix case 'reduce' warning('...) filt = ft_preproc_bandpassfilter(dat, Fs, Fbp, N-1, type, dir, instabilityfix); case 'split' % here comes the repeated filtering that is already in bandstopfilter otherwise retrow(ME) end % case end % try-catch
Vladimir Litvak - 2011-06-28 16:11:49 +0200
(In reply to comment #4) > YeP, reducing filter order reduces the attenuation. Designing a filter with a > very steep cutoff is non trivial. If needed, someone should spend some time on > extending the options for ft_preproc_xxxfilter to add another type of filter > which is able to have a steeper cutoff and stronger attenuation. But I don't > see the scientific reason for trying to make the attenuation so large. > There is one reason I know about. I talked to Burkhard Maess at last Biomag and he mentioned that when he does ERFs he always designs high-pass filters by hand to have very high attenuation because otherwise the baseline is not flat enough. So when you want to use filtering instead of baseline correction you need to make sure that low frequencies and DC are very strongly attenuated. > How should the interface look like for making this optional? Should we change > the functions to key-value style, or just have a few keyvals at the end of the > existing list, or just add an extra input argument like this > I don't mind. Whatever you think is most consistent with other interfaces in FT. Vladimir > function [filt] = ft_preproc_bandpassfilter(dat, Fs, Fbp, N, type, dir, > instabilityfix) > > if nargin<7 > instabilityfix = 'reduce'; % can be 'reduce' or 'split' > end > > ... > > try > filt = filter_with_correction(B,A,dat,dir); > catch ME > switch instabilityfix > case 'reduce' > warning('...) > filt = ft_preproc_bandpassfilter(dat, Fs, Fbp, N-1, type, dir, instabilityfix); > case 'split' > % here comes the repeated filtering that is already in bandstopfilter > otherwise > retrow(ME) > end % case > end % try-catch
Vladimir Litvak - 2011-06-28 16:15:39 +0200
Maybe also add a third option 'error' with the present behavior. Vladimir (In reply to comment #5) > (In reply to comment #4) > > YeP, reducing filter order reduces the attenuation. Designing a filter with a > > very steep cutoff is non trivial. If needed, someone should spend some time on > > extending the options for ft_preproc_xxxfilter to add another type of filter > > which is able to have a steeper cutoff and stronger attenuation. But I don't > > see the scientific reason for trying to make the attenuation so large. > > > > There is one reason I know about. I talked to Burkhard Maess at last Biomag and > he mentioned that when he does ERFs he always designs high-pass filters by hand > to have very high attenuation because otherwise the baseline is not flat > enough. So when you want to use filtering instead of baseline correction you > need to make sure that low frequencies and DC are very strongly attenuated. > > > How should the interface look like for making this optional? Should we change > > the functions to key-value style, or just have a few keyvals at the end of the > > existing list, or just add an extra input argument like this > > > > I don't mind. Whatever you think is most consistent with other interfaces in > FT. > > Vladimir > > > function [filt] = ft_preproc_bandpassfilter(dat, Fs, Fbp, N, type, dir, > > instabilityfix) > > > > if nargin<7 > > instabilityfix = 'reduce'; % can be 'reduce' or 'split' > > end > > > > ... > > > > try > > filt = filter_with_correction(B,A,dat,dir); > > catch ME > > switch instabilityfix > > case 'reduce' > > warning('...) > > filt = ft_preproc_bandpassfilter(dat, Fs, Fbp, N-1, type, dir, instabilityfix); > > case 'split' > > % here comes the repeated filtering that is already in bandstopfilter > > otherwise > > retrow(ME) > > end % case > > end % try-catch
Robert Oostenveld - 2013-02-20 15:12:31 +0100
I just implemented the suggestions for comment 5 and 6. mbp> git commit [bug780 b7cd0e8] enhancement - implemented the three variations for dealing with filter instabilities as suggested in http://bugzilla.fcdonders.nl/show_bug.cgi?id=780. Also added a test script. 1 files changed, 62 insertions(+), 0 deletions(-) create mode 100644 test/test_bug780.m I made the change in my (new) git version, so now have to figure out how to get them pushed to svn...
Vladimir Litvak - 2013-02-20 15:22:35 +0100
(In reply to comment #7) I don't see the commit in my SVN, but it's good you finally got to look at this. Please don't forget the scaling issues as well.
Robert Oostenveld - 2013-02-20 16:28:07 +0100
(In reply to comment #7) I forgot to commit some files mbp> git commit preproc/ [bug780 5c75c4f] enhancement - implemented the three variations for dealing with filter instabilities as suggested in http://bug illa.fcdonders.nl/show_bug.cgi?id=780 (part two of the commit) 5 files changed, 126 insertions(+), 53 deletions(-)
Robert Oostenveld - 2013-02-20 16:28:54 +0100
(In reply to comment #9) mbp> git checkout master mbp> git merge bug780 Updating 07622d0..5c75c4f Fast-forward preproc/ft_preproc_bandpassfilter.m | 32 +++++++++++++-- preproc/ft_preproc_bandstopfilter.m | 55 +++++++++++++------------- preproc/ft_preproc_highpassfilter.m | 48 ++++++++++++++++++----- preproc/ft_preproc_lowpassfilter.m | 30 ++++++++++++-- preproc/private/filter_with_correction.m | 14 +++---- test/test_bug780.m | 62 ++++++++++++++++++++++++++++++ 6 files changed, 188 insertions(+), 53 deletions(-) create mode 100644 test/test_bug780.m
Robert Oostenveld - 2013-02-20 16:30:52 +0100
(In reply to comment #10) mbp> git push Counting objects: 22, done. Delta compression using up to 2 threads. Compressing objects: 100% (13/13), done. Writing objects: 100% (13/13), 2.65 KiB, done. Total 13 (delta 11), reused 0 (delta 0) To git@github.com:oostenveld/fieldtrip.git 07622d0..5c75c4f master -> master this makes them visible on https://github.com/oostenveld/fieldtrip/commits/master PS sorry for the lengthy report on git, but I want to track all of this, as it is the first time to get some development from git back into svn. There will be more of these comments ...
Vladimir Litvak - 2013-02-20 16:31:07 +0100
(In reply to comment #10) Still don't see it. Have you switched to git now? Does it take time to sync with SVN?
Robert Oostenveld - 2013-02-20 16:53:44 +0100
(In reply to comment #12) > Does it take time to sync with SVN? Yes, because it takes some time for me to figure out with google how to do it ;-) We are (or I am) exploring git as alternative way of contributing. Some external collaborators prefer it. We are not planning to switch over completely, but if we can accept both (more loose) git and (more tight) svn contributions, then that might help in strengthening the development team. Since this bug is well defined and affects only some files, I considered it a good test case for the git->svn update direction.
Robert Oostenveld - 2013-02-20 18:00:48 +0100
(In reply to comment #11) After pushing them to github.com/oostenveld/fieldtrip, I had to get them into svn. mac001> git remote add oostenveld git@github.com:oostenveld/fieldtrip.git mac001> git pull oostenveld master this gets the changes from oostenveld/fieldtrip into my local fieldtrip-github-svn copy, i.e. the copy that is linked to the svn repository. mac001> git svn rebase mac001> git svn dcommit this pushes the changes to svn. http://code.google.com/p/fieldtrip/source/detail?r=7512 and http://code.google.com/p/fieldtrip/source/detail?r=7513 represent the two git commits in the svn version. And finally mac001> git push pushes the same changes to fieldtrip/fieldtrip. @Vladimir, you should now see them in svn.
Robert Oostenveld - 2013-03-08 12:00:38 +0100
there seems to be no further work related to this request.
Robert Oostenveld - 2013-03-17 09:25:54 +0100
I noticed defaults to be missing and the functionality not yet being available to the high-level user. allow the specification of instabilityfix for the low-level filters in the high-level cfg, provide documentation roboos@mentat001> svn commit Sending ft_preprocessing.m Sending preproc/ft_preproc_bandpassfilter.m Sending preproc/ft_preproc_bandstopfilter.m Sending preproc/ft_preproc_highpassfilter.m Sending preproc/ft_preproc_lowpassfilter.m Sending private/preproc.m Transmitting file data ...... Committed revision 7679. See also a few previous commits for the low-level code. Note that instead of 'none', the default behavior is now triggered by 'no'. That wording fits better in the documentation.
Robert Oostenveld - 2013-04-23 12:06:01 +0200
closed various bugs