Я новичок в scikit-learn, но я пытаюсь удалить моргание глаз (шумовые пики) внутри одного канала ЭЭГ. Я искал в Интернете, но вижу только более сложные показания с MNE, PyEEG или другими модулями Python. Я просто хочу что-то простое и зависящее только от sklearn. Вот мой код:
#The channel containing some eye-blinks
X = f1ep1_data[:,[4]]
#run ICA on signal
ica = FastICA(n_components=2)
ica.fit(X)
#reconstruct signal with independent components
components = ica.fit_transform(X)
X_restored = ica.inverse_transform(components)
fig1 = plt.figure()
plt.subplot(3,1,1)
plt.title("Original signal")
plt.plot(f1ep1_timescale, X)
plt.subplot(3,1,2)
plt.title("Components")
plt.plot(f1ep1_timescale, components)
plt.subplot(3,1,3)
plt.title("Signal Reconstructed")
plt.plot(f1ep1_timescale, X_restored)
plt.draw()
Я ожидал разделения на две составляющие: более чистый сигнал ЭЭГ и моргание глаз. Я не могу понять, в чем проблема. Может кто-нибудь протянуть руку помощи?
Вы заметили, что ваши «компоненты» - это точно масштабированный исходный сигнал, перевернутый вверх ногами? Это потому, что вы не можете получить больше компонентов, чем сигналов.
Вам необходимо выполнить следующие действия:
Давайте подробно рассмотрим шаг 2: зачем удалять компоненты вручную? ICA ничего не знает о моргании глаз. Он разделяет сигналы на компоненты на основе статистических показателей. Если вам повезет, некоторые из этих компонентов будут похожи на моргание глаз.
Пока что хорошо, но настоящая проблема в том, что порядок компонентов не определен. Запустите ICA, и вы можете обнаружить, что компонент 1 моргает. Запустите его снова, и они будут в компоненте 3. Снова, и они будут в обоих компонентах 2 и 5 ...
Нет возможности заранее узнать, какие и сколько компонентов нужно удалить. Вот почему вам нужно вручную сообщать об этом алгоритму при каждом его запуске.
В коде это выглядело бы примерно так:
# Use all channels - they will contain eye blinks to varying degrees
X = f1ep1_data[:, :]
# run ICA on signal
ica = FastICA(n_components=x.shape[1]) # we want *all* the components
ica.fit(X)
# decompose signal into components
components = ica.fit_transform(X)
# plot components and ask user which components to remove
# ...
remove_indices = [0, 1, 3] # pretend the user selected components 0, 1, and 3
# "remove" unwanted components by setting them to 0 - simplistic but gets the job done
components[:, remove_indices] = 0
#reconstruct signal
X_restored = ica.inverse_transform(components)
Скорее всего, вас не устраивает ручной шаг. Не повезло :) В 2013 году не было надежного автоматического алгоритма, который бы отмечал компоненты моргания. Я не думаю, что это изменилось, но если есть что-то, вы найдете это в одном из пакетов, специфичных для домена, например MNE или PyEEG.
Однако, если у вас есть записи EOG, есть надежда! Существует Полностью автоматизированный метод коррекции артефактов ЭОГ в записях ЭЭГ.. Этот подход основан на канонической корреляции или регрессии (я не помню подробностей), но вам необходимо записывать сигналы ЭОГ вместе с ЭЭГ.
Я создал рабочий пример с смоделированными данными «ЭЭГ». Данные состоят из трех каналов: лобного, центрального и теменного. Альфа-активность 10 Гц наиболее сильна сзади, а несколько импульсов, подобных миганию, наиболее сильны спереди.
Надеюсь, этот пример лучше иллюстрирует, как удалить компоненты из многоканальных данных.
import numpy as np
import scipy.signal as sps
from sklearn.decomposition import FastICA
import matplotlib.pyplot as plt
np.random.seed(42)
n = 1000
fs = 100
noise = 3
# simulate EEG data with eye blinks
t = np.arange(n)
alpha = np.abs(np.sin(10 * t / fs)) - 0.5
alpha[n//2:] = 0
blink = np.zeros(n)
blink[n//2::200] += -1
blink = sps.lfilter(*sps.butter(2, [1*2/fs, 10*2/fs], 'bandpass'), blink)
frontal = blink * 200 + alpha * 10 + np.random.randn(n) * noise
central = blink * 100 + alpha * 15 + np.random.randn(n) * noise
parietal = blink * 10 + alpha * 25 + np.random.randn(n) * noise
eeg = np.stack([frontal, central, parietal]).T # shape = (100, 3)
# plot original data
plt.subplot(3, 1, 1)
plt.plot(frontal + 50)
plt.plot(central + 100)
plt.plot(parietal + 150)
plt.yticks([50, 100, 150], ['Fz', 'Cz', 'Pz'])
plt.ylabel('original data')
# decompose EEG and plot components
ica = FastICA(n_components=3)
ica.fit(eeg)
components = ica.transform(eeg)
plt.subplot(3, 1, 2)
plt.plot([[np.nan, np.nan, np.nan]]) # advance the color cycler to give the components a different color :)
plt.plot(components + [0.5, 1.0, 1.5])
plt.yticks([0.5, 1.0, 1.5], ['0', '1', '2'])
plt.ylabel('components')
# looks like component #2 (brown) contains the eye blinks
# let's remove them (hard coded)!
components[:, 2] = 0
# reconstruct EEG without blinks
restored = ica.inverse_transform(components)
plt.subplot(3, 1, 3)
plt.plot(restored + [50, 100, 150])
plt.yticks([50, 100, 150], ['Fz', 'Cz', 'Pz'])
plt.ylabel('cleaned data')
Даже вручную выбирая моргание глаз, я получаю только (151L,1L) float64
от ica_data = ica.fit_transform(blink)
. Полагаю, я не использую функцию ICA, как предполагалось ...
ICA преобразует каналы N в самое большее количество компонентов N. если * N * = 1 канал -> 1 компонент; это не очень полезно.
@Pedro Я обновил ответ конкретным примером, с которым вы можете поиграть. Это поможет продемонстрировать, как работает декомпозиция ICA.
Это фантастика @kazemakase! Большое спасибо за обновленный ответ и правило N-N ICA. Я смог сделать это с 4 каналами. Еще пара вопросов: 1)
каждый раз, когда запускается код, компонент eyeblinks может быть другим. Как узнать, что это? Может свертка с блинк-шаблоном? Что ты предлагаешь? 2)
, что было бы лучшей альтернативой для жесткого кодирования компонента eyeblinks?
Вижу, у меня еще много работы ... Вы ведь эксперт по ЭЭГ? Как вы смоделировали форму волны моргания (не использовали фильтр Баттуорта)? Эта форма волны может быть использована для соответствия сигналу.
Я больше не работаю в этой сфере. Что касается фильтра, я использовал Баттерворта, потому что это была первая функция, которая пришла в голову. Все, что мне было нужно, - это полоса пропускания, которую можно было применять на импульсных всплесках. Затем я поигрался с частотами фильтра и порядком, пока все не выглядело как-то правильно;)
В случае, если вы имеете дело с одноканальной ЭЭГ, простой в реализации метод может быть следующим:
1) Обнаруживайте мигания в сигнале Икс, используя простое определение пиков на основе пороговых значений (возможно, вам придется выяснить это, посмотрев на свои сигналы для нескольких случаев мигания). Устройства от Neurosky, Muse и т. д. Поставляются с API для обнаружения моргания. Вы можете использовать это при необходимости.
2) Возьмите кадр, соответствующий миганию (xb). Нарисуйте на нем плавную линию (xbs).
3) Вычтите плавную линию (xbs) из мигания (xb) и добавьте к ней среднее значение сигнала. Назовем это xbc.
4) Заменить xbc вместо xb в исходном сигнале Икс.
Почему мне нужно передавать каналы все в ICA? Что делать, если я использую только один канал? ICA работать не будет?