diff --git a/doc/source/getting_started.rst b/doc/source/getting_started.rst
index 7ed6dc4b18d..dde14f99287 100755
--- a/doc/source/getting_started.rst
+++ b/doc/source/getting_started.rst
@@ -52,26 +52,6 @@ If you get a new prompt with no error messages, you should be good to go.
 Learning Python
 ---------------
 
-If you are new to Python here are some online courses:
+If you are new to Python here is a very good place to get started:
 
-    * http://docs.python.org/
-
-    * http://learnpythonthehardway.org/book/
-
-    * http://code.google.com/intl/fr/edu/languages/google-python-class/
-
-    * http://homepage.mac.com/s_lott/books/python.html
-
-More specific tutorials on scientific computing with Python using Numpy and Scipy see:
-
-    * http://showmedo.com/videotutorials/numpy (video)
-
-    * http://www.scipy.org/NumPy_for_Matlab_Users
-
-    * http://mathesaurus.sourceforge.net/matlab-numpy.html
-
-Some discussions online are:
-
-    * http://stackoverflow.com/questions/111857/what-did-you-use-to-teach-yourself-python
-
-    * http://stackoverflow.com/questions/17988/how-to-learn-python
+    * http://scipy-lectures.github.com
diff --git a/examples/time_frequency/plot_source_label_time_frequency.py b/examples/time_frequency/plot_source_label_time_frequency.py
new file mode 100644
index 00000000000..17621d49fc9
--- /dev/null
+++ b/examples/time_frequency/plot_source_label_time_frequency.py
@@ -0,0 +1,84 @@
+"""
+=========================================================
+Compute power and phase lock in label of the source space
+=========================================================
+
+Returns time-frequency maps of induced power and phase lock
+in the source space. The inverse method is linear based on dSPM inverse operator.
+
+"""
+# Authors: Alexandre Gramfort <gramfort@nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import numpy as np
+
+import mne
+from mne import fiff
+from mne.datasets import sample
+from mne.minimum_norm import read_inverse_operator, source_induced_power
+
+###############################################################################
+# Set parameters
+data_path = sample.data_path('..')
+raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
+fname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'
+label_name = 'Aud-lh'
+fname_label = data_path + '/MEG/sample/labels/%s.label' % label_name
+
+tmin, tmax, event_id = -0.2, 0.5, 1
+
+# Setup for reading the raw data
+raw = fiff.Raw(raw_fname)
+events = mne.find_events(raw)
+inverse_operator = read_inverse_operator(fname_inv)
+
+include = []
+exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053']  # bads + 2 more
+
+# picks MEG gradiometers
+picks = fiff.pick_types(raw.info, meg=True, eeg=False, eog=True,
+                                stim=False, include=include, exclude=exclude)
+
+# Load condition 1
+event_id = 1
+epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6),
+                    preload=True)
+
+# Compute a source estimate per frequency band
+frequencies = np.arange(7, 30, 2)  # define frequencies of interest
+label = mne.read_label(fname_label)
+power, phase_lock = source_induced_power(epochs, inverse_operator, frequencies,
+                            label, baseline=(-0.1, 0), baseline_mode='percent',
+                            n_cycles=2, n_jobs=1)
+
+power = np.mean(power, axis=0)  # average over sources
+phase_lock = np.mean(phase_lock, axis=0)  # average over sources
+times = epochs.times
+
+###############################################################################
+# View time-frequency plots
+import pylab as pl
+pl.clf()
+pl.subplots_adjust(0.1, 0.08, 0.96, 0.94, 0.2, 0.43)
+pl.subplot(2, 1, 1)
+pl.imshow(20 * power, extent=[times[0], times[-1],
+                                      frequencies[0], frequencies[-1]],
+          aspect='auto', origin='lower')
+pl.xlabel('Time (s)')
+pl.ylabel('Frequency (Hz)')
+pl.title('Induced power in %s' % label_name)
+pl.colorbar()
+
+pl.subplot(2, 1, 2)
+pl.imshow(phase_lock, extent=[times[0], times[-1],
+                              frequencies[0], frequencies[-1]],
+          aspect='auto', origin='lower')
+pl.xlabel('Time (s)')
+pl.ylabel('Frequency (Hz)')
+pl.title('Phase-lock in %s' % label_name)
+pl.colorbar()
+pl.show()
diff --git a/examples/time_frequency/plot_source_space_time_frequency.py b/examples/time_frequency/plot_source_space_time_frequency.py
index b5bbcc474ef..d34b9e2471c 100644
--- a/examples/time_frequency/plot_source_space_time_frequency.py
+++ b/examples/time_frequency/plot_source_space_time_frequency.py
@@ -17,7 +17,7 @@
 import mne
 from mne import fiff
 from mne.datasets import sample
