[MRG] Read_raw extension to extend sidecar json func.

[adam2392]

PR Description

Extends read_raw functionality to extract data from the sidecar json file. Closes: #340

Also includes check on data elements, such as:

• Reference
• SamplingFrequency

Merge checklist

Maintainer, please confirm the following before merging:

• All comments resolved
• This is not your own PR
• All CIs are happy
• PR title starts with [MRG]
• whats_new.rst is updated
• PR description includes phrase "closes <#issue-number>"
• Commit history does not contain any merge commits

[jasmainak]

any test to update?

[sappelhoff]

Thanks for the contributions. I made some comments.

Also:

• What happens if we CAN read a line_freq from the JSON, and the raw.info[line_freq] is None? --> In that case it'd probably be good to just set the raw.info[line_freq] according to what we read, instead of raising an error

[adam2392]

any test to update?

I checked test_read.py, but there were not specific unit tests for read_raw_bids. Do you want me to write an additional test? If so, what should it look like?

[jasmainak]

see how we do the rest of the tests. Check if the errors are raised etc.

also you need to rebase to make the CIs happy. Good luck! :)

[codecov-io]

Codecov Report

Merging #341 into master will increase coverage by 0.23%.
The diff coverage is 98.36%.

@@            Coverage Diff             @@
##           master     #341      +/-   ##
==========================================
+ Coverage    94.9%   95.14%   +0.23%     
==========================================
  Files          11       11              
  Lines        1218     1277      +59     
==========================================
+ Hits         1156     1215      +59     
  Misses         62       62

Impacted Files	Coverage Δ
mne_bids/write.py	97.27% <100%> (+0.2%)	⬆️
mne_bids/read.py	99.41% <100%> (+0.07%)	⬆️
mne_bids/utils.py	98.3% <97.43%> (-0.14%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a54d616...d82ccb7. Read the comment docs.

[jasmainak]

@adam2392 can you see the checklist above? You need to say "closes #340" or "fixes #340". I don't think "addresses" will do it.

Also codecov is failing because you didn't add a test.

[jasmainak]

seems like you need to rebase

[jasmainak]

I had a brief look at the automatic power line detection. That’s cool! Couple of comments:

• can we avoid making a big copy of the raw data? Computing PSD on one data channel is enough?

• is the peak detection robust? If not, you could just compare the power at 50 and 60 Hz.

Sorry I can’t look in more detail now ...

[adam2392]

I had a brief look at the automatic power line detection. That’s cool! Couple of comments:

• can we avoid making a big copy of the raw data? Computing PSD on one data channel is enough?
• is the peak detection robust? If not, you could just compare the power at 50 and 60 Hz.

I make a big copy, so that we can make use of all channels with sampling data. Assuming that there isn't severely corrupted data, then this should work robustly due to the fact that averaging many channels (with possible non-noise power at 50/60 Hz ranges) gets rid of the signal, and amplifies common noise.

I agree that probably "most" of the time, a single channel suffices, but I figured this could potentially cover essentially "all" cases, since it is a more robust estimate of what frequency is common across many electrodes (e.g. line noise). I trimmed the algorithm to only analyze the minimum of 10 seconds, or the whole dataset. I'm not sure if MNE-Python has the capability of slicing the dataset before loading?

Sorry I can’t look in more detail now ...

Enjoy your vacation :). I'm just putting these in the pipeline to get them out of my way before I forget what I wanted to do.

[adam2392]

    picks = pick_types(raw.info, meg=True, eeg=True,
                       stim=True, ecog=True, seeg=True,
                       exclude='bads')
    ch_names = [raw.ch_names[k] for k in picks]

    # first bandpass the signal to just between 45 and 65 Hz
    sfreq = raw.info['sfreq']
    l_freq = 45
    h_freq = min(sfreq // 2 - 1, 65)

    # only sample first 10 seconds, or whole dataset
    tmin = 0
    tmax = max(len(raw.times), 10 * sfreq)
    data = raw.get_data(start=tmin, stop=tmax, picks=picks, return_times=False)
    data = filter_data(data, sfreq=sfreq,
                       l_freq=l_freq, h_freq=h_freq)

    # run a multi-taper FFT
    psds, freqs = psd_array_multitaper(data, sfreq=sfreq,
                                       fmin=45, fmax=65,
                                       n_jobs=1)
    
    # find the maximum occurrence of frequency and if it is close to line
    psds = np.mean(np.abs(10 * np.log10(psds)), axis=0)
    if len(freqs) != len(psds):
        raise RuntimeError("{} != {}".format(len(freqs), len(psds)))
    possible_power_lines = [50, 60]
    max_freq = freqs[np.argmax(psds)]
    estimated_index = np.argmax([abs(max_freq - x)
                                 for x in possible_power_lines])
    powerlinefrequency = possible_power_lines[estimated_index]
    return powerlinefrequency

Could be used I suppose, but yeah you're right, probably better to just check frequencies around 50 and 60 Hz and compare average PSD for all channels.

[jasmainak]

Thanks @adam2392 !

[jasmainak]

Let's close the other two PRs as well ;-)