Test of aliasing in receptive field estimation

rf_aliasing.py

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Fri Dec 30 16:03:52 2016

@author: larsoner
"""

import numpy as np
from mne.decoding import ReceptiveField
import matplotlib.pyplot as plt

sfreq = 100.
n_signal = int(100 * sfreq)
n_system = int(sfreq)
X = np.random.RandomState(0).randn(n_signal, 1)
freqs = [50., 49.]  # alias to DC and 1 Hz
h_x = np.arange(0, -n_system, -1) / sfreq
fig, axes = plt.subplots(2, 2)
kinds = ['No', 'Input']
tmin, tmax = -2., 0.1
for fi, freq in enumerate(freqs):
    h = np.cos(2 * np.pi * freq * np.arange(n_system) / sfreq)
    y = np.convolve(X[:, 0], h)[:n_signal]

    for ki, kind in enumerate(kinds):
        if ki == 0:
            rf = ReceptiveField(tmin, tmax, sfreq)
            rf.fit(X, y)
        else:
            rf = ReceptiveField(tmin, tmax, sfreq / 2.)
            rf.fit(X[::2], y[::2])
        axes[ki, fi].plot(h_x, h, linewidth=2)
        axes[ki, fi].plot(np.arange(rf.coef_.shape[0]) / rf.sfreq + rf.tmin,
                          rf.coef_)
        if fi == 0:
            axes[ki, fi].set(ylabel='%s\ndecimation' % kind)
        if ki == 0:
            axes[0, fi].set(title='%0.0f Hz' % freq)
fig.tight_layout()