-from mne.minimum_norm import read_inverse_operator, source_induced_power
+from mne.minimum_norm import read_inverse_operator, source_band_induced_power
 
 ###############################################################################
 # Set parameters
@@ -48,8 +48,8 @@
 # Compute a source estimate per frequency band
 bands = dict(alpha=[9, 11], beta=[18, 22])
 
-stcs = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
-                            use_fft=False, n_jobs=-1)
+stcs = source_band_induced_power(epochs, inverse_operator, bands, n_cycles=2,
+                                 use_fft=False, n_jobs=1)
 
 for b, stc in stcs.iteritems():
     stc.save('induced_power_%s' % b)
diff --git a/mne/baseline.py b/mne/baseline.py
index 6cb7c11a4f0..16435d4f32d 100644
--- a/mne/baseline.py
+++ b/mne/baseline.py
@@ -29,7 +29,7 @@ def rescale(data, times, baseline, mode, verbose=True, copy=True):
         If baseline is equal ot (None, None) all the time
         interval is used.
 
-    mode: 'logratio' | 'ratio' | 'zscore' | 'mean'
+    mode: 'logratio' | 'ratio' | 'zscore' | 'mean' | 'percent'
         Do baseline correction with ratio (power is divided by mean
         power during baseline) or zscore (power is divided by standard
         deviatio of power during baseline after substracting the mean,
@@ -44,7 +44,7 @@ def rescale(data, times, baseline, mode, verbose=True, copy=True):
     if copy:
         data = data.copy()
 
-    valid_modes = ['logratio', 'ratio', 'zscore', 'mean']
+    valid_modes = ['logratio', 'ratio', 'zscore', 'mean', 'percent']
     if mode not in valid_modes:
         raise Exception('mode should be any of : %s' % valid_modes)
 
@@ -73,6 +73,9 @@ def rescale(data, times, baseline, mode, verbose=True, copy=True):
             std = np.std(data[..., imin:imax], axis=-1)[..., None]
             data -= mean
             data /= std
+        elif mode == 'percent':
+            data -= mean
+            data /= mean
 
     elif verbose:
         print "No baseline correction applied..."
diff --git a/mne/minimum_norm/__init__.py b/mne/minimum_norm/__init__.py
index 6a22998bfe9..3ee45ec01ce 100644
--- a/mne/minimum_norm/__init__.py
+++ b/mne/minimum_norm/__init__.py
@@ -1,4 +1,4 @@
 from .inverse import read_inverse_operator, apply_inverse, \
                      apply_inverse_raw, make_inverse_operator, \
                      apply_inverse_epochs
-from .time_frequency import source_induced_power
+from .time_frequency import source_band_induced_power, source_induced_power
diff --git a/mne/minimum_norm/tests/test_time_frequency.py b/mne/minimum_norm/tests/test_time_frequency.py
index b5fd78fc077..c6a1d9c5ce4 100644
--- a/mne/minimum_norm/tests/test_time_frequency.py
+++ b/mne/minimum_norm/tests/test_time_frequency.py
@@ -1,13 +1,14 @@
 import os.path as op
 
 import numpy as np
-from numpy.testing import assert_array_almost_equal, assert_equal
+from numpy.testing import assert_array_almost_equal
 from nose.tools import assert_true
 
 from ...datasets import sample
 from ... import fiff, find_events, Epochs
+from ...label import read_label
 from ..inverse import read_inverse_operator
-from ..time_frequency import source_induced_power
+from ..time_frequency import source_band_induced_power, source_induced_power
 
 
 examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples')
@@ -16,10 +17,11 @@
                                         'sample_audvis-meg-oct-6-meg-inv.fif')
 fname_data = op.join(data_path, 'MEG', 'sample',
                                         'sample_audvis_raw.fif')
