Change to a filtering with Morlet wavelet?

[TomDLT]

I would like to question your choice of filtering method.

In your firf method, you both specify a bandwidth (i.e. two cut-off frequencies) AND a number of cycles. It makes the use of this function quite difficult, since if the number of cycles is not high enough, the obtained bandwidth may not be the one expected:

Example: with f_range=[5, 15], fs=300.0, w=3, the obtained filter has a bandwidth (at -3dB) around 10 Hz (expected). with f_range=[45, 55], fs=300.0, w=3, the obtained filter has a bandwidth (at -3dB) around 30 Hz (unexpected).

This is the normal behavior since you are based on scipy.signal.firwin, so the obtained bandwidth is somehow around the maximum of (a) the expected bandwidth and (b) the bandwidth obtain with a Morlet wavelet with the same number of cycles. However, this is not documented and it is quite difficult to know which parameter to use without making mistakes (as in example above).

I wonder if it would not be better to use a filtering based on a Morlet wavelet. In this case, you would have to specify only the center frequency and the number of cycles, and the bandwidth would derive directly from these two parameters.

If you don't want to keep the f_range parameter (with two cutoff frequencies), you could deprecate the number of cycles. In this case, you could use Morlet wavelet filtering with:

# center frequency
f0 = 0.5 * (f_range[1] - f_range[0])
# bandwidth
bandwidth = f_range[1] - f_range[0]
# number of cycles
w = 1.7 * f0 / bandwidth
# where 1.7 is chosen to have the desired bandwidth at -3dB
x_filt = morletf(x, f0, fs, w)

It is also more intuitive to use bandwidth than number of cycles.

The bandwidth choice for PAC is a critical issue [1], so I believe it is important to help this choice by avoiding adding more pitfalls.

[1] Aru, J., Aru, J., Priesemann, V., Wibral, M., Lana, L., Pipa, G., ... & Vicente, R. (2015). Untangling cross-frequency coupling in neuroscience. Current opinion in neurobiology, 31, 51-61.

[srcole]

I agree with all of your points. The filtering could definitely be improved, as I've been realizing over the past few months. For the FIR filter, I should have just made the input to be the filter order. I think the approach you outlined to derive f0 and bandwidth from an f_range is a good way to go. We, @parenthetical-e , are hesitant to deprecate the current API (w in this case), but I personally think it would be worth it to improve the filter design, which as you pointed out, could work more intuitively.

[parenthetical-e]

Regardless of whether we shift to the Morlett by default (which could be good) It's clear we're going to need to revisit the filter design and filter API of pacpy. ASAP. Thanks for pointing this out, @TomDLT!

[srcole]

I applied the new morlet filtering in pacpy in the branch, newfilter. This will require writing new answers for all of the pac tests. I'm going through some of my own sanity checks now.

Sorry for the pacpy.egg files. I think those showed up while I was learning to use pypi on my computer and I accidentally pushed them with this branch.

[srcole]

The frequency responses of the new narrowband filters used in the comodulogram method look good. E.g. 96-100Hz:

However, the filtering for broadband frequency ranges, e.g. beta (13-30Hz) and high gamma (50-200Hz), don't look as good. Perhaps for the non-comodulogram functions, we should stick to window FIR filters and make the user declare a filter order (instead of number of cycles, which was admittedly a bad idea for this before).

Morlet wavelet 50-200Hz:

Window FIR filter 50-200Hz, order 240:

[TomDLT]

You should also plot the time response, since a good filter in frequency is usually not so good in time (i.e. very long). cf. recent tutorial in MNE.

[srcole]

Yup, sorry didn't mention that. The narrow band-pass morlet filters are quite long (e.g. 850ms for the 96-100Hz bandpass). I accidentally made that last FIR filter really long, but a bandpass 50-200Hz looks pretty good with order 240ms.

Right now I think the best default filtering would be Morlet wavelets (with length determined as Tom described) for comodulogram and FIR window for everything else in which the user inputs the order. However, I think this would now require 2 filter_kwargs keywords: one for the low frequency range and one for the high frequency range. I'm wondering if instead of having one 'PAC' function if we should break out the filtering to its own function and then feed the filtered signals into the PAC functions.

