Usage of parameter 'SigmaFreqs'

[manimoh]

First of all, I must congratulate you and your colleagues on this very elegant work on estimating phases in Real-Time! And Kudos to your team for making the code available as open source. :) I was trying to use the SSPE method for offline phase estimation and looked at the demo.m. My understanding is that freqs must be the central frequency of the oscillation that we are interested in, say assuming theta is 6 - 11Hz, the central frequency will be 8.5 Hz. But It is not clear to me how to initialize the SigmaFreqs. For instance, if one were interested in the high gamma range (60-100 Hz) then the freqs will be 80, but will SigmaFreqs be the same as what we might set for theta?

[wodeyara]

Thanks for the interest! So, my suggestion would be to set SigmaFreqs to the same as the theta SigmaFreqs as a first pass fitting attempt, and if there is no oscillator identified in the gamma range, attempt scaling the SigmaFreqs down an order of magnitude. So if you started out with SigmaFreqs_Gamma = 1, then try again at SigmaFreqs = 0.1 and then SigmaFreqs = 0.01. EM can be sensitive to the initial start conditions but it operates well once you're in the same magnitude ballpark. Let me know if that works!

On Wed, Apr 21, 2021 at 6:20 PM Manish Mohapatra @.***> wrote:

First of all, I must congratulate you and your colleagues on this very elegant work on estimating phases in Real-Time! And Kudos to your team for making the code available as open source. :) I was trying to use the SSPE method for offline phase estimation and looked at the demo.m https://github.com/wodeyara/stateSpacePhasePredictor/blob/master/demo.m. My understanding is that freqs must be the central frequency of the oscillation that we are interested in, say assuming theta is 6 - 11Hz, the central frequency will be 8.5 Hz. But It is not clear to me how to initialize the SigmaFreqs. For instance, if one were interested in the high gamma range (60-100 Hz) then the freqs will be 80, but will SigmaFreqs be the same as what we might set for theta?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/wodeyara/stateSpacePhasePredictor/issues/1, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMYNTOJ3OR6YZP3BHDYV65DTJ5FOFANCNFSM43LHMPAQ .

[manimoh]

Actually I couldn't get a good fit by varying the range of SigmaFreqs from 10 to 1e-9, for any of the frequency ranges. A message that was displayed in each case was the following: Low freq band limits incorrect OR there is no low freq signal; retaining initial params To give you a little more context, each trial for me is about 2.5 seconds and I am using half of the trial length as the window size. I tried this for data with sampling rate of 2.6 kHz and another one with a rate of 32 kHz, but found no success in either of them.

Do you think intializeParams.m can help? I am not sure about what num_components means and didn't quite understand how to interpret the output of that function.

[wodeyara]

That's unfortunate, and seems to suggest that the theta itself is difficult to identify in the data. Some thoughts:

1. I'd suggest checking if there's any drift in the data, and even if not bandpass filter the data between 0.5 to 200 Hz, this will help avoid issues with low frequencies power overwhelming the frequencies identified
2. Check to ensure that the segment that is being passed to fit parameters has a visible theta (which I assume is the low frequency band of interest?), I've found that when the peaks are low SNR (<2) then the algorithm struggles to identify an oscillator
3. Try using a longer data segment for fitting the parameters relative to the segment used for real-time analysis, will take a little futzing with the causalPhaseEM_MKmdl.m to ensure that it uses a longer segment for parameter fitting (I can help with this if you want)
4. Using more samples doesn't really help actually, and can actually hurt fitting; this has to do I think with the scaling parameter a which needs to be closer and closer to 1 to sustain a rhythm as there's more samples
5. Finally, it's actually alright if it chooses to retain initial params as long as you believe the parameters are in a reasonable range (to confirm I would plot the theta estimated in the state vector against the data to visually ensure it is capturing what you intended). So reasonable range in my mind means: Central freqs set at the peaks in PSD, sigmaFreqs set to the square root of the amplitude, with high freq (gamma) an order of magnitude smaller than low freq (theta), ampVec = 0.99 and sigmaObs =1.

I'll make a readme that has a little segment dedicated to parsing through this soon. Let me know if this helps!

On Thu, Apr 22, 2021 at 9:53 PM Manish Mohapatra @.***> wrote:

Actually I couldn't get a good fit by varying the range of SigmaFreqs from 10 to 1e-9, for any of the frequency ranges. A message that was displayed in each case was the following: Low freq band limits incorrect OR there is no low freq signal; retaining initial params To give you a little more context, each trial for me is about 2.5 seconds and I am using half of the trial length as the window size. I tried this for data with sampling rate of 2.6 kHz and another one with a rate of 32 kHz, but found no success in either of them.