+fname_label = op.join(data_path, 'MEG', 'sample', 'labels', 'Aud-lh.label')
 
 
 def test_tfr_with_inverse_operator():
-    """Test MNE inverse computation"""
+    """Test time freq with MNE inverse computation"""
 
     tmin, tmax, event_id = -0.2, 0.5, 1
 
@@ -33,27 +35,41 @@ def test_tfr_with_inverse_operator():
 
     # picks MEG gradiometers
     picks = fiff.pick_types(raw.info, meg=True, eeg=False, eog=True,
-                                    stim=False, include=include, exclude=exclude)
+                                stim=False, include=include, exclude=exclude)
 
     # Load condition 1
     event_id = 1
-    events = events[:3]  # take 3 events to keep the computation time low
-    epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
-                        baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6),
-                        preload=True)
+    events3 = events[:3]  # take 3 events to keep the computation time low
+    epochs = Epochs(raw, events3, event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6),
+                    preload=True)
 
     # Compute a source estimate per frequency band
     bands = dict(alpha=[10, 10])
+    label = read_label(fname_label)
 
-    stcs = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
-                                use_fft=False, pca=True)
+    stcs = source_band_induced_power(epochs, inverse_operator, bands,
+                            n_cycles=2, use_fft=False, pca=True, label=label)
 
     stc = stcs['alpha']
     assert_true(len(stcs) == len(bands.keys()))
     assert_true(np.all(stc.data > 0))
     assert_array_almost_equal(stc.times, epochs.times)
 
-    stcs_no_pca = source_induced_power(epochs, inverse_operator, bands, n_cycles=2,
-                                use_fft=False, pca=False)
+    stcs_no_pca = source_band_induced_power(epochs, inverse_operator, bands,
+                            n_cycles=2, use_fft=False, pca=False, label=label)
 
     assert_array_almost_equal(stcs['alpha'].data, stcs_no_pca['alpha'].data)
+
+    # Compute a source estimate per frequency band
+    epochs = Epochs(raw, events[:10], event_id, tmin, tmax, picks=picks,
+                    baseline=(None, 0), reject=dict(grad=4000e-13, eog=150e-6),
+                    preload=True)
+
+    frequencies = np.arange(7, 30, 2)  # define frequencies of interest
+    power, phase_lock = source_induced_power(epochs, inverse_operator,
+                            frequencies, label, baseline=(-0.1, 0),
+                            baseline_mode='percent', n_cycles=2, n_jobs=1)
+    assert_true(np.all(phase_lock > 0))
+    assert_true(np.all(phase_lock < 1))
+    assert_true(np.max(power) > 10)
diff --git a/mne/minimum_norm/time_frequency.py b/mne/minimum_norm/time_frequency.py
index b36adce5969..b596046c26c 100644
--- a/mne/minimum_norm/time_frequency.py
+++ b/mne/minimum_norm/time_frequency.py
@@ -6,47 +6,19 @@
 from scipy import linalg
 
 from ..fiff.constants import FIFF
-from ..source_estimate import SourceEstimate
 from ..time_frequency.tfr import cwt, morlet
 from ..baseline import rescale
-from .inverse import combine_xyz, prepare_inverse_operator
+from .inverse import combine_xyz, prepare_inverse_operator, _assemble_kernel, \
+                     _make_stc
 from ..parallel import parallel_func
 
 
