[FieldTrip] FIR filter default order

[Andreas Widmann]

Hi Lucas,

> When performing a spectral estimation of EEG data by applying a FIR bandpass filter and then doing a Hilbert transform (ft_specest_hilbert), the default filter order proposed by ft_preproc_bandpassfilter is :
>     N = 3*fix(Fs / Fbp(1)); (l.236).
> where Fs is the sampling frequency and Fbp(1) the lower boundary of the frequency band.
> 
> I was simply wondering, where does this value comes from ? Any explication, litterature references or even advice on how to choose N (different from default value) when performing hilbert spectral estimation is welcome.

This equation is explicitly wrong for almost any application in electrophysiology, including yours (as I understand it). To my understanding this is long known. We discussed this issue some years ago in the context of a bug report. As far as I remember it was decided to keep it as it is for backward compatibility. Possibly, a warning notice that 'fir'-type should not be used without order parameter could be added?

All concepts underlying FIR filter order are explained here:
https://home.uni-leipzig.de/biocog/eprints/widmann_a2015jneuroscimeth250_34.pdf

In a nutshell:
For FIR filters we need a so called transition band between passband (no attenuation) and stopband (full attenuation) in which attenuation gradually changes from no to full attenuation. The width of this transition band is a function of filter order (and window type but this doesn’t matter for now). The higher the filter order the sharper the transition from passband to stopband, the lower the filter order the shallower the transition from passband to stopband.

The cutoff frequency is in the center of the transition band at half amplitude (-6 dB). In the first step you define your requested band edges, for example a passband from 9 to 11 Hz and stopbands below 7 and above 13 Hz. The cutoff frequencies (to be given as parameters in 'Fbp' as [Fhp Flp]) would be 8 and 12 Hz. Required  transition band width is 2 Hz (note that with the current implementation all transition bands must have equal width; different widths are in principle possible but complex to implement; it is simpler to successively apply lowpass and highpass filters with different orders; this approach has no relevant drawbacks besides computation time).

The 'fir‘ type filter by default uses a hamming window. The equation to compute filter order from requested transition bandwidth is (see Table 1 incl. caption of the above paper):

m = 3.3 / (df / fs)

3.3 is a window type specific parameter. df is transition bandwidth and fs is sampling frequency. For fs = 500 Hz this would result in a required filter order of m = 825. Filter order must be even, thus m = 826.

The equation used in 'fir' would give you an order of m = 186 in this example. But, order m = 186 would give you a transition bandwidth of df = 8.8 Hz (just solve the equation for df). That is, the passband of the highpass part would start at 8 Hz plus half the transition bandwidth df / 2 = 4.4, thus at 12.4 Hz and the passband of the lowpass part at 7.6 Hz (12 Hz - 4.4 Hz). As the higher passband edge is below the lower passband edge this obviously cannot work out; this filter does not have a passband.

I would advise against using the Fieldtrip 'fir‘ type filter for two more reasons: (1) It is based on the MATLAB fir1 function. This function is again based on firls *fitting* a filter to a rectangular response function. This is numerically considerably less accurate than computing the correct sinc function (by orders of magnitude at particular points as for example DC attenuation). (2) It defaults to twopass filtering. Twopass filtering is unnecessary for linear-phase FIR filters (as here), doubles filter order (but without actually much improving transition bandwidth) and doubles computation time.

In general I would recommend to use 'firws‘ filter type not showing these problems. The default for order is different from the above equation used in 'fir‘-type and more appropriate for typical ERP and related applications (wide passbands). In your case (narrow passbands; as I understood it) you will still need to give transition bandwidth as separate parameter 'df‘ (or compute order and give it as parameter ’N’).

Hope this help! Best,
Andreas

> Am 29.05.2020 um 08:35 schrieb Lucas Struber <lucas.struber at univ-grenoble-alpes.fr>:
> 
> Hello,
> 
> When performing a spectral estimation of EEG data by applying a FIR bandpass filter and then doing a Hilbert transform (ft_specest_hilbert), the default filter order proposed by ft_preproc_bandpassfilter is :
>     N = 3*fix(Fs / Fbp(1)); (l.236).
> where Fs is the sampling frequency and Fbp(1) the lower boundary of the frequency band.
> 
> I was simply wondering, where does this value comes from ? Any explication, litterature references or even advice on how to choose N (different from default value) when performing hilbert spectral estimation is welcome.
> 
> For information, i'm working on EEG data acquired in human subjects in a visuo-motor adpatation task.
> 
> Kind regards,
> 
> -- 
> Lucas Struber
> Ingénieur de recherche
> +33 (0)6.33.91.93.52
> 
> Laboratoire TIMC-IMAG, équipe SPM
> UGA Site Santé
> 38700 La Tronche
> 
> _______________________________________________
> fieldtrip mailing list
> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
> https://doi.org/10.1371/journal.pcbi.1002202