Do you think intializeParams.m https://github.com/wodeyara/stateSpacePhasePredictor/blob/master/initializeParams.m can help? I am not sure about what num_components means and didn't quite understand how to interpret the output of that function.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/wodeyara/stateSpacePhasePredictor/issues/1#issuecomment-825326194, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMYNTOLAEUINOD272AM4WKLTKDHE7ANCNFSM43LHMPAQ .

[manimoh]

I think some more background about what I intend to do might be helpful here: I wanted to estimate the phase of various oscillations (delta, theta, beta, alpha, low gamma and high gamma) that might be ongoing in the LFP at the time of an incoming stimulus. The inter-stimulus interval is 3 sec, and therefore to avoid any contamination of the LFP due to the stimulus, I planned on using the time window of 'stimulus-2.5 sec: stimulus' to estimate the phase at the final bin of this window. This 2.5 sec value could be further increased to maybe 2.8 because the stimulus artifact might not be longer than 0.2 sec. Matlab's filtfilt or other acausal filters are not suitable for this, and the SSPE method seemed ideal for this purpose because of its forward prediction nature. So what I am currently doing is I run the SSPE for each of the frequency ranges of interest. Here's an example piece of how I set the parameters: f_list = {[2 5], [6 10],[15 25], [30 40],[40 60], [60 80]}; for iF = 1:length(f_list) cfg_sspe.freqs = mean(f_list{iF}); % initialization using above not great at identifying starting freq cfg_sspe.Fs = Fs; cfg_sspe.ampVec = 0.99; cfg_sspe.sigmaFreqs = 10^(1-iF); cfg_sspe.sigmaObs = 1; cfg_sspe.window = length(this_data)/2; cfg_sspe.lowFreqBand = flist{iF}; [phase,phaseBounds, fullX] = causalPhaseEM_MKmdl(this_data, cfg_sspe); % do rest of the stuff end I do not know apriori, if the oscillations that I am interested in are ongoing in the time window of my interest. From what you suggested in the previous reply, I guess I can do a PSD of this time window and look for peaks in each range, and if found initialize the freqs and sigmaFreqs accordingly. Do you think I am using SSPE appropriately, or am I missing something here?

[wodeyara]

So the issue here I think is that there is an attempt to fit an oscillator without checking where the rhythms are. If the rhythms are of sufficient SNR they will overcome the hurdle of an incorrect starting parameter, but if not, then parameters will fail to converge entirely. Also, if some bands don't have a fit, there are two potential interpretations:1) EM got stuck in a local minima, or, 2) the data simply doesn't have the evidence to support the existence of activity in that band.

I do think the optimal way to fit parameters is to guesstimate the location of the primary rhythms present (either from the PSD or a spectrogram) and to estimate the model with those, rather than the fully data agnostic approach based on predefined bands as you are trying right now. As for filtering - you could filter a larger segment of time that includes that stimulus periods and drop those periods before using that data for fitting the model/simply filter the entire time series before epoching. Finally, given your use case, I would suggest using the output of the fitted model directly (by using the stateVec from fit_MKModel_multSines) rather than the output of a real-time phase estimate from causalPhaseEM_MKmdl.

As for your particular use case, an extension I've been thinking about but not set up is one that minimizes error across trials so that you could potentially combine information about the parameters across trials but I haven't formalized the idea yet.

[manimoh]

Thank you for your suggestions, let me try that real quick. Meanwhile, I had a question about the stateVec- It is a 2 x windowsize array, where stateVec(1,:) looks very similar but not exactly the same as y(1:windowsize-1) and stateVec(2,:) looks like an estimated signal. Am I correct in my interpretation?

[wodeyara]

Yes, stateVeccan be interpreted the same way as x in the causalPhaseEM_MKmdl. The first index will be the real part of the rhythm and the second will be the imaginary part so that phase can be estimated (for the first rhythm) as angle(stateVec( 1, :) + 1i*stateVec( 2, :)).

[manimoh]

Unfortunately, I wasn't able to use fit_MKModel_multSines to estimate the oscillation even when I can visually detect theta in the raw LFP and using the parameters as you suggested. This got me thinking about whether I am doing the right thing by looking for a particular frequency band when there are multiple rhythms present. As a test, i modified the demo.m slightly to add another interfering oscillation of frequency 50 Hz and keeping the rest of the code the same