[srcole]

I didn't fully realize how much the length of the filter would mess up the PAC estimate. My new thoughts are below.

Since our new Morlet wavelets are so long, an increase in high-frequency amplitude at one time point gets smeared across many low-frequency cycles. Similarly, the long filter duration will cause issues with the phase estimate by overestimating the stationarity of these oscillations.

e.g. a Morlet filter with cutoff frequencies at 120Hz and 125Hz has a filter length of 680

And this causes a very low-frequency amplitude envelope of this high frequency component. Original signal is in red, the signal morlet-filtered 120-125Hz cutoff frequencies is shown in black

In terms of filtering for regular PAC calculation (outside of the comodulogram function), these may be broadband (e.g. beta 13-30Hz) or narrowband (e.g. 20-22Hz) frequency ranges. So the best time-frequency tradeoff likely changes in different situations.

If we want to use a Morlet wavelet, some options could be:

1. User chooses frequency range. We design a morlet wavelet with the cutoff frequencies at -3dB. We tell the user what the filter length is and warn that if it is too long then their PAC estimate will not be accurate. The user can suppress this output.
2. User chooses a center frequency and a number of cycles. We tell the user what the true -3dB cutoff frequency and filter length are. We may also tell them the filter length (even though they can calculate this) just as a warning.

If we decide to keep using the window FIR filter, some options could be:

1. Force the user to choose a filter order in addition to cutoff frequencies (which are both of the inputs to scipy.signal.firwin, so they have full control. Tell them the true -3dB cutoff frequencies.
2. Only force the user to choose cutoff frequencies and define the filter length in one of the ways below. 2a. Define filter order based on maintaining the -3dB cutoff frequencies. Warn them of the filter length. 2b. Treat low frequency and high frequency filters differently. Such as using (2a) for the low frequency filter but determining the order of the high frequency filter based on the low frequency signal (e.g. 1 low-frequency cycle). Then tell them the true cutoff frequencies.

We should further think specifically about the comodulogram method which should automate many filters.

1. Choose Morlet filter center frequencies and number of cycles (default, e.g. 7). We could optional output an array of the true cutoff frequencies for each range. For FIR filters, they cannot choose the order of each filter.
2. Hold the filter length constant for all amplitude frequencies so that bandwidth increases in a 'natural' manner.
3. Same as (2) but we could make the filter length dependent on the phase frequency for which PAC is being calculated (as in 2b above), but this would take much longer as we would have to do multiple filters for each amplitude frequency range.

For non-comodulogram PAC, I'm partial to the FIR filter option (2a) in which they define the frequency range and we warn them of the filter length. For comodulgoram PAC, I think the best bet will be (1) in which they choose the Morlet filter center frequencies and number of cycles. Any other suggestions?

[parenthetical-e]

@roemervandermeij I'd appreciate your input on this thread.

[TomDLT]

For clarity: HFO = High Frequency Oscillation LFO = Low Frequency Oscillation

Filtering I prefer the Morlet wavelet since the firwin-filter's limitations are more difficult to understand. Giving full control of the firwin-filter would be necessary, but only experienced signal processing users would know how to use it properly.

I also tend to prefer the frequency range parameter (1), since this is more intuitive than the number of cycles (2). It also guides the users toward a constant bandwidth when they change the center frequency, whereas specifying the number of cycles guides them toward a variable bandwidth. Since the HFO's bandwith has to be at least twice the LFO's frequency (!), I strongly recommand using a constant bandwidth, so I prefer the bandwidth parameter (or a frequency range).

Comodulogram As suggested by Chris in #45, the priority would be to reuse the filtered signals.

About how to parametrize the comodulogram's filtering, I again prefer your choice (2), i.e. having a constant bandwidth. The proper way would be (3), but as you said it is much slower.

[srcole]

Non-comodulogram Filtering I don't think morlet will work here because it is not optimal compared to firwin when an arbitrary window is chosen (e.g. 'beta' 13-30Hz or 'high gamma' 50-200Hz). It seems to work much better if you don't have a specific frequency range but rather a rough idea of its parameters, # cycles and center frequency. I agree with you, keeping the filter parameter as a frequency range is more intuitive and preferable. Here is a morlet filter that I designed using sp.signal.morlet with the w parameter defined as in your original post and the standard M = 2 * w * fs / f0.

High gamma, 50-200Hz (f0=125, bw=150, w=1.41666, M (order) = 22) Obviously, we need a longer filter to capture this frequency range better, but since a Morlet filter is designed to have a Gaussian frequency response, it cannot accomplish this band-pass like the FIR filter below (firwin, 50-200Hz, order=240)

Comodulogram I don't think a constant bandwidth is feasible (my option (2) was a constant filter order, not bandwidth). If we keep the bandwidth constant in a comodulogram ranging from 38-42Hz to 198-202Hz, the difference in filter length is increased 5-fold (whereas for the higher frequencies, the filter length should be decreased), and so there will be much more smoothing of the high frequency component across time. This will lead to different behaviors in PAC across these frequencies, which would be very bad. Therefore, I think the best option is using morlet filters with a constant number of cycles, stepping the center frequency.

[roemervandermeij]

Some thoughts from my side.

1) FIR filtering and Morlet filtering (or wavelet convolution) can be used to achieve different things, and which one is most desired depends on the application. I think choosing between the two as a default doesn't make a lot of sense (i.e. I'd argue for no default to force user choice/thought). If one has to be chosen, for reason number 2 below I think FIR is a more appropriate one.