-def _compute_power(data, K, sel, Ws, source_ori, use_fft, Vh):
-    """Aux function for source_induced_power"""
-    power = 0
-
-    for e in data:
-        e = e[sel]  # keep only selected channels
-
-        if Vh is not None:
-            e = np.dot(Vh, e)  # reducing data rank
-
-        for w in Ws:
-            tfr = cwt(e, [w], use_fft=use_fft)
-            tfr = np.asfortranarray(tfr.reshape(len(e), -1))
-
-            for t in [np.real(tfr), np.imag(tfr)]:
-                sol = np.dot(K, t)
-
-                if source_ori == FIFF.FIFFV_MNE_FREE_ORI:
-                    # print 'combining the current components...',
-                    sol = combine_xyz(sol, square=True)
-                else:
-                    np.power(sol, 2, sol)
-
-                power += sol
-                del sol
-
-    return power
-
-
-def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
-                         dSPM=True, n_cycles=5, df=1, use_fft=False,
-                         baseline=None, baseline_mode='logratio', pca=True,
-                         subtract_evoked=False, n_jobs=1):
-    """Compute source space induced power
+def source_band_induced_power(epochs, inverse_operator, bands, label=None,
+                              lambda2=1.0 / 9.0, dSPM=True, n_cycles=5, df=1,
+                              use_fft=False, baseline=None,
+                              baseline_mode='logratio', pca=True,
+                              n_jobs=1):
+    """Compute source space induced power in given frequency bands
 
     Parameters
     ----------
@@ -56,6 +28,8 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
         The inverse operator
     bands: dict
         Example : bands = dict(alpha=[8, 9])
+    label: Label
+        Restricts the source estimates to a given label
     lambda2: float
         The regularization parameter of the minimum norm
     dSPM: bool
@@ -83,23 +57,120 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
         If True, the true dimension of data is estimated before running
         the time frequency transforms. It reduces the computation times
         e.g. with a dataset that was maxfiltered (true dim is 64)