Vlo = (10).*cos(2*pi*(6).*[1/Fs:1/Fs:time]); % random movement in amplitude, creates a little bump in PSD V1 = (2.5).*cos(2*pi*(50).*[1/Fs:1/Fs:time] + pi/3); [pn] = make_pink_noise(1,1e4,1/Fs); pn = 10*pn; data = Vlo + V1; ...

Here's a snapshot from the two runs of demo (6 Hz + 50 Hz + pink noise on the left, 6 Hz + pink noise on the right):

The error estimates increased and as you can see in the zoomed in picture, the phase estimates get wildly different. I tried increasing the window size for this interference signal example, but that didn't improve the estimates.

[wodeyara]

Yes, so what is happening there is that without a 50 Hz oscillator, the state vector is compensating for the presence of a 50 Hz oscillations by oscillating the phase of the slower oscillation to a degree. Instead, if you include a second oscillator by using as initial parameters:

initParams.freqs = [4,45]; % initialization using above not great at identifying starting freq
initParams.Fs = 1000;
initParams.ampVec = [.99,.99]; % in a pinch this can be initialized to 0.99 to start
initParams.sigmaFreqs = [10,1]; % its important to use a value of this that is in the same ballpark scale
initParams.sigmaObs = 1;
initParams.window = 5000;
initParams.lowFreqBand = [4,8];

You should be able to recover the plot below.

Nonetheless, this is a fairly ideal situation since the error is kept to a very controlled degree . If the 50 Hz is below the level of the noise, the EM struggles to identify a genuine oscillator in that range.

[manimoh]

I wanted to estimate the phase of various oscillations (delta, theta, beta, alpha, low gamma and high gamma) that might be ongoing in the LFP at the time of an incoming stimulus. The inter-stimulus interval is 3 sec, and therefore to avoid any contamination of the LFP due to the stimulus, I planned on using the time window of 'stimulus-2.5 sec: stimulus' to estimate the phase at the final bin of this window. This 2.5 sec value could be further increased to maybe 2.8 because the stimulus artifact might not be longer than 0.2 sec. Matlab's filtfilt or other acausal filters are not suitable for this, and the SSPE method seemed ideal for this purpose because of its forward prediction nature. So what I am currently doing is I run the SSPE for each of the frequency ranges of interest. Here's an example piece of how I set the parameters: f_list = {[2 5], [6 10],[15 25], [30 40],[40 60], [60 80]}; for iF = 1:length(f_list) cfg_sspe.freqs = mean(f_list{iF}); % initialization using above not great at identifying starting freq cfg_sspe.Fs = Fs; cfg_sspe.ampVec = 0.99; cfg_sspe.sigmaFreqs = 10^(1-iF); cfg_sspe.sigmaObs = 1; cfg_sspe.window = length(this_data)/2; cfg_sspe.lowFreqBand = flist{iF}; [phase,phaseBounds, fullX] = causalPhaseEM_MKmdl(this_data, cfg_sspe); % do rest of the stuff end I do not know apriori, if the oscillations that I am interested in are ongoing in the time window of my interest. From what you suggested in the previous reply, I guess I can do a PSD of this time window and look for peaks in each range, and if found initialize the freqs and sigmaFreqs accordingly. Do you think I am using SSPE appropriately, or am I missing something here?

So in that case what I did above was incorrect because if there were oscillations in all the present, because I need to have a single run with all the oscillators and not separate runs for each. But then for which rhythm are the phase values returned when you pass multiple oscillators in the parameters? Based on the parameters you mentioned in your previous comment, the plots look like it is the phase of the 6Hz oscillator. What if I want to recover the phase of a particular rhythm (say maybe that of the 50 Hz one)?

[wodeyara]

Thanks for all the questions and for following through on this by the way! Helping me figure out how to improve on the current setup considerably :)

Ok so. About your first question, you're right, I misread your code and assumed you were trying to fit a distinct oscillator in each band rather than each one separately. So yes, you would need to set up initialparamswith all the rhythms you expect to want to track. The algorithm returns the phase values for the oscillator that exists in the lowFreqBand (if there's multiple, then the one closest to the center of the band gets chosen).

So if you wanted the phase for the 50 Hz oscillator, then you could set: initParams.lowFreqBand = [45,55];

Alternatively, the algorithm returns fullX which tells you the estimated real/imag parts for each oscillator, but I'm realizing we would want to know their center frequencies so I will change the code to return the results for the estimated parameters as well.