Merge pull request #4 from larsoner/csd

ENH: Restore default reg

mne/simulation/raw.py

@@ -185,7 +185,8 @@ def simulate_raw(raw, stc, trans, src, bem, cov='simple',
     if not isinstance(stc, _BaseSourceEstimate):
         raise TypeError('stc must be a SourceEstimate')
     if not np.allclose(info['sfreq'], 1. / stc.tstep):
-        raise ValueError('stc and info must have same sample rate')
+        raise ValueError('stc (%s) and info (%s) must have same sample rate'
+                         % (1. / stc.tstep, info['sfreq']))
     if len(stc.times) <= 2:  # to ensure event encoding works
         raise ValueError('stc must have at least three time points')
 

tutorials/plot_dics.py

@@ -33,9 +33,6 @@
 from mne.time_frequency import csd_epochs
 from mne.beamformer import make_dics, apply_dics_csd
 
-# Suppress irrelevant output
-mne.set_log_level('ERROR')
-
 # We use the MEG and MRI setup from the MNE-sample dataset
 data_path = sample.data_path(download=False)
 subjects_dir = op.join(data_path, 'subjects')
@@ -45,9 +42,9 @@
 meg_path = op.join(data_path, 'MEG', 'sample')
 raw_fname = op.join(meg_path, 'sample_audvis_raw.fif')
 trans_fname = op.join(meg_path, 'sample_audvis_raw-trans.fif')
-src_fname = op.join(mri_path, 'bem/sample-oct-6-src.fif')
-bem_fname = op.join(mri_path, 'bem/sample-5120-5120-5120-bem-sol.fif')
-fwd_fname = op.join(meg_path, 'sample_audvis-meg-eeg-oct-6-fwd.fif')
+src_fname = op.join(mri_path, 'bem', 'sample-oct-6-src.fif')
+bem_fname = op.join(mri_path, 'bem', 'sample-5120-bem-sol.fif')
+fwd_fname = op.join(meg_path, 'sample_audvis-meg-oct-6-fwd.fif')
 cov_fname = op.join(meg_path, 'sample_audvis-cov.fif')
 
 # Seed for the random number generator
@@ -59,10 +56,15 @@
 #
 # The following function generates a timeseries that contains an oscillator,
 # whose frequency fluctuates a little over time, but stays close to 10 Hz.
-# We'll use this function to generate our two signals.
+# We'll use this function to generate our two signals. Each signal will be
+# zero for negative time, and have an oscillator in positive time.
 
-sfreq = 50.  # Sampling frequency of the generated signal
-times = np.arange(10. * sfreq) / sfreq  # 10 seconds of signal
+sfreq = 100.  # Sampling frequency of the generated signal
+base_freq = 10.  # Base frequency of the oscillators in Hertz
+n_trials = 5
+tmin = -0.5  # amount of zero signal to produce (sec)
+tmax = 2.  # amount of oscillating signal to produce (sec)
+times = np.arange(int(round(tmax * sfreq)) + 1) / sfreq
 
 
 def coh_signal_gen():
@@ -73,23 +75,21 @@ def coh_signal_gen():
     signal : ndarray
         The generated signal.
     """
-    t_rand = 0.001  # Variation in the instantaneous frequency of the signal
+    t_rand = 1e-1 / sfreq  # Variation in the instantaneous frequency / phase
     std = 0.1  # Std-dev of the random fluctuations added to the signal
-    base_freq = 10.  # Base frequency of the oscillators in Hertz
     n_times = len(times)
 
     # Generate an oscillator with varying frequency and phase lag.
-    iflaw = base_freq / sfreq + t_rand * rand.randn(n_times)
-    signal = np.exp(1j * 2.0 * np.pi * np.cumsum(iflaw))
-    signal *= np.conj(signal[0])
-    signal = signal.real
+    signal = np.sin(2.0 * np.pi *
+                    (base_freq * np.arange(n_times) / sfreq +
+                     np.cumsum(t_rand * rand.randn(n_times))))
 
     # Add some random fluctuations to the signal.
     signal += std * rand.randn(n_times)
 
-    # Scale the signal to be in the right order of magnitude (~500 nAm)
+    # Scale the signal to be in the right order of magnitude (~250 nAm)
     # to achieve a SNR of 1 with our noise covariance matrix.
-    signal *= 500e-9
+    signal *= 250e-9
 
     return signal
 
@@ -111,18 +111,18 @@ def coh_signal_gen():
 ax.set(xlabel='Time (s)', xlim=times[[0, -1]], title='Signal 2')
 
 # Power spectrum of the first timeseries
-f, p = welch(signal1, fs=sfreq, nperseg=128, nfft=256)
+freq_kwargs = dict(fs=sfreq, nperseg=50, noverlap=25, detrend=False)
+f, p = welch(signal1, **freq_kwargs)
 ax = axes[1][0]
-# Only plot the first 100 frequencies
-ax.plot(f[:100], 20 * np.log10(p[:100]), lw=1.)
-ax.set(xlabel='Frequency (Hz)', xlim=f[[0, 99]],
+ax.plot(f, 20 * np.log10(p), lw=1.)
+ax.set(xlabel='Frequency (Hz)', xlim=f[[0, -1]],
        ylabel='Power (dB)', title='Power spectrum of signal 1')
 
 # Compute the coherence between the two timeseries
-f, coh = coherence(signal1, signal2, fs=sfreq, nperseg=100, noverlap=64)
+f, coh = coherence(signal1, signal2, **freq_kwargs)
 ax = axes[1][1]
-ax.plot(f[:50], coh[:50], lw=1.)
-ax.set(xlabel='Frequency (Hz)', xlim=f[[0, 49]], ylabel='Coherence',
+ax.plot(f, coh, lw=1.)
+ax.set(xlabel='Frequency (Hz)', xlim=f[[0, -1]], ylabel='Coherence',
        title='Coherence between the timeseries')
 fig.tight_layout()
 
@@ -140,18 +140,18 @@ def coh_signal_gen():
 source_vert1 = 146374
 source_vert2 = 33830
 
-# The timeseries at each vertex: one part signal, one part silence
-timeseries1 = np.hstack([signal1, np.zeros_like(signal1)])
-timeseries2 = np.hstack([signal2, np.zeros_like(signal2)])
+# The timeseries at each vertex: one part zeros, one part signal
+n_noise = int(round(-tmin * sfreq))
+timeseries1 = np.hstack([np.zeros(n_noise), signal1])
+timeseries2 = np.hstack([np.zeros(n_noise), signal2])
 
 # Construct a SourceEstimate object that describes the signal at the cortical
 # level.
+tstep = 1. / sfreq
 stc = mne.SourceEstimate(
     np.vstack((timeseries1, timeseries2)),  # The two timeseries
     vertices=[[source_vert1], [source_vert2]],  # Their locations
-    tmin=0,
-    tstep=1. / sfreq,
-    subject='sample',  # We use the brain model of the MNE-Sample dataset
+    tmin=tmin, tstep=tstep, subject='sample',
 )
 
 ###############################################################################
@@ -172,31 +172,32 @@ def coh_signal_gen():
 cov['data'] /= snr ** 2
 
 # This is the raw file we're going to use as template for the simulated data.
-raw = mne.io.read_raw_fif(raw_fname, preload=True)
-raw = raw.crop(0, 20).resample(sfreq)  # Trim to 20 seconds at 50 Hz.
-raw = raw.pick_types(meg='grad')  # Use only gradiometers
-
-# Simulate the raw data
+# Here we will restrict to gradiometers for speed.
+info = mne.io.read_info(raw_fname)
+info.update(sfreq=sfreq, bads=[])
+picks = mne.pick_types(info, meg='grad', stim=True, exclude=())
+mne.pick_info(info, picks, copy=False)
+raw = mne.io.RawArray(np.zeros((len(info['ch_names']),
+                                stc.data.shape[1] * n_trials)), info)
+
+# Simulate the raw data, with a lowpass filter on the noise
 raw = simulate_raw(raw, stc, trans_fname, src_fname, bem_fname, cov=cov,
-                   random_state=rand)
+                   random_state=rand, iir_filter=[4, -4, 0.8])
 
 ###############################################################################
-# We create an :class:`mne.Epochs` object containing two trials: one with
-# both noise and signal and one with just noise
-
-t0 = raw.first_samp  # First sample int the data
-t10 = t0 + int(10 * sfreq) - 1  # Sample just before the 10 second mark
-epochs = mne.Epochs(
-    raw,
-    events=np.array([[t0, 0, 1], [t10, 0, 2]]),
-    event_id=dict(signal=1, noise=2),
-    tmin=0, tmax=10,
-    preload=True,
-)
+# We create an :class:`mne.Epochs` object containing trials where positive
+# time contains both noise and signal, and negative time only noise.
+
+events = mne.find_events(raw, 'STI 014')
+assert len(events) == n_trials
+epochs = mne.Epochs(raw, events, tmin=tmin, tmax=tmax, preload=True)
 
 # Plot the simulated data
 epochs.plot()
 
+# crop to just our region with signal
+epochs.crop(0., None)
+
 ###############################################################################
 # Power mapping
 # -------------
@@ -213,53 +214,49 @@ def coh_signal_gen():
 
 # Estimating the noise covariance on the trial that only contains noise.
 fwd = mne.read_forward_solution(fwd_fname)
-inv = make_inverse_operator(epochs.info, fwd, cov)
+inv = make_inverse_operator(info, fwd, cov)
 
 # Apply the inverse model to the trial that also contains the signal.
-s = apply_inverse(epochs['signal'].average(), inv)
+s = apply_inverse(epochs.average(), inv)
 
 # Take the root-mean square along the time dimension and plot the result.
-s_rms = (s ** 2).mean()
+s_rms = np.sqrt((s ** 2).mean())
 brain = s_rms.plot('sample', subjects_dir=subjects_dir, hemi='both', figure=1,
-                   size=400)
-
-# Plot the result
-brain = s.plot('sample', subjects_dir=subjects_dir, hemi='both', figure=1,
-               size=400)
+                   size=600)
 
 # Indicate the true locations of the source activity on the plot.
 brain.add_foci(source_vert1, coords_as_verts=True, hemi='lh')
 brain.add_foci(source_vert2, coords_as_verts=True, hemi='rh')
 
 # Rotate the view and add a title.
 mlab.view(0, 0, 550, [0, 0, 0])
-mlab.title('MNE-dSPM inverse (RMS)', height=0.9)
+title = mlab.title('dSPM inverse\n(RMS)', height=0.9)
 
 ###############################################################################
-# Computing a cortical power map at 10 Hz. using a DICS beamformer:
+# Computing a cortical power map at the base frequency using a DICS beamformer:
 
 # Estimate the cross-spectral density (CSD) matrix on the trial containing the
 # signal.
-csd_signal = csd_epochs(epochs['signal'], mode='cwt_morlet', frequencies=[10])
+csd_signal = csd_epochs(epochs, mode='cwt_morlet',
+                        frequencies=[base_freq])
 
 # Compute the DICS powermap. An important parameter for this is the
 # regularization, which is set quite high for this toy example. For real data,
 # you may want to lower this to around 0.05.
-filters = make_dics(epochs.info, fwd, csd_signal, reg=0.5,
-                    pick_ori='max-power')
+filters = make_dics(info, fwd, csd_signal, reg=0.05, pick_ori='max-power')
 power, f = apply_dics_csd(csd_signal, filters)
 
 # Plot the DICS power map.
 brain = power.plot('sample', subjects_dir=subjects_dir, hemi='both', figure=2,
-                   size=400)
+                   size=600)
 
 # Indicate the true locations of the source activity on the plot.
 brain.add_foci(source_vert1, coords_as_verts=True, hemi='lh')
 brain.add_foci(source_vert2, coords_as_verts=True, hemi='rh')
 
 # Rotate the view and add a title.
 mlab.view(0, 0, 550, [0, 0, 0])
-mlab.title('DICS power map at %.1f Hz' % f[0], height=0.9)
+title = mlab.title('DICS power map\n%.1f Hz' % f[0], height=0.9)
 
 ###############################################################################
 # Excellent! Both methods found our two simulated sources. Of course, with a