2) Morlet filtering should not be possible to apply with two-pass filtering. Morlet's are most commonly thought of in the context of wavelet convolution, and arguable will always be applied in a one-pass manner. Unexpected two-pass (with additional smearing and an oddly shaped amplitude profile) will be very surprising, and not something to think of for users.

3) Arguable by far the most important parameter in spectral analysis, especially in the context for cross-frequency analyses, is the length of the time-domain information that is used. This is the reason IIR filtering should never be used in this context unless with great care. The length of the filter/wavelet should be a key parameter the user has to think about, or should be a sensible default (like 3 cycles, 5 cycles, etc). I strongly recommend this to never be set automatically based on the desired frequency response.

4) Tom brought up an important issue, which, since the toolbox is pretty mature, can be handled more cleverly. The nice thing about Morlet wavelets for filtering/convolution is that they have a very predictable frequency response, due to the fact they are specified by the number of cycles and center frequency only (usually at least), and are then nothing more than a sine/cosine multiplied with a Gaussian. A FIR filter is more difficult, and more flexible, because it allows you to modify the shape of the response to be e.g. more broadband/narrowband than a Morlet. What you specify depends a bit on the implementation (I'm not familiar with the one you use) but it's usually the time-domain length combined with the cut-off frequencies combined with the transition width (specified as the width at which the amplitude becomes X dB). (These parameters are dependent on each other). The problem is then, that not all of the input space is reasonable. E.g. you cannot perform a 15-30Hz bandpass with a filter that is 10ms long. But you can request it. And not everybody knows that this won't work, so some user protection is warranted.

There are different options, but I'd encourage you to do the following one (reasoned from point 3, and one that's broadly accepted): 4a) specify the time-domain length and the cut-off frequencies, and default the transition width 4b) give a warning/error when the resulting transition width result is undesirable

An example of this defaulting and warnings can be found in the following code, the filter being 'firws' (line 150). This is nearly identical to the way EEGLab does it currently (implemented by same person). It creates a FIR filter using the window approach, but that doesn't matter, the key part is the checking of the transition width at line 167 and the warning when a reasonable one cannot be achieved. https://github.com/fieldtrip/fieldtrip/blob/master/preproc/ft_preproc_bandpassfilter.m#L150 https://github.com/fieldtrip/fieldtrip/blob/master/preproc/private/invfirwsord.m https://github.com/fieldtrip/fieldtrip/blob/master/preproc/private/fir_df.m

[TomDLT]

About choosing or not choosing FIR(scipy.signal.firwin) vs Morlet:

• First, I don't know if forcing the user to choose is really helping, since it does not immediately change users into experts. And FIR are hard to set, even for experienced users.
• Second, in this PAC package, filtering is only used before applying a Hilbert Transform to estimate the phase/amplitude of the real-valued signal. As the Hilbert Transform should ONLY be applied on narrow-band signals _[1], I would not consider the wide-band limitations of the Morlet-based filtering relevant. Computing the PAC on the frequency range [50, 200] Hz arguably makes no sense. So I tend to prefer a Morlet-based filtering by default.

_[1] Chavez, M., et al. "Towards a proper estimation of phase synchronization from time series." Journal of neuroscience methods 154.1 (2006): 149-160.

About two-pass filtering with Morlet-based filters: I don't understand what is wrong with two-pass filtering, and if it's a problem, I don't see why we can't just have one-pass filtering. If I understand correctly, the main purpose of two-pass filtering is to have a zero-phase filter (which is crucial to PAC analysis). But if we make sure than the Morlet-based filter is symmetric and odd in length, even the one-pass filter is zero-phase.

About the time length of the filter: I totally agree that the time length is very important (and should not be too long), which is why I tend to avoid FIR filters that have very good precision in frequency but which are very long in time. This is also why I proposed to use the bandwidth as a intuitive and reasonnable parameter:

• The bandwidth is easy to understand for a non-expert in signal processing.
• The length of the Morlet filter (aka the order) is fixed with the bandwidth: order = 2 * w * fs / fa = 2 * fs * 1.7 / bandwidth is independant of fa the frequency of the amplitude signal. In other words, a fixed bandwidth and a fixed filter length are equivalent.
• To be able to see any PAC, we have the condition bandwidth > 2 * fp where fp is the low modulation frequency _[2]. A bandwidth parameter makes this condition easier to set. Moreover, if the comodulogram is not computed with a fixed bandwidth, it adds a bias in the low fa, since the condition is not fullfiled everywhere in the comodulogram.

_[2] Aru, J., Aru, J., Priesemann, V., Wibral, M., Lana, L., Pipa, G., ... & Vicente, R. (2015). Untangling cross-frequency coupling in neuroscience. Current opinion in neurobiology, 31, 51-61.

@srcole

Comodulogram I don't think a constant bandwidth is feasible (my option (2) was a constant filter order, not bandwidth). If we keep the bandwidth constant in a comodulogram ranging from 38-42Hz to 198-202Hz, the difference in filter length is increased 5-fold (whereas for the higher frequencies, the filter length should be decreased), and so there will be much more smoothing of the high frequency component across time. This will lead to different behaviors in PAC across these frequencies, which would be very bad. Therefore, I think the best option is using morlet filters with a constant number of cycles, stepping the center frequency.

This is exactly my point, up to the difference that as far as I understand it, a fixed number of cycles makes the filter length change with the center frequency. Whereas a fixed bandwidth makes the filter length be fixed. A constant filter order IS a constant bandwidth.

[roemervandermeij]

@TomDLT

About two-pass filtering with Morlet-based filters: I don't understand what is wrong with two-pass filtering, and if it's a problem, I don't see why we can't just have one-pass filtering. If I understand correctly, the main purpose of two-pass filtering is to have a zero-phase filter (which is crucial to PAC analysis). But if we make sure than the Morlet-based filter is symmetric and odd in length, even the one-pass filter is zero-phase.

The problem with two-pass filtering, for the purpose of extracting a signal to further analyze whose time-domain smoothing is important, is that the signal is attenuated in amplitude twice and temporally 'spread' twice as far. Wether the former is undesirable is debatable. But the latter is in nearly all usage scenario's undesirable, as this effectively doubles the filter length without an increase in spectral precision. Because of this, signal coming from samples a full filter length away on both sides impact every sample, instead of only half. You are correct in that the filters used here can easily be zero-phase filters with one pass. I think Scott is implementing this right now as well.

[srcole]

tl;dr I'll implement the windowed sinc filter with specifying order and cutoff frequency parameters

But if we make sure than the Morlet-based filter is symmetric and odd in length, even the one-pass filter is zero-phase.

Right, I'm doing one-pass filtering now, but I didn't make sure the filter length was odd. I'll change that if we choose to use morlet.

4a) specify the time-domain length and the cut-off frequencies, and default the transition width 4b) give a warning/error when the resulting transition width result is undesirable

