[FieldTrip] FIR filter default order

Fri May 29 16:01:32 CEST 2020

Hi Lucas,

> Thank you for your answer, which confirm what I understood about FIR filters (I read your paper before you mentionned it).
> 
> However I still can't explain myself why this *wrong* default parameters gives me spectral amplitudes very similar to spectral amplitudes estimated with a Morlet wavelet transform (with n cycles equal to 7). While using FIR and setting m = 3.3/(df/fs) (m is then the same for all frequencies), gives me a very different result...
> 
> Any insights or suggestions on this ?
There are various intricacies and versions of amplitude scaling in wavelet analysis (e.g. normalized frequency as in textbook equations vs. sampling frequency in the application etc.). My suggestion would be to start debugging with synthetic test signals where you know the ground truth, for example impulses, sine waves and white noise. 

Morlet wavelet analysis can be implemented using convolution, that is, Morlet wavelets can in principle be (ab-)used as FIR filters (as STFT can be implemented via convolution; the only major difference between Gabor STFT transform and Morlet wavelet transform is the definition of time zero).

And you can compare your FIR filter coefficients (returned by ft_preproc_*passfilter.m in B) and your Morlet wavelets (real+imag) in the frequency domain (the example code in the supplementary materials of the paper should contain all necessary parts if I remember correctly).

Best,
Andreas

> In any case i'll try with firws.
> 
> Best regards,
> 
> Lucas
> 
> Le 29/05/2020 à 13:11, Andreas Widmann a écrit :
>> 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
>> 
>> _______________________________________________
>> fieldtrip mailing list
>> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>> https://doi.org/10.1371/journal.pcbi.1002202
> 
> -- 
> 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