-    subtract_evoked: bool
-        If True, the evoked component (average of all epochs) if subtracted
-        from each epochs.
     n_jobs: int
         Number of jobs to run in parallel
     """
+    frequencies = np.concatenate([np.arange(band[0], band[1] + df / 2.0, df)
+                                 for _, band in bands.iteritems()])
+
+    powers, _, lh_vertno, rh_vertno = _source_induced_power(epochs,
+                                      inverse_operator, frequencies,
+                                      label=label,
+                                      lambda2=lambda2, dSPM=dSPM,
+                                      n_cycles=n_cycles,
+                                      use_fft=use_fft, pca=pca, n_jobs=n_jobs,
+                                      with_plv=False)
+
+    Fs = epochs.info['sfreq']  # sampling in Hz
+    stcs = dict()
+
+    for name, band in bands.iteritems():
+        idx = [k for k, f in enumerate(frequencies) if band[0] <= f <= band[1]]
+
+        # average power in band + mean over epochs
+        power = np.mean(powers[:, idx, :], axis=1)
 
-    parallel, my_compute_power, n_jobs = parallel_func(_compute_power, n_jobs)
+        # Run baseline correction
+        power = rescale(power, epochs.times, baseline, baseline_mode,
+                        verbose=True, copy=False)
+
+        stc = _make_stc(power, epochs.times[0], 1.0 / Fs, lh_vertno, rh_vertno)
+        stcs[name] = stc
+
+        print '[done]'
+
+    return stcs
+
+
+def _compute_pow_plv(data, K, sel, Ws, source_ori, use_fft, Vh, with_plv,
+                     pick_normal):
+    """Aux function for source_induced_power"""
+    n_times = data.shape[2]
+    n_freqs = len(Ws)
+    n_sources = K.shape[0]
+    is_free_ori = False
+    if source_ori == FIFF.FIFFV_MNE_FREE_ORI:
+        is_free_ori = True
+        n_sources /= 3
+
+    shape = (n_sources, n_freqs, n_times)
+    power = np.zeros(shape, dtype=np.float)  # power
+    if with_plv:
+        shape = (n_sources, n_freqs, n_times)
+        plv = np.zeros(shape, dtype=np.complex)  # phase lock
+    else:
+        plv = None
+
+    for e in data:
+        e = e[sel]  # keep only selected channels
+
+        if Vh is not None:
+            e = np.dot(Vh, e)  # reducing data rank
+
+        for f, w in enumerate(Ws):
+            tfr = cwt(e, [w], use_fft=use_fft)
+            tfr = np.asfortranarray(tfr.reshape(len(e), -1))
+
+            # phase lock and power at freq f
+            if with_plv:
+                plv_f = np.zeros((n_sources, n_times), dtype=np.complex)
+            pow_f = np.zeros((n_sources, n_times), dtype=np.float)
 
+            for k, t in enumerate([np.real(tfr), np.imag(tfr)]):
+                sol = np.dot(K, t)
+
+                sol_pick_normal = sol
+                if is_free_ori:
+                    sol_pick_normal = sol[2::3]
+
+                if with_plv:
+                    if k == 0:  # real
+                        plv_f += sol_pick_normal
+                    else:  # imag
+                        plv_f += 1j * sol_pick_normal
+
+                if is_free_ori:
+                    # print 'combining the current components...',
+                    sol = combine_xyz(sol, square=True)
+                else:
+                    np.power(sol, 2, sol)
+                pow_f += sol
+                del sol
+
+            power[:, f, :] += pow_f
+            del pow_f
+
+            if with_plv:
+                plv_f /= np.abs(plv_f)
+                plv[:, f, :] += plv_f
+                del plv_f
+
+    return power, plv
+
+
+def _source_induced_power(epochs, inverse_operator, frequencies, label=None,
+                          lambda2=1.0 / 9.0, dSPM=True, n_cycles=5,
+                          use_fft=False, pca=True, pick_normal=True,
+                          n_jobs=1, with_plv=True):
+    """Aux function for source_induced_power
+    """
+    parallel, my_compute_pow_plv, n_jobs = parallel_func(_compute_pow_plv,
+                                                         n_jobs)
     #
     #   Set up the inverse according to the parameters
     #
     epochs_data = epochs.get_data()
 
-    if subtract_evoked:  # subtract with a copy not to touch epochs
-        epochs_data = epochs_data - np.mean(epochs_data, axis=0)
-
     nave = len(epochs_data)  # XXX : can do better when no preload
 
     inv = prepare_inverse_operator(inverse_operator, nave, lambda2, dSPM)
@@ -116,10 +187,7 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
     #   This does all the data transformations to compute the weights for the
     #   eigenleads
     #
-    K = inv['reginv'][:, None] * reduce(np.dot,
-                                           [inv['eigen_fields']['data'],
-                                           inv['whitener'],
-                                           inv['proj']])
+    K, noise_norm, lh_vertno, rh_vertno = _assemble_kernel(inv, label, dSPM)
 
     if pca:
         U, s, Vh = linalg.svd(K)
@@ -130,61 +198,92 @@ def source_induced_power(epochs, inverse_operator, bands, lambda2=1.0 / 9.0,
     else:
         Vh = None
 
-    #
-    #   Transformation into current distributions by weighting the
-    #   eigenleads with the weights computed above
-    #
-    if inv['eigen_leads_weighted']:
-        #
-        #     R^0.5 has been already factored in
-        #
-        # print '(eigenleads already weighted)...',
-        K = np.dot(inv['eigen_leads']['data'], K)
-    else:
-        #
-        #     R^0.5 has to factored in
-        #
-        # print '(eigenleads need to be weighted)...',
-        K = np.sqrt(inv['source_cov']['data'])[:, None] * \
-                             np.dot(inv['eigen_leads']['data'], K)
-
     Fs = epochs.info['sfreq']  # sampling in Hz
 
-    stcs = dict()
-    src = inv['src']
+    print 'Computing source power ...'
 
-    for name, band in bands.iteritems():
-        print 'Computing power in band %s [%s, %s] Hz...' % (name, band[0],
-                                                             band[1])
+    Ws = morlet(Fs, frequencies, n_cycles=n_cycles)
 
-        freqs = np.arange(band[0], band[1] + df / 2.0, df)  # frequencies
-        Ws = morlet(Fs, freqs, n_cycles=n_cycles)
+    n_jobs = min(n_jobs, len(epochs_data))
+    out = parallel(my_compute_pow_plv(data, K, sel, Ws,
+                                      inv['source_ori'], use_fft, Vh,
+                                      with_plv, pick_normal)
+                        for data in np.array_split(epochs_data, n_jobs))
+    power = sum(o[0] for o in out)
+    power /= len(epochs_data)  # average power over epochs
 
-        power = sum(parallel(my_compute_power(data, K, sel, Ws,
-                                            inv['source_ori'], use_fft, Vh)
-                            for data in np.array_split(epochs_data, n_jobs)))
+    if with_plv:
+        plv = sum(o[1] for o in out)
+        plv = np.abs(plv)
+        plv /= len(epochs_data)  # average power over epochs
+    else:
+        plv = None
 
-        if dSPM:
-            # print '(dSPM)...',
-            power *= inv['noisenorm'][:, None] ** 2
+    if dSPM:
+        power *= noise_norm.ravel()[:, None, None] ** 2
 
-        # average power in band + mean over epochs
-        power /= len(epochs_data) * len(freqs)
+    return power, plv, lh_vertno, rh_vertno
 
-        # Run baseline correction
-        power = rescale(power, epochs.times, baseline, baseline_mode,
-                        verbose=True, copy=False)
 
-        stc = SourceEstimate(None)
-        stc.data = power
-        stc.tmin = epochs.times[0]
-        stc.tstep = 1.0 / Fs
-        stc.lh_vertno = src[0]['vertno']
-        stc.rh_vertno = src[1]['vertno']
-        stc._init_times()
+def source_induced_power(epochs, inverse_operator, frequencies, label=None,
+                         lambda2=1.0 / 9.0, dSPM=True, n_cycles=5,
+                         use_fft=False, pick_normal=False, baseline=None,
+                         baseline_mode='logratio', pca=True, n_jobs=1):
+    """Compute induced power and phase lock
 
