Устранение моргания глаз из сигнала ЭЭГ с помощью ICA

Я новичок в 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()

Устранение моргания глаз из сигнала ЭЭГ с помощью ICA

Я ожидал разделения на две составляющие: более чистый сигнал ЭЭГ и моргание глаз. Я не могу понять, в чем проблема. Может кто-нибудь протянуть руку помощи?

Почему в Python есть оператор "pass"?
Почему в Python есть оператор "pass"?
Оператор pass в Python - это простая концепция, которую могут быстро освоить даже новички без опыта программирования.
Некоторые методы, о которых вы не знали, что они существуют в Python
Некоторые методы, о которых вы не знали, что они существуют в Python
Python - самый известный и самый простой в изучении язык в наши дни. Имея широкий спектр применения в области машинного обучения, Data Science,...
Основы Python Часть I
Основы Python Часть I
Вы когда-нибудь задумывались, почему в программах на Python вы видите приведенный ниже код?
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
Алиса и Боб имеют неориентированный граф из n узлов и трех типов ребер:
Оптимизация кода с помощью тернарного оператора Python
Оптимизация кода с помощью тернарного оператора Python
И последнее, что мы хотели бы показать вам, прежде чем двигаться дальше, это
Советы по эффективной веб-разработке с помощью Python
Советы по эффективной веб-разработке с помощью Python
Как веб-разработчик, Python может стать мощным инструментом для создания эффективных и масштабируемых веб-приложений.
2
0
2 119
2
Перейти к ответу Данный вопрос помечен как решенный

Ответы 2

Ответ принят как подходящий

Вы заметили, что ваши «компоненты» - это точно масштабированный исходный сигнал, перевернутый вверх ногами? Это потому, что вы не можете получить больше компонентов, чем сигналов.

Вам необходимо выполнить следующие действия:

  1. подать все каналы ЭЭГ в ICA
  2. вручную удалите компоненты, содержащие моргание глаз или другие артефакты
  3. реконструировать с помощью обратного преобразования

Давайте подробно рассмотрим шаг 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')

Почему мне нужно передавать каналы все в ICA? Что делать, если я использую только один канал? ICA работать не будет?

Pedro 06.11.2018 14:14

Даже вручную выбирая моргание глаз, я получаю только (151L,1L) float64 от ica_data = ica.fit_transform(blink). Полагаю, я не использую функцию ICA, как предполагалось ...

Pedro 06.11.2018 14:17

ICA преобразует каналы N в самое большее количество компонентов N. если * N * = 1 канал -> 1 компонент; это не очень полезно.

MB-F 06.11.2018 14:34

@Pedro Я обновил ответ конкретным примером, с которым вы можете поиграть. Это поможет продемонстрировать, как работает декомпозиция ICA.

MB-F 06.11.2018 15:10

Это фантастика @kazemakase! Большое спасибо за обновленный ответ и правило N-N ICA. Я смог сделать это с 4 каналами. Еще пара вопросов: 1) каждый раз, когда запускается код, компонент eyeblinks может быть другим. Как узнать, что это? Может свертка с блинк-шаблоном? Что ты предлагаешь? 2), что было бы лучшей альтернативой для жесткого кодирования компонента eyeblinks?

Pedro 06.11.2018 18:14
1) вы не знаете, поэтому я изначально предложил мануал «спросить шаг пользователя». Итак, вы запускаете алгоритм, затем появляется диалоговое окно с надписью пожалуйста, выберите мигающий компонент (вам нужно будет найти способ кодировать это), а затем алгоритм продолжается. Если у вас есть шаблон мигания, отлично - используйте свертку и порог; это может сработать. 2) Я бы записал моргания глаз. Поместите запасной электрод ЭЭГ возле глаза и запишите его вместе с ЭЭГ. Затем удалите компонент, который больше всего коррелирует с этим каналом.
MB-F 07.11.2018 10:31

Вижу, у меня еще много работы ... Вы ведь эксперт по ЭЭГ? Как вы смоделировали форму волны моргания (не использовали фильтр Баттуорта)? Эта форма волны может быть использована для соответствия сигналу.

Pedro 07.11.2018 20:59

Я больше не работаю в этой сфере. Что касается фильтра, я использовал Баттерворта, потому что это была первая функция, которая пришла в голову. Все, что мне было нужно, - это полоса пропускания, которую можно было применять на импульсных всплесках. Затем я поигрался с частотами фильтра и порядком, пока все не выглядело как-то правильно;)

MB-F 08.11.2018 09:19

В случае, если вы имеете дело с одноканальной ЭЭГ, простой в реализации метод может быть следующим:

1) Обнаруживайте мигания в сигнале Икс, используя простое определение пиков на основе пороговых значений (возможно, вам придется выяснить это, посмотрев на свои сигналы для нескольких случаев мигания). Устройства от Neurosky, Muse и т. д. Поставляются с API для обнаружения моргания. Вы можете использовать это при необходимости.

2) Возьмите кадр, соответствующий миганию (xb). Нарисуйте на нем плавную линию (xbs).

3) Вычтите плавную линию (xbs) из мигания (xb) и добавьте к ней среднее значение сигнала. Назовем это xbc.

4) Заменить xbc вместо xb в исходном сигнале Икс.

Другие вопросы по теме