This sounds like a good choice to me. Thanks for pointing to the fieldtrip code!

First, I don't know if forcing the user to choose is really helping, since it does not immediately change users into experts.

Good point, but I think this is something the user should pick since there is no right answer for us to choose. I think having this as a required parameter will make the user appreciate how the filter length can affect the calculation. And maybe we can provide some simple function they can call to visualize the kernel and frequency response. I think it's reasonable to stick with FieldTrip and EEGLAB's filter standard here. I'm not familiar with how MNE handles this though.

As the Hilbert Transform should ONLY be applied on narrow-band signals _[1], I would not consider the wide-band limitations of the Morlet-based filtering relevant. Computing the PAC on the frequency range [50, 200] Hz arguably makes no sense.

I get what you're saying, but there are people in this field (including us) who use broadband ranges when calculating coupling (e.g. High gamma, ~50-200Hz [1,2,3]) partially since broadband power has been associated with local spiking activity [4]. Another example would be 'beta' which is not very narrowband at its typical quoted 13-30Hz. While choosing these semi-arbitary frequency ranges is arguably not the best approach for this (and may change in the future), I think it would be best for this package to accommodate these. [1] Canolty, R.T., Edwards, E., Dalal, S.S., Soltani, M., Nagarajan, S.S., Berger, M.S., Barbaro, N.M., Knight, R.T., 2006. High Gamma Power is Phase-Locked to Theta Oscillations in Human Neocortex 313, 1626–1628. doi:10.1126/science.1128115.High [2] de Hemptinne, C., Ryapolova-Webb, E.S., Air, E.L., Garcia, P.A., Miller, K.J., Ojemann, J.G., Ostrem, J.L., Galifianakis, N.B., Starr, P.A., 2013. Exaggerated phase-amplitude coupling in the primary motor cortex in Parkinson disease. Proc. Natl. Acad. Sci. U. S. A. 110, 4780–5. doi:10.1073/pnas.1214546110 [3] Voytek, B., Kramer, M. A., Case, J., Lepage, K. Q., Tempesta, Z. R., Knight, R. T., & Gazzaley, A. (2015). Age-related changes in 1/f neural electrophysiological noise. The Journal of Neuroscience, 35(38), 13257-13265. [4] Manning, J.R., Jacobs, J., Fried, I., Kahana, M.J., 2009. Broadband Shifts in Local Field Potential Power Spectra Are Correlated with Single-Neuron Spiking in Humans. J. Neurosci. 29, 13613–13620. doi:10.1523/JNEUROSCI.2041-09.2009

A constant filter order IS a constant bandwidth.

Right. My bad. There was a bug in my code and I was getting a change in order.

To be able to see any PAC, we have the condition bandwidth > 2 * fp where fp is the low modulation frequency _[2]

So given this, in the comodulogram method when we are scrolling across fp, should we be changing the order/bandwidth of fa? e.g. for fp=30-32Hz, all of the amplitude bandwidths should be >=60Hz? In that case, the cutoff frequencies for fa would need to change as a function of fp. This could be a good thing, but I'm not familiar with anyone calculating the comodulogram in this manner.

[roemervandermeij]

To be able to see any PAC, we have the condition bandwidth > 2 * fp where fp is the low modulation frequency _[2]

So given this, in the comodulogram method when we are scrolling across fp, should we be changing the order/bandwidth of fa? e.g. for fp=30-32Hz, all of the amplitude bandwidths should be >=60Hz? In that case, the cutoff frequencies for fa would need to change as a function of fp. This could be a good thing, but I'm not familiar with anyone calculating the comodulogram in this manner.

Although I agree with you Tom that, in general, the modulating frequency should be twice as low as the modulated frequency, there is a lot of value to visualize PAC metrics into the "bad" ratios. Usually, they are a good sign of either mistakes having been made, or certain artifacts being present in the time-series. A common example in intracranial data is a small (easy to miss) fast spike riding on a slow wave. This will lead to high PAC between the phase of high frequencies and the amplitude of lower frequencies, and also PAC at the reverse ratio. As such spikelets are usually not the phenomenon of interest, it is good to know they occurred so that we can remove them.