-        stcs[name] = stc
+    Computation can optionaly be restricted in a label.
 
-        print '[done]'
+    Parameters
+    ----------
+    epochs: instance of Epochs
+        The epochs
+    inverse_operator: instance of inverse operator
+        The inverse operator
+    label: Label
+        Restricts the source estimates to a given label
+    frequencies: array
+        Array of frequencies of interest
+    lambda2: float
+        The regularization parameter of the minimum norm
+    dSPM: bool
+        Do dSPM or not?
+    n_cycles: int
+        Number of cycles
+    use_fft: bool
+        Do convolutions in time or frequency domain with FFT
+    pick_normal: bool
+        If True, rather than pooling the orientations by taking the norm,
+        only the radial component is kept. This is only implemented
+        when working with loose orientations.
+    baseline: None (default) or tuple of length 2
+        The time interval to apply baseline correction.
+        If None do not apply it. If baseline is (a, b)
+        the interval is between "a (s)" and "b (s)".
+        If a is None the beginning of the data is used
+        and if b is None then b is set to the end of the interval.
+        If baseline is equal ot (None, None) all the time
+        interval is used.
+    baseline_mode: None | 'logratio' | 'zscore'
+        Do baseline correction with ratio (power is divided by mean
+        power during baseline) or zscore (power is divided by standard
+        deviatio of power during baseline after substracting the mean,
+        power = [power - mean(power_baseline)] / std(power_baseline))
+    pca: bool
+        If True, the true dimension of data is estimated before running
+        the time frequency transforms. It reduces the computation times
+        e.g. with a dataset that was maxfiltered (true dim is 64)
+    n_jobs: int
+        Number of jobs to run in parallel
+    """
+    power, plv, lh_vertno, rh_vertno = _source_induced_power(epochs,
+                            inverse_operator, frequencies,
+                            label, lambda2, dSPM, n_cycles,
+                            use_fft, pick_normal=pick_normal, pca=pca,
+                            n_jobs=n_jobs)
+
+    # Run baseline correction
+    if baseline is not None:
+        power = rescale(power, epochs.times, baseline, baseline_mode,
+                        verbose=True, copy=False)
 
-    return stcs
+    return power, plv