Математические методы обработки многомерных временных рядов в применении к анализу электрофизиологических сигналов в реальном времени тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Сметанин Николай Михайлович

• Сметанин Николай Михайлович
• кандидат наук
• 2021

Оглавление диссертации кандидат наук Сметанин Николай Михайлович

1.1 Объект исследования

1.2 Цели и задачи исследования

1.3 Основные идеи, результаты и выводы диссертации

1.4 Теоретическая и практическая значимость приведенных результатов исследований

1.5 Вклад автора в проведенное исследование

1.6 Публикации и апробация работы

2 Содержание работы

2.1 NFBLab - платформа для проведения экспериментов с замкнутым контуром

2.2 Разработка низко-латентных методов оценки параметров ритмической активности мозга в режиме реального времени

2.3 Исследование влияния задержки на эффективность парадигмы нейрообрат-

ной связи

3 Заключение

3.1 Список выносимых на защиту результатов

3.2 Дальнейшие исследования

Список литературы

Приложения

А Статья "NFBLab—A Versatile Software for Neurofeedback and Brain-Computer

Interface Research"

B Статья "Digital filters for low-latency quantification of brain rhythms in real time" 58 C Статья "Short-delay neurofeedback facilitates training of the parietal alpha rhythm"

1 Введение

1.1 Объект исследования

Методы обработки многомерных временных рядов в применении к анализу электрофизиологических сигналов в реальном времени широко используются в исследованиях, относящихся к области нейронаук. Такие подходы необходимы в экспериментах, в которых возникает задача мгновенной оценки в реальном времени состояния активности центральной нервной системы (ЦНС). Регистрация данной активности как правило осуществляется при помощи набора сенсоров чувствительных к изменениям электромагнитного поля, создаваемого в результате работы ЦНС. Такими способами регистрации являются электроэнцефалография (ЭЭГ), магнитоэнцефалография (МЭГ), электрокортикография (ЭКоГ). В настоящей работе в качестве электрофизиологических сигналов будут подразумеваться сигналы ЭЭГ, однако все приведенные замечания справедливы также для МЭГ и ЭКоГ. В общем виде эксперименты, требующие мгновенной оценки состояния ЦНС, могут быть представлены в виде замкнутого контура, изображенного на рис. 1. А именно: (1) активность ЦНС участника эксперимента измеряется в режиме реального времени несколькими сенсорами, (2) многоканальные сигналы обрабатываются и выделяются некоторые целевые характеристики состояния ЦНС, (3) данные характеристики используются для формирования стимула для испытуемого, либо для управления программой или внешним устройством. Замыкание контура происходит в момент восприятия испытуемом стимула, либо в результате выполнения программой или устройством действия, каузально зависящего от результатов обработки измерений активности ЦНС.

В англоязычной литературе такого рода экспериментальные парадигмы обозначаются термином Closed-loop neuroscience. Примерами данной парадигмы с замкнутым контуром являются:

1. Интерфейс мозг-компьютер (ИМК), в котором активность мозга используется для прямого управления внешней программой или устройством [1]. Данная парадигма находит своё применение для обеспечения коммуникации с полностью парализованными пациентами, для замещения утраченной двигательной функции, а также в постинсультной нейро-реабилитации.

2. Нейро-обратная связь (НОС) - парадигма, в которой целевые характеристики активности мозга преобразуются в интерпретируемый человеком зрительный, слуховой или тактильный стимул [2]. Перед человеком ставится задача поддерживать стимул в определенном состоянии, например, стараться увеличивать высоту столбика отображенного на экране монитора, что соответствует поддержанию в нужном диапазоне

целевой характеристики активности мозга. В результате человек, проходящий сессию нейро-обратной связи, учится регулировать активность собственной ЦНС. Данная парадигма находит применение как для коррекции психоэмоционального состояния, так и для тренировки с целью повышения эффективности когнитивных функций, а также терапии широкого спектра нейродегенеративных заболеваний, в том числе эпилепсии.

3. Стимуляция (транскраниальная магнитная, постоянным, переменным током) активности мозга, зависящая от текущего состояния ЦНС [3, 4]. В такой парадигме отслеживается активность мозга в режиме реального времени и на основании оценки ее параметров принимается решение о моменте включения стимуляции, которая может быть направлена на подавление патологической активности или вызов определенной поведенческой реакции, что используется, например, для подавления тремора у больных Паркинсонизмом или сокращения вероятности эпилептического припадка.

4. Поведенческие эксперименты с онлайн мониторингом активности мозга, в которой стимулы или задания предъявляются в моменты, когда ЦНС находится в определенном состоянии [5]. Такого рода парадигмы являются альтернативой экстенсивному подходу, при котором стимулы предъявляются в случайные моменты времени, и затем при пост-обработке отбираются предъявления попадающие во временные интервалы, на которых присутствует нужный тип активности. Использование онлайн-мониторинга позволяет существенно сократить эксперимент и упростить интерпре-

Рис. 1: Схематическое изображение парадигм с замкнутым контуром, применяемых в исследованиях области нейронаук (Closed-loop neuroscience)

тацию результатов.

Следует отметить, что способы предъявления обратной связи в перечисленных парадигмах обладают значительным разнообразием - от стимула на экране до электромагнитного импульса или фазовой синхронизации прямой электрической стимуляции. При этом этап обработки сигналов часто является схожим в данных парадигмах и представляет собой основной объект исследования данной работы.

Пространственная и временная фильтрация представляют собой две базовые операции, используемые в обработке многоканального ЭЭГ сигнала.

Пространственная фильтрация заключается в приведении многомерного многоканального сигнала к одномерному временному ряду. Такое преобразование может быть осуществлено тривиально, при помощи выбора одного "реального" канала (или отведения, как его называют в ЭЭГ литературе) для последующего анализа. Более общий способ подразумевает построение "виртуального" канала, которое как правило выполняется линейным объединением нескольких каналов в один и называется линейной пространственной фильтрацией. Далее одномерный временной ряд после этапа пространственной фильтрации будет называться виртуальным отведением. Перед исследователем ставится задача определения коэффициентов пространственного фильтра исходя из требований экспериментальной парадигмы. Известно, что мозг человека можно разделить на области по функциональному признаку, например, зрительная система расположена в окци-питальной (затылочной) зоне, моторный контроль относится к областям у центральной извилины, фронтальная кора отвечает за различные когнитивные функции. Целью пространственной фильтрации как правило является выделение сигнала активности требуемой области мозга. Также следует заметить что пространственная фильтрация позволяет перейти от сигналов ЭЭГ-сенсоров к сигналам нейрональных источников, тем самым снизить эффект объемного проводника (отображение одного источника на нескольких сенсорах), который затрудняет анализ данных с сенсоров напрямую.

Пространственный фильтр может быть получен в результате решения обратной задачи [6] или использования пространственного разложения многоканальных электрофизиологических данных [7]. В настоящий момент разработан широкий спектр методов построения пространственных фильтров. В основном данные методы применяются для постобработки данных. Для применения данных методов онлайн, т.е. непосредственно в экспериментах, необходимо создание платформы, позволяющей производить гибкую настройку параметров методов и получаемых пространственных фильтров. Пространственная фильтрация в экспериментах с замкнутым контуром в сравнении с прямым использованием данных с физических ЭЭГ отведений позволяет сфокусироваться на функционально специфических источниках нейрональной активности, а также отстроится от артефактов и

фоновой активности мозга. В случае парадигм с замкнутым контуром пространственная фильтрация позволяет повысить пространственную специфичность сигнала обратной связи.

Временная фильтрация позволяет при помощи вычисления скользящей линейной комбинации временных отсчетов сигнала одного отведения выделить определенный временной паттерн активности ЦНС. Одним из наиболее исследованных компонентов нейро-нальной активности являются ритмы головного мозга являются [8]. Первым из открытых ритмов мозга является альфа-ритм - колебания с центральной частотой 8-12 Гц, обнаруженные Гансом Бергером в 1924 году. Например, локализованный в окципитальной коре альфа-ритм отражает состояние зрительных отделов мозга: при закрытых глазах, когда на зоны зрительной коры фактически не поступает визуальная информация, мощность данного ритма возрастает, при открытых уменьшается. Другим важным ритмом, широко используемым, например, в парадигмах мозг-компьютерного интерфейса является сен-сомоторный ритм (или мю ритм), локализованный в районе центральной борозды. Центральная частота данного ритма лежит в диапазоне 10-14 Гц и по аналогии со зрительной системой данный ритм отражает состояние сенсомоторной системы: в состоянии покоя и моторного бездействия данный ритм увеличивается, во время выполнения, наблюдения и воображения движений определенной части тела амплитуда ритма уменьшается в зонах представительства данного органа в сенсомоторном отделе коры головного мозга.

Ритмическая активность может быть представлена как узкополосный, частотно - модулированный сигнал, основными характеристиками которого являются мгновенная фаза и амплитуда осцилляций [9]. В соответствии с современными представлениями перепад амплитуд и фаза ритмической активности головного мозга являются фундаментальными параметрами, связанными с состоянием ЦНС [8, 10]. Для парадигм с замкнутым контуром ставится задача оценки мгновенной фазы и/или амплитуды (огибающей) узкополосного сигнала по сырому, в общем случае широкополосному сигналу в режиме реального времени. Следует отметить тот факт, что любая фильтрация в режиме реального времени вносит задержку, связанную с применением каузальных фильтров. Такая задержка имеет сугубо фундаментальную природу и является проявлением принципа неопределенности Габора [11] в задачах обработки сигналов, в соответствии с которым невозможно с одинаковой точностью локализовать сигнал на частотной и временной осях. Эту задержку не стоит путать с техническими задержкой, включающей время трансфера данных между записывающим устройством и компьютером, время выполнения вычислений и время генерации стимула исполнительным устройством, например монитором. Как правило, с использованием современных программно-аппаратных средств, технические задержки могут сведены к минимуму и принимать значения значительно меньшие в сравнении с фундаментальной

задержкой, возникающей вследствие каузальной временной фильтрации сигнала.

Таким образом в системах реального времени существует задержка между возникновением целевой активности мозга и моментом, когда это событие отражается в сигнале обратной связи. Ряд приложений в рамках парадигмы с замкнутом контуром требуют достижения минимально возможной задержки, что зачастую приводит к ухудшению качества оценки параметров ритмической активности. В зависимости от типа экспериментальной парадигмы данная задержка должна быть либо минимизирована, либо учтена при генерации стимула для обеспечения требуемой временной специфичности.

Снижение временной специфичности ведет к низкой эффективности всей парадигмы. Например, увеличение задержки между ментальной инициацией (началом воображения) движения и началом движения курсора программы в реализации ИМК ведет к снижению ощущения управляемости [12]. В случае парадигмы нейрообратной связи анализ симу-ляционных данных [13] показывает, что задержка или случайное смещение во времени обратной связи отрицательно влияют на скорость обучения. Кроме того, как будет показано нами, сокращение временной задержки предъявления сигнала НОС ведет к повышению эффективности тренировки в парадигме нейрообратной связи [14]. Таким образом перед разработчиком систем с замкнутым контуром ставится задача сокращения математической части задержки к фундаментально-возможному минимуму, при этом сохраняя качество оценки, огибающей или фазы. Помимо классических подходов к этой задаче [9] к настоящему моменту развиты специализированные методы, например [15, 16]. Однако данные подходы представляют собой сложные в реализации эвристики, зависящие от большого числа параметров и не предусматривающие возможности эксплицитного регулирования параметра задержки системы в сочетании с обеспечением лучшего возможного качества оценки параметров ритмической активности головного мозга.

1.2 Цели и задачи исследования

Неоптимальность тракта обработки данных в парадигме с замкнутым контуром ведет к снижению эффективности различных реализаций этой парадигмы, затрудняя ее применение и ограничивая потенциальные возможности парадигмы. Основная цель настоящей работы заключается в разработке методов и программных средств обработки электроэнцефалографических данных для применения в режиме реального времени и нацеленных на повышение эффективности реализации парадигм с замкнутым контуром.

Следует отставить несколько замечаний про способы оценки эффективности реализации парадигм с замкнутым контуром. Эффективность парадигмы НОС как правило оценивается приростом целевого сигнала во время состояния покоя сразу после тренировки относительно состояния покоя до тренировки. Также показателем эффективности

является скорость роста целевого сигнала во время тренировки. Кроме того, может быть произведено измерение прироста различных поведенческих характеристик после эксперимента. Эффективность ИМК может быть измерена точностью распознавания состояния испытуемого, точностью выполнения испытуемым заданий, требующих управление через ИМК, а так же, например, при помощи опроса испытуемого для оценки чувства контроля ИМК. Если экспериментальная парадигма заключается в мониторинге и детекции моментов для стимуляции то эффективность реализации тракта обработки данных может быть оценена точностью определения целевых состояний. Кроме того, во многих экспериментах стимуляция происходит в случайные моменты времени, а затем на этапе пост-обработки выделяются стимулы предъявленные в целевые моменты времени. Правильная реализация алгоритмов онлайн детекции целевых состояний ведет к сокращению длительности эксперимента, поскольку необходимое для анализа число предъявлений задачи или стимула набирается быстрее. Сокращение времени эксперимента уменьшает утомляемость испытуемого и как следствие повышает качество полученных данных. Таким образом в качестве метрики эффективности в такого рода экспериментах можно оценивать сокращение времени эксперимента при использовании онлайн детекции в сравнении с классическим экспериментом в которых стимул не привязан к активности мозга.

Как было отмечено ранее эффективность реализации перечисленных парадигм во многом зависит от точности методов оценки активности целевых нейрональных источников. А именно, методы обработки данных должны обладать двумя ключевыми свойствами: (1) высокая пространственная специфичность - высокая точность пространственной локализации целевого источника и отстройка от фоновой активности мозга и внешних артефактов, (2) высокая временная специфичность - низкая задержка между изменением в активности целевого источника и моментом, когда это изменение появилось в сигнале обратной связи.

Диссертационное исследование разделено на 3 связанных проекта - программный, методологический и экспериментальный. Программный проект посвящен разработке платформы для реализации парадигм с замкнутым контуром с возможностью имплементации методов удовлетворяющим вышеперечисленным свойствам (1) и (2). Методологический проект направлен на непосредственное развитие подходов к оценке параметров ритмической активности мозга в режиме реального времени, отвечающих свойству (2). Задачей экспериментального проекта является исследование влияния временной специфичности -свойство (2) - на эффективность реализации парадигм с замкнутым контуром.

Перечисленные проекты объединены общей тематикой «Математические методы обработки многомерных временных рядов в применении к анализу электрофизиологических сигналов в реальном времени» и представляют собой законченное исследование, в

результате которого были разработаны программно-алгоритмические средства обработки многоканальных ЭЭГ сигналов и продемонстрирована их эффективность при построении идеомоторных ИМК и в парадигме НОС. В настоящее время все разработанные средства и алгоритмы используются в научно-исследовательской деятельности Центра Биоэлектрических Интерфейсов НИУ ВШЭ.

1.3 Основные идеи, результаты и выводы диссертации

В рамках программного проекта была разработана программная платформа КЕБЬаЬ для реализации парадигм с замкнутым контуром. Это программное обеспечение позволяет: (1) настраивать тракт обработки данных, включая пространственную и временную фильтрацию, в том числе с возможностью индивидуализированной настройки параметров фильтров по записанным функциональным пробам, (2) гибко формировать дизайн эксперимента, а именно последовательность блоков с указанием их типов и параметров обработки сигнала в каждом их блоков, (3) проводить эксперимент в парадигме с замкнутым контуром, обеспечивая подключение, прием, запись и обработку многоканальных электрофизиологических данных, а также генерацию стимулов в том числе и с минимально-возможной задержкой генерации сигнала ОС. Платформа включает в себя специально разработанный скриптовый язык, позволяющий полностью описывать параметры тракта обработки сигналов и дизайн эксперимента, состоящий из последовательности блоков, включая возможность рандомизации их последовательности. Разработанная платформа включает в себя как традиционные, так и вновь разработанные методы обработки данных и использовался в последующих проектах в качестве основной средства для тестирования разрабатываемых методов и проведения экспериментов в парадигмах с замкнутым контуром в Центре Биоэлектрических Интерфейсов НИУ ВШЭ.

В рамках методологического проекта было разработано семейство методов для низколатентной оценки параметров ритмической активности мозга в режиме реального времени. Разработанные подходы обеспечивают возможность эффективной оценки огибающей и фазы узкополосной компоненты широкополосного электрофизиологического сигнала. При этом полная задержка тракта обработки является независимым параметром, для каждого из значений которого разработанные подходы обеспечивают оптимальную точность оценки огибающей и мгновенной фазы узкополосного процесса. Проектирование разработанных фильтров основано на решении оптимизационной задачи поиска каузального комплекснозначного фильтра с конечной импульсной характеристикой (КИХ), аппроксимирующего некаузальный идеальный фильтр оценки аналитического узкополосного сигнала. В сравнении со стандартными методами, применяемыми в парадигмах с замкнутым контуром [15], разработанное семейство методов обеспечивает меньшее значение задержки

выделения огибающей и фазы при сохранении точности оценок. При этом у пользователя появляется возможность эксплицитно настраивать соотношение задержка-точность в контексте конкретного приложения.

В рамках экспериментального проекта разработанная методология и программное обеспечение было применено к исследованию влияния временной специфичности на эффективность реализации парадигм с замкнутым контуром. В частности, проверялась гипотеза о негативном влиянии задержки обратной связи на эффективности обучения в парадигме нейрообратной связи (НОС). Проведенный эксперимент включал 4 группы испытуемых. Испытуемым из первых трех групп предъявлялся стимул обратной связи с полной задержкой в 250, 500 и 750 мс. Четвертая группа испытуемых получала ложную обратную связь и использовалась в качестве контрольной группы. По результатам эксперимента была выявлена статистически достоверная взаимосвязь эффективности и задержки НОС системы, при этом чем меньше задержка НОС, тем быстрее происходило обучение, а также было более выражено сохранение результатов тренировки по ее окончанию.

1.4 Теоретическая и практическая значимость приведенных результатов исследований

Платформа КЕБЬаЬ, разработанная в рамках программного проекта, позволяет проектировать и проводить эксперименты с замкнутым контуром и основана на разработанном формате файла-описания тракта обработки сигналов и дизайна эксперимента. Основным методом пространственной фильтрации в КЕБЬаЬ является создание пространственных фильтров на основе декомпозиции записанных во время эксперимента функциональных проб. В сравнении со стандартными методами основанными на решении обратной задачи [6], такой подход является более индивидуальными и позволяет выделять зоны мозговой активности согласно их функциональному поведению и при этом не требует выполнения достаточно трудоемких вычислений, связанных с построением индивидуализированной прямой электромагнитной модели и требующих сегментации МРТ (магнитно-резонансная томография) испытуемого. Тем не менее, КЕБЬаЬ предоставляет пользователю воспользоваться и такой общепринятой методологией посредством эффективного интерфейса с пакетом МКЕ-Ру^оп [17]. Тракт обработки сигналов использует механизм композитных сигналов, при помощи которого оказывается возможным рассчитывать в реальном времени меры функционального взаимодействия участков коры головного мозга. Кроме того, большое внимание на стадии разработки КЕБЬаЬ уделялась проблеме сокращения латентности в контуре обработки сигнала. Математическая задержка применяемых методов выступает

в роли независимого параметра, значение которого задается пользователем при проектировании тракта обработки сигнала. Присутствует возможность вносить дополнительную искусственную задержку для экспериментов по изучению влияния латентности системы на эффективность парадигмы замкнутого контура, а также реализовывать парадигмы идеомоторного ИМК. Программная платформа также предоставляет возможность проектировать эксперименты с практически произвольным дизайном, последовательностью, длительностью блоков, включая рандомизацию их порядка и использование методов статистической нормировки сигналов специфичных каждому из таких блоков. Разработанная ПО является платформой тестирования новых разработанных в рамках диссертации методов низколатентной оценки параметров ритма. Данная программная платформа является проектом с открытым исходным кодом, написанный на языке python [18] и продолжает развиваться в настоящее время командой разработчиков и пользователей по всему миру.

Разработанная платформа NFBLab является уникальным программным обеспечением, которое открывает международному сообществу возможность проводить воспроизводимые эксперименты в рамках парадигмы замкнутого контура, при этом предоставляя возможность гибкой настройки параметров тракта обработки сигналов и экспериментального дизайна. Особое внимание в данном ПО уделено контролю за задержкой в контуре обработке сигнала. Как было показано авторами данного исследования [14] такая задержка оказывает существенное влияние на эффективность тренировки в парадигме нейрообрат-ной связи и минимизация данной задержки открывает ранее недоступные возможности по формированию петли взаимодействия с человеческим мозгом, оперирующей на скорости работы самого мозга.

Это стало возможным при помощи разработанного в данном исследовании семейства низко-латентных методов оценки параметров ритмической активности головного мозга, которые позволяют сокращать задержку оценки сигнала обратной связи с сохранением качества оценки в сравнении с классическими методами. Подходы из предлагаемого семейства позволяют учесть нестационарность сигналов ритмической активности головного мозга, а также по построению обладают способностью адаптироваться к произвольной форме спектра таких сигналов. Кроме того, новые методы позволяют в явном виде управлять задержкой тракта обработки сигналов, которая является независимым параметром методов из разработанного семейства. Таким образом предлагаемые методы могут быть использованы в качестве основных подходов для мгновенной оценки фазы или огибающей целевой ритмической активности в широком спектре парадигм с замкнутым контуром, в которых минимизация задержки является критическим требованием. Автором данного диссертационного исследования было продемонстрировано, что предлагаемое семейство методов обеспечивает меньшую задержку и более высокую точность оценки параметров

ритмической активности мозга [19], чем существующие подходы [15]. Повышение временной специфичности позволяет существенно повысить эффективность реализации данных экспериментальных парадигм.

Проект по исследованию влияния временной специфичности обратной связи на эффективность реализации парадигм с замкнутым контуром является первой в мире попыткой систематического изучения этого явления в парадигме нейрообратной связи (НОС). Ранее, факт влияния задержки сенсорной обратной связи на эффективность обучения, был подтвержден в широком спектре поведенческих исследований. Например, еще в 1948 году Опсе показал, что эффективность обучения в задаче дискриминации сложных визуальных паттернов существенно зависит от задержки сигнала обратной связи [20]. В работе [21] показано, что поведенческие корреляты обучения ухудшаются, при неизвестной задержке предъявления обратной. Увеличение задержки обратной связи значительно ухудшает ощущение вовлеченности и уменьшает ощущение авторства во время управления внешними устройствами при помощи интерфейса мозг-компьютер [12], который по сути своей является одной из реализаций парадигмы с замкнутым контуром. Однако практически во всех реализациях парадигм нейрообратной связи до настоящего времени уделялось недостаточно внимания контролю за значением параметра задержки при предъявлении сигнала обратной связи. Данное исследование по сути первая в мире попытка разработки методологии обратной связи с эксплицитно управляемой и контролируемой задержкой в предъявлении сигнала обратной связи. Предпосылки необходимости такой работы были получены ранее.

В недавнем исследовании [22] были изучены изменения во временной структуре ЭЭГ, вызванные НОС. Альфа-ритм ЭЭГ регистрировался на теиенном отведении Р4, и средняя мощность этого сигнала предъявлялась испытуемым посредством визуальной обратной связи. Задачей испытуемого было повысить уровень средней мощности альфа-ритма. Анализ эпизодов высокой мощности альфа-ритма показал, что испытуемые не могли модулировать амплитуду ритма или продолжительность поддержания состояния высоко амплитудного альфа-ритма. Вместо этого повышение средней мощности альфа-ритма в сессии достигалось исключительно за счет увеличения удельного числа эпизодов вхождения в целевое состояние, характеризующееся высоким значением мгновенной амплитуды альфа-колебания. В результате была выдвинута гипотеза, что для повышения эффективности НОС, вместо того чтобы использовать среднюю мощность ритма на всем временном промежутке, эпизоды вхождения в состояние высокой синхронизации альфа-ритма должны рассматриваться как дискретные события. Важность дискретной компоненты в интерпретации нейрональной активности была также продемонстрирована в [23], где авторы показали, что количество веретен бета-ритма в единицу времени (а не амплитуда и не

длительность) определяет эффективность выполнения моторной задачи. В соответствии с полученными результатами, можно сделать вывод о том, что начало и конец вхождения в целевое состояние являются специфически значимыми событиями, подкрепление воспроизведения которых может привести к более специфичной тренировке в парадигме НОС. Учитывая, что характерная длина всплесков альфа-активности лежит в диапазоне 200-300 мс можно высказать предположение о важности временной специфичности как одного из факторов, влияющих на эффективность реализации парадигмы НОС. Программное обеспечение и методология обработки сигналов, разработанные в данной работе позволили провести первое в мире систематическое исследование [14], которое подтвердило исключительную важность параметра задержки предъявления сигнала обратной связи на эффективность тренировки в парадигме НОС.

1.5 Вклад автора в проведенное исследование

Автор настоящего исследования является главным разработчиком платформы NFBLab [24]. Семейство низколатентных алгоритмов было сформулировано и исследовано автором работы [25]. В экспериментальном исследовании [14] автор обеспечивал методологическую и техническую поддержку системы НОС, а также принимал активное участие в обработке результатов эксперимента. По результатам работы опубликованы три статьи в Q1 журналах, в двух из которых диссертант является первым автором, в третьей статье вклад диссертанта равен вкладу первого автора.

Введение диссертации (часть автореферата) на тему «Математические методы обработки многомерных временных рядов в применении к анализу электрофизиологических сигналов в реальном времени»

1.6 Публикации и апробация работы

Публикации повышенного уровня

1. Smetanin N, Volkova K, Zabodaev S, Lebedev MA and Ossadtchi A (2018) NFBLab—A Versatile Software for Neurofeedback and Brain-Computer Interface Research. Front. Neuroinform. 12:100. doi: 10.3389/fninf.2018.00100

2. Smetanin N, Belinskaya A, Lebedev M, Ossadtchi A. Digital filters for low-latency quantification of brain rhythms in real time. J Neural Eng. 2020;17(4):046022. Published 2020 Aug 4. doi:10.1088/1741-2552/ab890f

3. Belinskaia, A., Smetanin, N., Lebedev, M., Ossadtchi, A. (2020). Short-delay neurofeedback facilitates training of the parietal alpha rhythm. Journal of Neural Engineering, 17(6), 066012.

Конференции

1. Limitless: Augmenting brain function, г. Лозанна, Швейцария, 19 - 21 октября 2018 г. Доклад: "Augmenting the brain with temporally-structured neurofeedback". Осадчий А.Е., Сметанин Н.М., Волкова К.В. Лебедев М.А.

2. Cell-NERF Symposium: Neurotechnologies, г. Лёвен, Бельгия, 29 сентября - 02 октября

2018 г. Доклад: "Online and offline modulation of sensorimotor components following focal vibration". Булгакова В.О., Сметанин Н.М., Волкова К.В., Осадчий А. Е

3. Международная конференция Нейрокомпьютерный Интерфейс: Наука и практика, г. Самара, Россия, 11-12 ноября 2018 г. Доклад: "Towards zero-latency neurofeedback". Осадчий А.Е., Сметанин Н.М., Белинская А.А.

4. 5-я ежегодная конференция СССР Коды мозга: управление и восприятие, Москва, Россия, 29 - 30 ноября 2018 г. Доклад: "К нейрообратной связи с нулевой задержкой". Осадчий А.Е., Сметанин Н.М, Белинская А.А.

5. Форум Skolkovo robotics, Россия, Москва, 16 апреля 2019 г. Доклад: "Интерфейсы мозг - компьютер: декодирование ходьбы". Лебедев М.А., Сметанин Н.М

6. OHBM Annual Meeting, г. Рим, Италия, 09 - 13 июня 2019 г. Доклад: "From low-latency to predictive neurofeedback: methods and feasibility check". Сметанин Н.М, Осадчий А.Е.

7. OHBM Annual Meeting, г. Рим, Италия, 09 - 13 июня 2019 г. Доклад: "The effect of feedback latency on the effectiveness of training in neurofeedback paradigm". Белинская А.А., Сметанин Н.М., Осадчий А.Е.

8. The Second Neuroadaptive Technology Conference NAT 2019, г. Ливерпуль, Великобритания The Second Neuroadaptive Technology Conference NAT 2019, 16 - 18 июля

2019 "Differential effects of neurofeedback latency on the incidence rate, amplitude and duration of alpha-bursts". Сметанин Н.М., Белинская А. А , Осадчий А. Е. - третий приз за лучший доклад.

9. Международная конференция Нейрокомпьютерный Интерфейс: Наука и практика, г. Самара, Россия, 02 - 05 ноября 2019 г. Доклад: "Software Platform for MEG-Based Neurofeedback Training". Сметанин Н.М.,Осадчий А.Е.

10. Международная конференция Нейрокомпьютерный Интерфейс: Наука и практика, г. Самара, Россия, 02 - 05 ноября 2019 г. Доклад: "Foot motor imagery triggered

locomotion in exoskeleton: first results with paraplegic patients". Сметанин Н.М, Кузнецова А.А., Осадчий А.Е.

11. Открытый семинар проекта #CNBR_Open, Инновационный центр "Сколково", г. Москва, Россия, 12 февраля 2020 г. Открытая лекция и демонстрация: "Demystifying brain-computer interfaces". Алексей Осадчий, Николай Сметанин.

Патенты и авторские права

1. Изобретение. RU 2713110 C1. Способ оценки различий мощности осцилляторных компонент сигналов электроэнцефалограммы в психофизиологических состояниях на основе квантильного анализа. 03.02.2020

2. Программа для ЭВМ. RU 2020613238. Когниграф. 12.03.2020

3. Программа для ЭВМ. RU 2020618826. Программная реализация парадигмы низколатентной нейрообратной связи с визуальной стимуляцией на основе регистрации электроэнцефалограммы человека "Низколатентная НОС". 05.08.2020

4. Программа для ЭВМ. RU2018611760. Среда для проведения тренинга в парадигме нейрообратной связи. 06.02.2018

2 Содержание работы

2.1 NFBLab - платформа для проведения экспериментов с замкнутым контуром

В данном разделе приведено краткое содержание опубликованной работы [25], в которой дано подробное описание разработанного программного обеспечения NFBLab. На ранних этапах построения данного ПО прежде всего ставилась задача поддержки экспериментов в парадигме нейрообратной связи с учетом обеспечения воспроизводимости и гибкости формирования тракта обработки сигнала и экспериментального дизайна. Однако, по мере развития эта разработка превратилась в полноценную программную платформу, позволяющую реализовывать широкий спектр парадигм с замкнутым контуром.

Назначение и свойства платформы Для реализации парадигм с замкнутым контуром на практике используются различные программные решения, реализующие подключение к ЭЭГ/МЭГ устройству, обработку получаемых многоканальных данных с целью отстройки от артефактов и выделения целевого сигнала с последующим предъявлением сигнала ОС посредством одной из сенсорных модальностей человека или использования такого сигнала для управления внешними устройствами. Одним из таких решений является разработанная авторами программная платформа NFBLab, краткому описанию которой посвящена данная часть реферата диссертационного исследования. NFBLab можно охарактеризовать как программное обеспечение (ПО) для проведения экспериментов в парадигмах с замкнутым контуром на основе стандартных и оригинальных низколатентных алгоритмов обработки многоканальных биоэлектрических сигналов. Ключевыми особенностями данного ПО являются:

• поддержка низколатентного подключения к большинству распространенных в настоящий момент ЭЭГ/МЭГ устройств посредством протокола на основе сокетной технологии Lab Streaming Layer [26];

• наличие внутреннего псевдоязыка позволяющего реализовать гибкую настройку сценария эксперимента при помощи xml-псевдокода, обеспечивающего воспроиводимость и автоматическую документацию проведенных экспериментов, а также позволяющего реализовать автоматический запуск программы всего эксперимента

• возможность изменения сценария эксперимента посредством редактирования xml-псевдокода или при помощи графического интерфейса, а также методов графического программирования;

• наличие интерфейса для проведения эксперимента, включая интерактивный модуль для формирования тракта обработки сигналов на основе функциональных проб и пространственно-временных разложений;

• возможность воспроизведения выполненного эксперимента;

• возможность визуализации извлеченных целевых сигналов, а также наличие программного интерфейса для подключения к извлеченным сигналам сторонних программ предъявления ОС, управления исполняющими устройствами и игр;

• наличие возможности гибкой настройки экспериментального дизайна, как последовательности блоков, механизма рандомизации последовательности экспериментальных блоков и формирования поддельного сигнала обратной связи (mock feedback) для проведения научных экспериментов с использованием контрольных групп

• наличие реализованных новых авторских алгоритмов для низколатентной оценки мгновенной мощности узкополосных сигналов [19]

• открытый код на языке Python, кроссплатформенный интерфейс

Аналоги NFBLab Наиболее популярными проектами для проведения экспериментов с ЭЭГ в реальном времени, развивающимися и поддерживаемыми в настоящее время, являются OpenVIBE [27] и BCI2000 [28]. Область применения данных платформ - эксперименты с применением различного рода парадигм обработки и визуализации био-сигналов в режиме реального времени. В отличие от OpenVIBE и BCI2000, NFBLab позволяет не только конфигурировать тракты выделения целевых сигналов ограниченной сложности, но при этом содержит в себе модуль управления всем экспериментом и переключает и при необходимости повторяет, а также рандомизирует экспериментальные блоки автоматически. Так же в отличие от перечисленных платформ проект NFBLab в качестве основных методов пространственной фильтрации использует индивидуализированный и наиболее пространственно специфичный подход - метод декомпозиции функциональных проб. Так же в NFBLab внедрены методы, повышающие временную специфичность тракта обработки данных. Кроме того, NFBLab распространяется с открытым кодом на языке Python [18], что позволяет продвинутым пользователям реализовать новые протоколы и модули обработки сигнала и предполагает дальнейшее развитие проекта усилиями мирового сообщества.

Архитектура NFBLab состоит из трёх основных модулей. Первый модуль "Редактор протокола эксперимента" (Experiment protocol editor) позволяет формировать сценарий

эксперимента. Полученный дизайн, включающий в себя описание трактов обработки сигналов, виртуальных отведений, целевых сигналов, а также параметров вычисления сигнала обратной связи для каждого из экспериментальных блоков и последовательность этих блоков, включая схему рандомизации и параметры нормировочной статистики, представляет собой программу эксперимента в форме псевдокода и сохраняется в .xml файле, который может быть загружен для дальнейшего повторного использования и проведения стереотипного эксперимента. Второй модуль "Модуль проведения эксперимента" (Experiment module) запускается при старте эксперимента, обрабатывает и отображает сырые и целевые сигналы в режиме реального времени, вычисленные в соответствии с построенным выше псевдокодом, управляет последовательностью блоков, а также предъявляет различного роды стимулы. Третий модуль "Модуль настройки фильтров по собранным данным" (Data-driven filter designer) является интерактивным модулем для редактирования свойств тракта обработки сигналов и построения пространственно-временных фильтров на основе анализа записанных данных, в том числе и функциональных проб, посредством их частотного анализа и пространственного разложения различными методами (см. разделы 4.1 и 6.3 в приложении A). Как правило, этот модуль выполняется непосредственно во время эксперимента, приостанавливает работу предыдущего модуля и служит для создания индивидуализированных трактов обработки сигнала и обновления xml-псевдокода эксперимента. Получение данных с ЭЭГ/МЭГ устройства осуществляется с помощью технологий Lab streaming layer (LSL) или FieldTripBuffer (FTB) (см. раздел 3.1 в приложении A). Записанные данные эксперимента, включая все целевые сигналы, вычисленные на основании описанных в xml-пседокоде трактов обработки сигнала, сохраняются в файл формата hdf5 (см. раздел 6.2 в приложении A). Также, вычисленные в режиме реального времени целевые сигналы отправляются в LSL Outlet для связи с внешними программами/устройствами (см. раздел 3.1 приложения A). Схема связи модулей представлена на рисунке 2A. Более подробное описание модулей рассмотрено в приложении A и статье [25].

Обработка сигналов в NFBLab В NFBLab реализована возможность обработки в режиме реального времени сырых энцефалографических сигналов, вычисления виртуальных отведений, фильтрации последних в заданном частотном диапазоне с последующей оценкой мгновенной мощности ритмов мозга (целевые сигналы) и произвольных функций целевых сигналов, задаваемых математическим выражением (композитные сигналы). На рисунке 2B схематически отображены основные этапы выделения целевых/композитных сигналов.

Под виртуальным отведением здесь и далее будем понимать линейную комбинацию сигналов сенсоров. Вектор коэффициентов линейной комбинации может носить триви-

Рис. 2: Схема связи основных модулей в NFBLab (А) и этапы обработки сигналов (В

альный характер и состоять из всех нулей за исключением одной позиции, в которой находится единица. В этом случае виртуальное отведение полностью совпадает с реальным отведением от электрода с номером, соответствующем позиции ненулевого элемента вектора весов. Вычисление линейной комбинации с весами, полученными на основании, например, решения обратной задачи и соответствующими определенной строке обратного оператора, даст виртуальное отведение, отражающее активность соответствующей области коры. Альтернативным методом нахождения весов для формирования виртуального отведения может служить использование функциональных проб в сочетании с математическими методами обработки сигналов для выделения проекций максимального контраста или проекций максимальной мощности. Так или иначе, при вычислении виртуального отведения используются весовые коэффициенты для пространственного выделения компоненты по анатомическому или функциональному признаку. Часто, компоненты, выделяемые по функциональному признаку, обладают, также, и пространственной специфичностью. Поэтому вычисление виртуального отведения в терминах фильтрации, может быть названо пространственной фильтрацией.

Более формально, пространственная фильтрация выражается следующим соотношением: уЩ = х[£], где х[£] - вектор-столбец многоканальных измерений в момент времени £, уЩ - значение сигнала на виртуальном отведении в момент времени £, ш - вектор-столбец

коэффициентов пространственного фильтра, количество элементов которого равняется количеству каналов записи. Пространственный фильтр может быть представлен в виде произведения двух компонент w = Ru. Режекторная матрица R представляет собой как правило идемпотентную матрицу проекции и используется для отстройки от некоторого паттерна физиологической активности (например, для исключения глазных артефактов). Вектор-столбец u действует противоположно R и служит для выделения требуемой активности, формируя сигнал виртуального отведения. Для нахождения режекторных матриц и пространственных фильтров в NFBLab применяются различные методы пространственного разложения:

1. Метод ICA (Independent Component Analysis) применяется для разложения сигнала на независимые компоненты и используется для выделения и удаления различного рода артефактов [29].

2. Метод CSP (Common Spatial Pattern) позволяет выделять компоненты с максимальным отношением мощности сигнала для двух окон (например, первое окно может соответствовать первой половине записи, в котором испытуемый находится с закрытыми глазами, второе - с открытыми). Ключевой частью алгоритма является решение обобщенной задачи на собственные значения [30].

3. В методе SSD (Spatio-Spectral Decomposition) разложение происходит также с помощью решения обобщенной задачи на собственные значения. Данный метод позволяет выделять компоненты с максимальным отношением мощности сигнала для двух различных частотных полос (центральная полоса и две прилежащие частоты), что позволяет выделять узкополосные осцилляторные компоненты [31]

Как правило, исследователь заинтересован в отслеживании динамики ритмической активности мозга, которая находит отражение в ЭЭГ в виде узкополосного случайного процесса. Ритмы мозга имеют характерные частоты и следующим элементом тракта обработки сигналов является цифровой частотный фильтр, коэффициенты которого вычисляются на основании указанной пользователем частотной полосы и порядка фильтра. Ритмическая активность мозга носит нестационарный характер, и может быть охарактеризована как последовательность всплесков [32]. Мгновенная интенсивность такой активности описывается огибающей узкополосного процесса, вычисляемой, например, при помощи преобразования Гильберта. Таким образом, следующим элементом тракта обработки сигнала является вычисление огибающей сигнала. В NFBLab помимо классических подходов реализован новый метод оценки огибающей узкополосного процесса опубликованный в работе [19] и описанный в следующем разделе диссертационного исследования.

Новый метод оценки огибающей позволяет сократить латентность предъявления сигнала обратной связи при сохранении качества. В результате описанных действий получается целевой сигнал, который представляет собой последовательность значений, отражающую мгновенную мощность сигнала в виртуальном отведении в заданной полосе частот.

Такие сигналы будем называть целевыми (Derived signal) и для их описания в XML файле предусмотрена определенная структура. Свойства используемых фильтров сигналов можно определить в модуле "Дизайн эксперимента" до начала эксперимента или изменить их на основе записанных функциональных проб в интерактивном модуле "Настройка фильтров по собранным данным и функциональным пробам" во время эксперимента. Кроме того, для продвинутых пользователей возможно простое редактирование текста XML файла.

Часто, в качестве сигнала обратной связи выступает соотношение мощности в различных реальных или виртуальных отведениях. В NFBLab реализован класс композитного сигнала (Composite signal), который определен как произвольная математическая функция двух целевых сигналов. Вид функции задается экспериментатором. Например, для построения протокола тренировки, в котором обуславливается отношение фронтального бета-ритма к затылочному альфа-ритму, необходимо создать два целевых сигнала, соответствующих фронтальному бета-ритму и затылочному альфа-ритму, а затем создать композитный сигнал, объединяющий два целевых сигнала при помощи функции деления. Помимо вычисления произвольных математических функций механизм композитных сигналов позволяет вычислять в реальном времени оценки функциональной взаимосвязи участков коры головного мозга, соответствующих паре виртуальных отведений.

Протокол эксперимента Сценарий эксперимента, как правило, состоит из нескольких блоков. В соответствии с идеологией NFBLab настройки сигналов и стимула на протяжении блока остаются неизменными. По окончании каждого блока записанные данные добавляются в HDF5 файл с результатами и, как правило, запускается и обрабатывается одно или несколько событий из списка: обновление z-score статистик сигнала (среднее и стандартное отклонение или максимальное и минимальное значение), которые потом используются для стандартизации или нормировки сигнала обратной связи; запуск интерактивного модуля настройки фильтров по собранным данным для изменения или создания фильтров целевых сигналов на основе записанных ранее функциональных проб; приостановление эксперимента подача звукового сигнала сигнализирующего о завершении блока.

Блоки эксперимента различаются по типу визуализации. В данном параграфе описаны основные блоки - Baseline и Feedback. Блок Baseline заключается в предъявлении

текстового сообщения с инструкцией для испытуемого. Такие блоки как правило нужны для записи состояний испытуемого, например для записи фонового состояния, состояния с закрытыми глазами, моторных состояний, а также для сбора данных, используемых далее модулем настройки фильтров либо для обновления z-score статистик целевых сигналов. В блоке Feedback происходит визуализация и предъявление стимула НОС. Помимо предъявления реального НОС-сигнала в данном блоке есть возможность генерировать ложный сигнал из записанных ранее данных во время текущего эксперимента или из других экспериментов. Основным применением данной возможности является проведение экспериментов в контрольной группе. При необходимости типы визуализации можно расширять при помощи LSL интерфейса для подключения внешних визуализирующих программ к извлеченным целевым сигналам в режиме реального времени. При этом, однако, надо учитывать дополнительную задержку порядка десятков миллисекунд, вносимую LSL интерфейсом.

Сценарий эксперимента формируется из набора сконфигурированных блоков в виде последовательности (см. пример на рис. 7 приложение A). Некоторое подмножество блоков (или все блоки) может быть добавлено в группу блоков (Blocks group), внутри которой есть возможность повторять и перемешивать блоки. В конце каждого из блоков есть возможность реализовывать одно или несколько событий, описанных в предыдущем разделе. Примеры экспериментов и xml-скрипты можно найти в разделе 8 в приложении A.

2.2 Разработка низко-латентных методов оценки параметров ритмической активности мозга в режиме реального времени

Постановка задачи Важной характеристикой методов обработки сигналов в реальном времени является временное разрешение и задержка тракта обработки. Как было описано во введении низкая временная специфичность может быть причиной низкой эффективности реализации парадигм с замкнутым контуром. Настоящая работа направлена на решение проблемы временной специфичности. Современные программные решения, реализующие петлю обратной связи и применяемые как в клинике (BrainMaster, NeuroRT Training, Cygnet и т.д.) так и для исследований (OpenVibe, BCI2000) позволяют оценивать мощность осцилляторных активностей мозга с задержкой превышающей 500 мс. Данная задержка измеряется от момента получения данных с устройства ЭЭГ до момента передачи сигнала модулю визуализации стимула НОС. Дополнительная задержка порядка 100 мс появляется из-за технических причин, а именно вследствие организации передачи между ЭЭГ/МЭГ устройством и из-за временных затрат на генерацию стимула НОС исполнительным устройством, например монитором. Таким образом общая латентность

НОС системы как правило превышает 600 мс. Вследствие присутствия такой задержки стимуляция в парадигмах с замкнутом контуром может происходить в моменты, когда целевой паттерн активности уже прошел. Например, такой паттерн активности как короткий всплеск в альфа диапазоне (8-14 Гц) длится около 200-300 мс. Для детекции такого рода активности необходимо снизить задержку системы с замкнутым контуром как минимум до 100-200 мс. Данное исследование, проведенное в рамках диссертационной работы, посвящено разработке методов оценки мгновенной амплитуды и фазы узкополосных сигналов по регистрируемой в режиме реального времени ЭЭГ/МЭГ. Подробные результаты работы опубликованы в [19] и приведены в приложении B.

Со стороны компьютера, производящего прием и обработку ЭЭГ, сигнал, поступающий с электроэнцефалографа является многоканальным временным рядом с заданной частотой дискретизации, например равной fs = 500 Гц. Первым этапом обработки ЭЭГ является приведение многоканального сигнала к одноканальному виду посредством пространственной фильтрация. Без ограничения общности будем предполагать, что на вход разрабатываемому алгоритму поступает одноканальный сигнал х[п]. Далее одноканаль-ный сигнал х[п] может быть представлен в виде суммы двух сигналов:

х[п] = s[n] + rq[a]

, где s[n] целевой узкополосный сигнал, мощность и фазу которого необходимо оценить, и r¡[n] - широкополосный шум, влияние которого на оценку интересующей мощности требуется минимизировать. В такой модели предпологается что все нерелевантные узкополосные источники отфильтрованы заранее с пренебрежимо малым искажением сигнала в целевой области частот. Например, источники шума 50Hz или 60Hz связанные с наводками от электрической сети могут быть отфильтрованы режекторными либо гребенчатыми фильтрами без значительных искажений сигнала в спектральном диапазоне, к примеру, альфа-ритма. Нерелевантные нейрональные источники, источники связанные с работой мышц и общий фоновый ЭЭГ шум со спектром вида 1// формируют шумовую широкополосную компоненту r¡[n] сигнала x[n].

Далее, s[n] может быть преобразован в комплекснозначный аналитический сигнал с помощью преобразования Гильберта. Полученный сигнал у[п] представим в виде:

у[п] = а[п]езф[п]

где а[п] - мгновенная амплитуда сигнала (квадратный корень из мгновенной мощности), ф[п] - мгновенная фаза, j - мнимая единица. Оценка а[п] и ф[п] при известных значениях у[п] получается вычислением абсолютного значения и угла комплексного числа у[п], а именно:

а[п] = (Ке(у[п])2 + 1т(у[п])2)2 ф[п] = агсЪд(у[а]/х[а])

, Ке(у[п\) - вещественная часть у[п], /т(у[и]) - мнимая часть у[п].

Следует заметить, что операции вычисления модуля и аргумента аналитического сигнала не вносят дополнительной фундаментальной задержки в тракт обработки, так как вычисляются для каждого момента времени с использованием значений сигнала только в этот момент времени. Однако, вычисление преобразования Гильберта в точке п' в идеале требует бесконечного окна, центрированного вокруг точки п' на временной оси. Воспользовавшись классической аппроксимацией фильтров с бесконечной импульсной характеристикой возможно представить вычисление преобразования Гильберта как свертки сигнала и конечной импульсной характеристики фильтра Гильберта. Однако такая аппроксимация предполагает некаузальную обработку, требуя для оценки преобразованного сигнала в момент времени п' знаний входных значений 5 [и] для и как из прошлого по отношению к текущему значению времени (п < п') так и значений из будущего (п > п'). Такое преобразование не может быть выполнено в режиме реального времени. Однако, применение такого алгоритма для известной полной записи ЭЭГ позволяет выделить огибающую а[п] и фазу ф[п], которые далее будут называться идеальной огибающей и идеальной фазой, за исключение значений на краях соответствующей записи.

Таким образом задача данного раздела формулируется как построение каузального алгоритма, оценивающего амплитуду а[п] и фазу ф[п] в режиме реального времени по одноканальному сигналу х[п]. При этом в качестве одного из параметров разрабатываемого метода должна выступать явно задаваемая задержка метода И. Наилучшим методом является тот метод, который позволяет с уменьшением параметра И производить оценку а[п] и ф[п] в режиме реального времени с наилучшим качеством, которое будет определено ниже.

Существующие методы Классическим методом (далее обозначается как тесЬ) выделения мгновенной мощности является метод, основанный на амплитудной демодуляции сигнала и схож с принципом работы простейшего радиоприемника, детектирующего амплитудно модулированный сигнал [9]. Данный метод включает в себя три последовательных шага: узкополосная фильтрация в заданном диапазоне, "выпрямление" узкополосного сигнала (вычисление абсолютного значения) и сглаживание сигнала фильтром низких частот (ФНЧ). Выходом данного алгоритма является оценка мгновенной амплитуды а[п]. Задержка данного алгоритма складывается из задержек узкополосного фильтра и фильтра низких частот. В случае, когда в качестве фильтров используются симметричные

фильтры с конечной импульсной характеристикой (КИХ) данная задержка составляет половину суммы длин импульсных характеристик фильтров. В данной работе, в качестве фильтров используются фильтры с КИХ, а в качестве параметров используются длина М1 импульсной характеристики узкополосного фильтра и величина задержки И. Соответственно, длина импульсной характеристики ФНЧ может определяется из значений желаемой полной задержки И и длины КИХ первого фильтра N1 как:

Список литературы диссертационного исследования кандидат наук Сметанин Николай Михайлович, 2021 год

Список литературы

[1] J. Wolpaw. Brain-computer interfaces: principles and practice. Oxford Univ. Press, 2012.

[2] R. Sitaram, T. Ros, L. Stoeckel, S. Haller, F. Scharnowski, J. Lewis-Peacock, N. Weiskopf, M. Blefari, M. Rana, E. Oblak, N. Birbaumer, and J. Sulzer. Closed-loop brain training: the science of neurofeedback. Nature Reviews Neuroscience, 2016.

[3] Ahmed El Hady. Closed Loop Neuroscience. Elsevier, 2016.

[4] C. Zrenner, P. Belardinelli, F. Muller-Dahlhaus, and U. Ziemann. Closed-loop neuroscience and non-invasive brain stimulation: A tale of two loops. Front Cell Neurosci., 2016.

[5] Mark Lawrence Andermann, Jaakko Kauramaki, Tapio Palomaki, Christopher I. Moore, Riitta Hari, Iiro P. Jaaskelainen, and Mikko Sams. Brain state-triggered stimulus delivery: An efficient tool for probing ongoing brain activity. Open journal of neuroscience, 2, 2012.

[6] M. Congedo, J.F. Lubar, and D. Joffe. Low-resolution electromagnetic tomography neurofeedback. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2004.

[7] D.J. White, M. Congedo, and J. Ciorciari. Source-based neurofeedback methods using EEG recordings: training altered brain activity in a functional brain source derived from blind source separation. Front. Behav. Neurosci, 2014.

[8] Buzsaki G. and Watson B. Brain rhythms and neural syntax: implications for efficient coding of cognitive content and neuropsychiatric disease. Dialogues Clin Neurosci, 2012.

[9] Oppenheim A.V. and Schafer R.W. Discrete-Time Signal Processing. 3rd Edition. Prentice Hall, 2010.

[10] G Pfurtscheller and A Aranibar. Evaluation of event-related desynchronization (ERD) preceding and following voluntary self-paced movement. Electroencephalogr Clin Neurophysiol., 46(2):138-46, 1979.

[11] Dennis Gabor. Theory of communication. part 1: The analysis of information. Journal of the Institution of Electrical Engineers-Part III: Radio and Communication Engineering, 93(26):429-441, 1946.

[12] N. Evans, S. Gale, A. Schurger, and O. Blanke. Visual feedback dominates the sense of agency for brain-machine actions. PLoS ONE, 2015.

[13] E.F. Oblak, J.A. Lewis-Peacock, and J.S. Sulzer. Self-regulation strategy, feedback timing and hemodynamic properties modulate learning in a simulated fmri neurofeedback environment. PLoS Comput Biol, 2017.

[14] Anastasiia Belinskaia, Nikolai Smetanin, Mikhail Lebedev, and Alexei Ossadtchi. Short-delay neurofeedback facilitates training of the parietal alpha rhythm. Journal of Neural Engineering, 17(6):066012, December 2020.

[15] Chen L.L., Madhavan R., Rapoport B.I., and Anderson W.S. Real-time brain oscillation detection and phase-locked stimulation using autoregressive spectral estimation and time-series forward prediction. IEEE Trans Biomed Eng, 60:753-762, 2013.

[16] J R Mcintosh and P Sajda. Estimation of phase in EEG rhythms for real-time applications. Journal of Neural Engineering, 17(3):034002, June 2020.

[17] Alexandre Gramfort, Martin Luessi, Eric Larson, Denis Engemann, Daniel Strohmeier, Christian Brodbeck, Roman Goj, Mainak Jas, Teon Brooks, Lauri Parkkonen, and Matti Hamalainen. Meg and eeg data analysis with mne-python. Frontiers in Neuroscience, 7:267, 2013.

[18] N. Smetanin. NFBLab software. https://github.com/nikolaims/nfb. GitHub, 2020.

[19] Nikolai Smetanin, Anastasia Belinskaya, Mikhail Lebedev, and Alexei Ossadtchi. Digital filters for low-latency quantification of brain rhythms in real time. Journal of Neural Engineering, 17(4):046022, August 2020.

[20] G.R. Grice. The relation of secondary reinforcement to delayed reward in visual discrimination learning. Journal of experimental psychology, 1948.

[21] H. Rahmandad, N. Repenning, and J. Sterman. Effects of feedback delay on learning. System Dynamics Review, 2009.

[22] A. Ossadtchi, T. Shamaeva, E. Okorokova, V. Moiseeva, and M. Lebedev. Neurofeedback learning modifies the incidence rate of alpha spindles, but not their duration and amplitude. Scientific Reports, 2017.

[23] H. Shin, R. Law, S. Tsutsui, C. Moore, and S. Jones. The rate of transient beta frequency events predicts behavior across tasks and species. eLife, 2017.

[24] N. Smetanin and A. Ossadtchi. Express estimation of brain rhythm power for low-latency neurofeedback. The First Biannual Neuroadaptive Technology Conference, 2017.

[25] Nikolai Smetanin, Ksenia Volkova, Stanislav Zabodaev, Mikhail A. Lebedev, and Alexei Ossadtchi. NFBLab—a versatile software for neurofeedback and brain-computer interface research. Frontiers in Neuroinformatics, 12, December 2018.

[26] Swartz Center for Computational Neuroscience. Lab streaming layer. https://github.com/sccn/labstreaminglayer. GitHub, 2020.

[27] Y. Renard, F. Lotte, G. Gibert, M. Congedo, E. Maby, V. Delannoy, O. Bertrand, and A. Lecuyer. Openvibe: An open-source software platform to design, test and use brain-computer interfaces in real and virtual environments. Presence : teleoperators and virtual environments, 2010.

[28] G. Schalk, D.J. McFarland, T. Hinterberger, N. Birbaumer, and J.R. Wolpaw. BCI2000: A general-purpose brain-computer interface (BCI) system. IEEE Trans Biomed Eng, 2004.

[29] A.J. Bell and T.J. Sejnowski. An information-maximization approach to blind separation and blind deconvolution. Neural Comput., 1995.

[30] Z.J. Koles, M.S. Lazaret, and S.Z. Zhou. Spatial patterns underlying population differences in the background EEG. Brain topography, 1990.

[31] V.V. Nikulin, G. Nolte, and G. Curio. A novel method for reliable and fast extraction of neuronal EEG/MEG oscillations on the basis of spatio-spectral decomposition. Neuroimage, 2011.

[32] Freek van Ede, Andrew J. Quinn, Mark W. Woolrich, and Anna C. Nobre. Neural Oscillations: Sustained Rhythms or Transient Burst-Events? Trends in Neurosciences, 41(7):415-417, July 2018.

[33] Natalie Schaworonkow, Pedro Caldana Gordon, Paolo Belardinelli, Ulf Ziemann, Til Ole Bergmann, and Christoph Zrenner. ^-rhythm extracted with personalized eeg filters correlates with corticospinal excitability in real-time phase-triggered eeg-tms. Frontiers in Neuroscience, 12:954, 2018.

[34] Til Ole Bergmann, Anke Karabanov, Gesa Hartwigsen, Axel Thielscher, and Hartwig Roman Siebner. Combining non-invasive transcranial brain stimulation with neuroimaging and electrophysiology: Current approaches and future perspectives. Neuroimage, 140:4-19, oct 2016.

[35] Aqsa Shakeel, Toshihisa Tanaka, and Keiichi Kitajo. Time-series prediction of the oscillatory phase of eeg signals using the least mean square algorithm-based ar model. Applied Sciences, 10(10), 2020.

Приложения

А Статья "NFBLab—A Versatile Software for Neurofeedback and Brain-Computer Interface Research"

Authors. Nikolai Smetanin, Ksenia Volkova, Stanislav Zabodaev, Mikhail A. Lebedev and Alexei Ossadtchi

Abstract. Neurofeedback (NFB) is a real-time paradigm, where subjects learn to volitionally modulate their own brain activity recorded with electroencephalographs (EEG), magneto-encephalographic (MEG) or other functional brain imaging techniques and presented to them via one of sensory modalities: visual, auditory or tactile. NFB has been proposed as an approach to treat neurological conditions and augment brain functions. Although the early NFB studies date back nearly six decades ago, there is still much debate regarding the efficiency of this approach and the ways it should be implemented. Partly, the existing controversy is due to suboptimal conditions under which the NFB training is undertaken. Therefore, new experimental tools attempting to provide optimal or close to optimal training conditions are needed to further exploration of NFB paradigms and comparison of their effects across subjects and training days. To this end, we have developed open-source NFBLab, a versatile, Python-based software for conducting NFB experiments with completely reproducible paradigms and low-latency feedback presentation. Complex experimental protocols can be configured using the GUI and saved in NFBLab's internal XML-based language that describes signal processing tracts, experimental blocks and sequences including randomization of experimental blocks. NFBLab implements interactive modules that enable individualized EEG/MEG signal processing tracts specification using spatial and temporal filters for feature selection and artifacts removal. NFBLab supports direct interfacing to MNE-Python software to facilitate source-space NFB based on individual head models and properly tailored individual inverse solvers. In addition to the standard algorithms for extraction of brain rhythms dynamics from EEG and MEG data, NFBLab implements several novel in-house signal processing algorithms that afford significant reduction in latency of feedback presentation and may potentially improve training effects. The software also supports several standard BCI paradigms. To interface with external data acquisition devices NFBLab employs Lab Streaming Layer protocol supported by the majority of EEG vendors. MEG devices are interfaced through the Fieldtrip buffer.

frontiers

in Neuroinformatics

METHODS

published: 24 December 2018 doi: 10.3389/fninf.2018.00100

®

Check for updates

OPEN ACCESS

Edited by:

Antonio Fernandez-Caballero, University of Castilla La Mancha, Spain

Reviewed by:

Alexander K. Kozlov, Royal Institute of Technology, Sweden Pyatin Vasiliy Fedorovitch, Samara State Medical University, Russia

Correspondence:

Alexei Ossadtchi [e-mail]

Received: 14 September 2018 Accepted: 12 December 2018 Published: 24 December 2018

Citation:

Smetanin N, Volkova K, Zabodaev S, Lebedev MA and Ossadtchi A (2018) NFBLab—A Versatile Software for Neurofeedback and Brain-Computer Interface Research. Front. Neuroinform. 12:100. doi: 10.3389/fninf.2018.00100

NFBLab—A Versatile Software for Neurofeedback and Brain-Computer Interface Research

Nikolai Smetanin1, Ksenia Volkova1, Stanislav Zabodaev2, Mikhail A. Lebedev1 and Alexei Ossadtchi1*

1 Center for Bioelectric Interfaces, National Research University Higher School of Economics, Moscow, Russia,2 Medical Computer Systems, Zelenograd, Russia

Neurofeedback (NFB) is a real-time paradigm, where subjects learn to volitionally modulate their own brain activity recorded with electroencephalographic (EEG), magnetoencephalographic (MEG) or other functional brain imaging techniques and presented to them via one of sensory modalities: visual, auditory or tactile. NFB has been proposed as an approach to treat neurological conditions and augment brain functions. Although the early NFB studies date back nearly six decades ago, there is still much debate regarding the efficiency of this approach and the ways it should be implemented. Partly, the existing controversy is due to suboptimal conditions under which the NFB training is undertaken. Therefore, new experimental tools attempting to provide optimal or close to optimal training conditions are needed to further exploration of NFB paradigms and comparison of their effects across subjects and training days. To this end, we have developed open-source NFBLab, a versatile, Python-based software for conducting NFB experiments with completely reproducible paradigms and low-latency feedback presentation. Complex experimental protocols can be configured using the GUI and saved in NFBLab's internal XML-based language that describes signal processing tracts, experimental blocks and sequences including randomization of experimental blocks. NFBLab implements interactive modules that enable individualized EEG/MEG signal processing tracts specification using spatial and temporal filters for feature selection and artifacts removal. NFBLab supports direct interfacing to MNE-Python software to facilitate source-space NFB based on individual head models and properly tailored individual inverse solvers. In addition to the standard algorithms for extraction of brain rhythms dynamics from EEG and MEG data, NFBLab implements several novel in-house signal processing algorithms that afford significant reduction in latency of feedback presentation and may potentially improve training effects. The software also supports several standard BCI paradigms. To interface with external data acquisition devices NFBLab employs Lab Streaming Layer protocol supported by the majority of EEG vendors. MEG devices are interfaced through the Fieldtrip buffer.

Keywords: neurofeedback, low-latency, software, real-time EEG, brain-computer interface, flexible experiment design, LSL-protocol

1. INTRODUCTION

In the biological feedback (BF) paradigm, subjects gain the ability to monitor and control the physiological parameters of their own bodies, recorded with artificial sensors. In the case of neurofeedback (NFB), subjects learn to control the activity of their brains (Kamiya, 1969; Sterman et al., 1974; Sitaram et al., 2016; Ossadtchi et al., 2017). In NFB settings, neural activity is recorded with such methods as electroencephalography (EEG), magnetoencephalography (MEG), electrocorticography (ECoG), functional magnetic resonance imaging (fMRI) or multielectrode implants and converted into stimuli that subjects can perceive. NFB processing stages include recording of neural activity, extracting features of interest from the recorded signals, feature transformation, and delivering the feedback signal to the subject via one of the sensory modalities, most often visual, auditory or tactile. A very similar processing scheme is utilized in brain-computer interfaces (BCIs), the systems that connect the brain to external assistive devices (Wolpaw, 2012). In BCIs, the feedback is implemented in the form of response of the external device to the commands derived from the neural activity of the user; this signal can be considered as a type of NFB.

NFB and BCI paradigms are utilized in neurophysiological research, clinical practice, and consumer applications, like video games, meditation and performance training gadgets. In the clinical world, NFB methods often complement conventional approaches to treatment of neurological disorders (Coben et al., 2010; Lofthouse et al., 2012) and cognitive enhancement training (Zoefel et al., 2011). EEG-based BCIs can be used for post-stroke neurorehabilitation (Ang et al., 2015; Frolov et al., 2018), where motor intentions are extracted from brain activity and directed to exoskeleton devices or functional electrical stimulators that evoke limb movements, which in turn generate sensory re-afferent signals that flow back to the brain and mediate restorative brain plasticity.

EEG-based NFB systems have become widespread because EEG recordings are noninvasive, relatively simple to implement, and can be informative of many neural functions. In modern implementations, multichannel EEG signals are amplified, digitized, and transmitted to a computer for denoising and feature extraction. EEG features are often derived from oscillatory components of the multichannel EEG data, such as the theta (3-7 Hz), alpha (8-12 Hz), mu (9-14 Hz), and beta (15-25 Hz) rhythms. The neural mechanisms underlying different EEG features depend on the rhythm frequency and cortical source. NFB therapy usually attempts to correct a targeted physiological function by attempting to increase or decrease particular EEG rhythm power. In clinical settings, quantitative EEG (QEEG) technology is utilized to detect the difference between the EEG of a patient and database values for healthy controls. Many QEEG methods incorporate statistical analyses of the EEG spectral bands. Additionally, the group independent component analysis (gICA) has been developed, where EEG dynamics of a patient is compared to a database of independent components (Ponomarev etal., 2012).

The efficiency of NFB training critically depends on the appropriate selection of the spatial, spectral and temporal

characteristics of the feedback signal. To ensure that these characteristics reflect in a timely fashion the neural functions of interest, several spatial and temporal filtering methods have been proposed that isolate signal ofinterest from noise and concurrent brain activity present in multichannel EEG/ECoG/MEG data, with minimal possible latency. Spatial filtering is implemented based on the signal analysis by means of popular algorithms such as principal component analysis (PCA), common spatial pattern (CSP) (Koles et al., 1990), independent component analysis (ICA) (Bell and Sejnowski, 1995), and spatio-spectral decomposition (SSD) (Nikulin et al., 2011) etc. Temporal filtering is then applied to the extracted spatial component in order to extract the signal from specific frequency band and estimate its instantaneous power (Lee et al., 2006; Smetanin and Ossadtchi, 2017). Individual variability of EEG rhythms is an important factor that affects NFB efficiency. Ideally, each subject's anatomical and physiological features should be taken into consideration and NFB-generating algorithms adjusted accordingly. Such adjustment is, however, rarely done, which causes sub-optimal performance of NFB. Overall, there is still no established procedure that ensures a reliable and efficient NFB. Unimpressive performance of several NFB implementations, particularly the ones with double blind controls, has caused criticism toward the NFB paradigm as a whole (Sitaram et al., 2016).

Among the obstacles hindering the achievement of better results with NFB, are the limitations of the existing software tools, which in many cases are cumbersome to use, are not flexible and do not enable certain desirable features. To address this range of problems, here we report a new software, the NFBLab, suitable for conducting many NFB and BCI experiments. The platform incorporates both standard and novel algorithms for real-time processing of EEG/MEG/ECoG activity and generating NFB. The key features of this software include:

• Support of the most common EEG/MEG devices using the lab streaming layer (LSL)1 and FieldTrip buffer2;

• Internal language that allows implementing different experimental protocols using XML-formatted descriptions;

• Graphical user interface for design of signal processing tracts and to form flexible sequences of experimental blocks;

• Interactive module for configuring signal processing pipelines based on the test data collected in individual subjects;

• Support of EEG/MEG inverse problem solvers (e.g., wMNE and LCMV beamformers) using direct interaction with MNE-Python software;

• Continuous visualization of signal features;

• Support of mock feedback and experimental blocks randomization;

• New in-house algorithms for low-latency processing of narrow-band signals (Smetanin and Ossadtchi, 2017);

• Open-source Python code3 and cross-platform compatibility.

1 https://github.com/sccn/labstreaminglayer

2 http://www.fieldtriptoolbox.org/development/realtime/buffer

3 https://github.com/nikolaims/nfb

2. COMPARABLE SOFTWARE

Several software platforms are currently available that furnish flexibility needed for conducting experiments that involve realtime processing and visualization of multichannel bioelectric signals. Among them, OpenVIBE (Renard et al., 2010) is a very popular, actively supported project. It incorporates visualprogramming tools and scripts for configuring complex signal processing pipelines. An experimental task is usually subdivided into several parts specified in separate XML files (e.g., a motor-imagery BCI4). Each part is launched manually, which may adversely affect reproducibility of the conducted experiments. Additionally, this software is not open-source.

BCI2000 is another popular software solution for implementing NFB and BCI paradigms (Schalk et al., 2004). This is a closed-source C/C++ software with several built-in algorithms for processing neural activity and extracting signals of interest. Like OpenVIBE, different experimental tasks are launched manually in BCI2000. BCI2000 supports integration with external programs and MATLAB environment, which allows users to incorporate additional data processing modules.

The NFBLab software proposed here advances the functionality of both OpenVIBE and BCI2000 platforms. One innovation is the NFBLab's capacity to configure signal processing pipelines based on the data collected in individual subjects. This functionality incorporates a rich set of decompositions (PCA, ICA, CSP, SSD) for extracting essential neural modulations while discarding noise. Furthermore, the NFBLab software contains an advanced experiment control module that enables tasks consisting of multiple blocks. The blocks can be automatically launched, switched, repeated and randomized. Mock NFB is supported, as well. NFBLab is written in Python. The code is open-source, which facilitates sharing and improving these software tools by the developers of NFB and BCI systems. The code and installation instructions can be found at https://github.com/nikolaims/nfb repository.

In addition to Open VIBE, BCI2000 and our NFBLab platform, several commercial clinical software platforms should be mentioned that have been designed for NFB therapy. These are BrainMaster, NeuroRT Training, Cygnet and several others. These systems should be clearly delineated from the research oriented software mentioned earlier. They operate under a very harsh requirement for exclusively high usability within the clinical environment which often adversely affects their flexibility and the fidelity of training protocols. While some of these platforms support spatial filtering of neural signals, the parameters of such filters are obtained from solving the generic inverse problem without taking into account the anatomical and functional parameters of individual subjects. Thus, sLORETA (Pascual-Marqui, 2002) with a generic head model is commonly utilized within popular NFB software. With this approach, the output signal could be different from the individual subject's true activity of the brain region of interest (ROI). Additional inaccuracies are introduced to the performance of forward models by the unknown tissue conductivity profiles.

4http://openvibe.inria.fr/motor-imagery-bci/

Furthermore, NFB systems suffer from the presence of delay between neural modulations and NFB delivery. This delay is typically more than 500 ms, which could be too long for NFB to be efficient (Larsen and Sherlin, 2013).

NFBLab offers solutions to both spatial and temporal NFB specificity problems. In the spatial domain, there is an option to use a subject-specific spatial filter. The filter parameters are derived based on the spatial-temporal decomposition (PCA, ICA, CSP, SSD) applied to the EEG/MEG/ECoG data recorded during functional tests, for example opening/closing the eyes, blinking, and performing hand movements. These functional tests allow for building subject-specific spatial filters without the need to perform laborious calculations, such as constructing forward models from the individual MRI data. Additionally, NFBLab supports conventional ROI-based approaches by importing, via the MNE-Python package, inverse operators calculated using a broad range of methods (Darvas et al., 2004) based on the generic or individual anatomy.

In the temporal domain, NFBLab incorporates several advanced signal processing algorithms that effectively decrease NFB latency. These algorithms process EEG signals and extract envelopes of their narrow-band components with minimal possible or even negative latency.

3. ARCHITECTURE

The NFBLab platform consists of three main modules, as illustrated in Figure 1. The first module, called "Experiment protocol editor," specifies an experimental task. It describes the sequence of experimental blocks, signal processing procedures to compute target signals and the parameters used to calculate NFB within each experimental block. The block descriptions are saved as XML files, which are then loaded into the second module, called "Experiment module," which needs to be launched at the start of the experiment. This XML-language based structure allows for reproducibility of experiments. Experimental module processes and displays the raw and transformed neural signals. The Experiment module also controls the sequence of experimental blocks and the presentation of stimuli and NFB. The third module, Data-driven filter designer, is an interactive module for editing signal processing flow and constructing spatial and temporal filters. The filters are based on the frequency analysis and spatial decomposition of the EEG data obtained from the functional tests. Short latency data exchange with EEG/MEG/ECoG recording devices is mediated by the LSL or FieldTrip buffer technologies. The experimental data along with the synchronized feedback signal is saved as HDF5 files, and output signals are streamed to the LSL outlet for communication with the external programs and devices. The relationships between NFBLab modules is shown in Figure 1 .

3.1. Support of External Devices

The main instrument for working with external devices in NFBLab is the Lab Streaming Layer (LSL) protocol. LSL is a socket based data transfer protocol that controls synchronized

FIGURE 1 | NFBLab architecture. Three main modules comprise NFBLab. Experiment protocol editor details experimental task. It describes the sequence of experimental blocks, properties of the individual blocks, sequence of blocks within the experiment, signal processing procedures used NFB presentation parameters. Data-driven filter designer allows to interactively edit signal processing flow and construct spatial and temporal filters. The filters are based on the frequency analysis and spatial decomposition of the EEG data obtained from the functional tests. Experiment module does actual real-time processing and displays the raw and transformed neural signals. It also controls the sequence of experimental blocks and the presentation of stimuli and NFB.

collection and transmission of multi-channel time series. The protocol can also add meta-information, for example, names of devices and recording channels. The central element of LSL is a multichannel data stream whose input/output operations are supported by LSL Inlet and LSL Outlet objects, correspondingly.

Data streaming from the EEG/MEG/ECoG devices to the NFBLab is carried out through the LSL protocol typically implemented within a separate independent module that reads data following the specifications of the specific acquisition device and converts it into a standard LSL stream. Acquisition software of some manufacturers of EEG recording devices includes LSL streaming feature, for example NeoRec software (Medical Computer Systems, Ltd.) and EEG Studio (Mitsar Co., Ltd). NFBlab intercepts such streams by name and allows selecting/excluding channels, and defining a new EEG reference.

NFBLab supports multiple LSL streams running simultaneously, for example simultaneous data acquisition and processing with several external devices. These are typically EEG/ECoG/MEG recording systems, but other types of data can be streamed, as well, such as multichannel electromyography, thermometry, eye-tracking and limb kinematics acquired using motion capture systems.

The NFBLab interface allows to set up a broadcast of input and output signals, and their visualization. For test purposes, previously recorded data could be replayed as an LSL stream at a specified sampling rate, which allows designing experimental tasks without the need to connect to any external equipment. In addition to the LSL streams, NFBLab supports FieldTrip buffer, an alternative data transfer protocol, which allows for nearly real-time communication with magnetoencephalographic acquisition devices.

4. SIGNAL PROCESSING PIPELINE SPECIFICATION

NFBLab implements virtual derivations (or virtual leads) and applies narrow-band filtering followed by envelope extraction to obtain derived signals representing instantaneous power within specific EEG band. It also allows for mathematical operations on pairs of derived signals to obtain composite signals. Appropriately scaled derived and composite signals can be then presented to a subject as NFB. Figure 2 schematically illustrates this processing structure.

4.1. Virtual Derivations and Derived Signals

Here and below, virtual derivation (virtual lead) is defined as an instantaneous linear combination of the signals from different channels stored in vector x(t), which formally can be written as y(t) = wTx(t). Although it is not exactly technically correct, in the EEG signal processing community this operation is referred to as spatial filtering. The weight vector w of this linear combination can be trivial, that is it can consist of all zeros except for "1" in a single position. In this case, the virtual derivation is simply the signal from the channel corresponding to the index of the non-zero element of the weight vector.

In the source-space NFB paradigm (Congedo et al., 2004), NFB represents activity of a cortical ROI. This signal is usually calculated by solving the EEG/MEG inverse problem, where an M x N inverse operator matrix H = {hij}, (i,j) e [(1, M), (1, N)] is calculated with the i—th row hT corresponding to the i—th ROI. The ROI activity yi(t) is derived from the N-channel signal x(t) as a dot product yi(t) = hTx(t). Thus, the spatial filter vector wT = hT . While there exist a plethora of methods for solving the EEG/MEG inverse problem (Darvas et al., 2004), most of the implementations of source space neurofeedback

FIGURE 2 | Derived and composite signals processing pipeline. Virtual lead Is obtained by means of spatial filtering operation applied to EEG channels, this corresponds to computing a linear combination of EEG channel signals with coefficients aimed to emphasize features of interest and suppress the interference. Derived signal usually represents envelope of the narrow-band filtered virtual lead signal. Composite signal can be computed by applying an arbitrary user-specified mathematical function to a pair of derived signals. Coherence represents a special case when the composite signal is obtained based on the pair of virtual signals (not derived signals).

rely on the generic sLoreta (Pascual-Marqui, 2002) solver based on the standardized head model. NFBLab implements sLoreta-based neurofeedback, as well, but also provides an interface to a full-blown inverse problem solver implemented in MNE-Python software, including several versions of the Minimum Norm Estimates (Hamalainen and Ilmoniemi, 1994) and beamforming approaches (Greenblatt et al., 2005).

The accuracy of the inverse modeling approaches suffers from the individual variability in head-shapes and spatial profiles of head tissue conductivity. The conductivity effects are especially strong when EEG-based NFB is used because tissue conductivity significantly affects the electrical potentials recorded from the scalp. To circumvent this problem, NFBLab offers an alternative approach for building spatial filters. This approach is based on the use of functional tests. One version of this approach was previously reported (White et al., 2014). In a general case, this approach consists of collecting segments of EEG/ECoG/MEG data during functional tasks, such as eyes opening/closing, rotating the hand, blinking, relaxing and several others. These tasks yield EEG/ECoG/MEG data that can be used for identifying patterns of interest, which are localized in space-time-frequency domain and reflect functionally relevant modulations of the underlying neuronal activity. Meaningful task/state related components are decomposed from the multichannel data using mathematical signal processing techniques, such as ICA, PCA, SSD, and CSP. Each component Zi(t) is obtained from the channel signals x(t) by applying a spatial filter w;: Z;(t) = wTx(t). This way, functionally relevant NFB can be generated without the need to solve the complete inverse problem. Individualized spatial filters constructed by NFBLab reflect activity of specific neuronal populations corresponding to the physiological functions of interest. Thus, for example, using the spatial filter obtained by means of the common spatial patterns (CSP) analysis contrasting the eyes-open and eyes-closed conditions, one can derive a spatial filter for extracting activity of the occipital alpha-rhythm generators. A spatial filter for extraction of the sensory-motor rhythm can be obtained using SSD analysis of the data collected during motor relaxation.

Spatial filtering can be used not only to extract functionally relevant signals but also to suppress the irrelevant signals ubiquitously present in the neural recordings. These are artifacts resulting from the heartbeat, eye blinks and other sources, and/or neural modulations unrelated to the function of interest.

The irrelevant signals are found after computing a spatial decomposition of the EEG/MEG data and examining the component time-series z(t) = [z1(t),z2(t),...,zn(t)]T = Bx(t), where B is spatial decomposition matrix (e.g., unmixing matrix in case of ICA). For example, the k-th component reflects an artifact and needs to be removed. This could be done by simply zeroing the k-th component Zk(t) = 0 and reconstructing the data from the rest of the components. This is equivalent to excluding the k-th row in the analysis matrix B and removing the k-th column of matrix C = B-1. The denoised data can be represented as xc(t) = B-kB_kx(t) where the subscript "—k" indicates the removal of the k-th column and the k-th row from B—1 and B, respectively. Such data cleaning can be combined with the spatial filtering operation described above. The spatial filter takes the form W; = (B—^ k B—kW; and then the interference-free virtual derivation signal can be computed as y;(t) = WTx(t).

4.2. Narrow-Band Component Envelope

Many NFB applications are based on the estimation of dynamical changes in oscillatory brain activity, such as EEG rhythms in specific frequency bands. Brain rhythms of interest are extracted from the ongoing EEG/MEG activity using digital band-pass filters. Filter parameters are calculated based on the frequency range and filter order specified by the user. Generally, rhythmic brain activity is non-stationary and can be quantified as the instantaneous power of the oscillatory component envelope. In what follows we will denote as d;(t) the envelope of the narrowband virtual derivation signaly;(t). Signals d;(t) are then referred to as derived signals.

4.2.1. Basic Methods

The existing NFB implementations use two major methods to estimate the envelope. The first method consists of performing low-pass filtering of the rectified narrow-band signal. This method is represented by the diagram in Figure 3A. This simple and ubiquitously used approach, however, introduces considerable delays that make it difficult for the subject to utilize the resultant NFB as it often arrives when the reinforced activity pattern has been already abandoned. Consequently, the desired pattern of brain activity cannot be efficiently reinforced. Attempts to reduce the delay using lower-order filters (including their minimum phase versions) result in a rapid deterioration of performance as shown in Figure 4 by the dashed black curve.

FIGURE 3 | Techniques that can be used to estimate the envelope of a narrow-band signal. (A) Narrow-band filtering with an Infinite impulse response (IIR) filter followed by the rectification results in a significant delay. (B) STFT with subsequent temporal smoothing of the absolute values of the STFT coefficient. Combined with mirror-reflection, this technique allows to slightly reduce the delay as compared to (A). (C) Complex demodulation based technique that has the lowest delay and allows to obtain narrow-band envelope estimation with the latency of 100 ms. The typical delay values reported for each of the techniques correspond to the combination of parameters that deliver comparable accuracy of narrow-band envelope estimation with all the three methods.

The second commonly practiced way to estimate instantaneous band power modulation is illustrated in Figure 3B. This method is based on the use of the short-time Fourier transform (STFT) with subsequent temporal smoothing of the absolute values of the STFT coefficients. We suggest that the temporal performance of the STFT-based technique can be improved with the approach where the data segment is mirror-reflected against the current time, followed by windowing and Fourier transform. Therefore, the FFT is applied to the symmetric sequence created from the segment of the last T samples reflected around the last acquired time-sample, which yields Xt—T,Xt—t+i,...,Xt,...,Xt—t+1,Xt—T as a sequence to perform the FFT on. This way, the subsequent application of a symmetric scaling window (Hann, Hanning, Blackman, etc.) effectively emphasizes contribution of the most recent samples. This method, on average, shortens the delay to around 250 ms while retaining sufficient accuracy of the envelope profile estimation. As it is the case with the rectification-based method, further attempts to reduce the delay lead to a significant drop in estimation accuracy, as illustrated by the dot-dashed black curve in Figure 4A.

4.2.2. Advanced Envelope Estimation Techniques

In addition to these two methods, NFBLab implements several advanced algorithms for decreasing the delay. An optimal latency-accuracy trade-off can be achieved with these methods. The general idea behind these techniques is based on the use of

a finite impulse response filter (FIR) approximation of narrowband Hilbert transform. This approach is schematically presented in Figure 3C. The accuracy of envelope estimation is inversely proportional to the length of the analysis interval while the longer intervals lead to longer delays. In contrast to the traditional methods, these advanced techniques allow to explicitly specify the desired delay value and then achieve the best possible accuracy for the specified delay.

In the Hilbert transform-based methods, narrow-band signal s[n] is represented as the real part of the analytic signal:

y[n] = s[n] + jsh[n] = A[n]ej(wc"+^[n])

where sh[n] is the imaginary part of the analytic signal, often called "second quadrature" of the original signal s[n], wc is the central band frequency, A[n] is the instantaneous power of the narrowband process, and <p[n]—is the instantaneous phase. Thus, the envelope power A[n] can be computed as the squared norm of the analytic signal.

To build the analytic signal y[n] corresponding to the original narrow-band signal s[n] derived from the noisy broad-band signal x[n], one can apply a narrow-band Hilbert transform filter. Frequency response of the filter can be defined as:

H (jw) í 2e-jwd, w € [wc - Sw, wc + Sw] c [0, n ] (1)

Hd(e ) =1 0, otherwise (1)

where Sw is the half of the bandwidth and d is the group delay in the samples. For any finite d this filter is non-causal and cannot be

applied in real time. To reconstruct analytical signal causally, one can use a complex-valued finite impulse response filter (cFIR) that approximates frequency response Hd (ejw). The desired filter b can be found by solving the least squares optimization problem defined in different ways.

The first and the most straightforward approach (denoted F-cFIR) that implements the least squares filter design strategy is to find the complex valued vector of cFIR filter weights b that minimizes the L2 distance between the cFIR filter frequency response and the ideal response Hd in the frequency domain:

bF = argmin ||Fb - Hd||l2 b

where F is the discrete Fourier transform operator. This simple approach does not take into account the EEG temporal structure and can be further improved.

One improvement can be achieved using optimal filter design ideas. In this approach, spectral density X of the input signal x[n] provides weights. We thus formulate the weighted frequency domain least squares design technique (denoted WF-cFIR) as the following optimization problem :

bwF = arg min || (Fb - Hd) • X||l2 b

This method allows to exploit the patterns and hidden relationships between rhythmic components in the EEG signal.

The last approach from this family (denoted T-cFIR) is based on minimization of the squared distance in the time domain between the complex delayed ground truth signal y[n — d] and the filtered signal x[n] * b[n]:

by = arg min ||x[n] * b[n] — y[n — d] | |l2 b

where * is the convolution operator. Here, during the training stage, the ground truth signal yd is obtained non-causally from the training dataset via ideal zero-phase Hilbert transformer (1). According to Parseval's theorem, this approach is equivalent to the WF-cFIR approach. However, in contrast to the frequency domain formulation, this method uses a straightforward way to add amplitude dependent weights in order to achieve better accuracy of the envelope shape estimation.

As illustrated in Figures 4A,C the advanced techniques described above allow for envelope estimation with a significant improvement in latency and accuracy as compared to standard methods (colored vs. black lines) used in the majority of the NFB software. When applied to a white noise sequence, see Figure 4B, the performance of envelope estimation methods tends to decrease.

It has been recently emphasized that both beta-rhythm (Shin et al., 2017) and alpha-rhythm activity (Ossadtchi et al., 2017) can be described as a sequence of oscillatory episodes, i.e., discrete events rather than continuous modulations of these rhythms. During alpha NFB training, the duration of each oscillatory episode stays constant, and the training effect consists of the increase in the number of episodes per unit time (Ossadtchi et al., 2017). These events (e.g., beta or alpha episodes) have

characteristic duration of 300 ms or shorter. Given such a short duration, the decreased NFB latency achieved with the advanced techniques will result in reinforcement arriving on time, during the rhythmic episode and not when it has finished (Oblak et al., 2017).

4.3. Normalization and Standardization

It is customary for NFB applications to apply various normalization methods prior to presenting NFB to the subject (Enriquez-Geppert et al., 2017). In NFBLab, derived and composite signals are standardized based on their statistics estimated for a baseline data segment. Standardization consists of subtracting the mean of the derived/composite signal y(t)

y(t)_^^

and division by the standard deviation: y(t) = —y. This operation ensures that y(t) is centered around zero with maximal amplitude around 3. The use of the standardized signal allows to easily control the fraction of outliers and adjust the feedback threshold (Thatcher and Lubar, 2008). Additionally, median based statistics can be used for robust normalization of outliers.

Some applications require operations with only positive signals. In these cases, robust estimates of ymin and ymax statistics can be used, where normalized signal is computed as y(t) = y(t)—ymir

ymax ymin

The normalized signal changes in the [0,1] range.

4.4. Motor Imagery BCI Paradigm Support

The derived signals-based structure fits well into the processing framework of modern motor-imagery BCIs. Accordingly, a separate module of NFBLab software offers a solution for generating the derived signals based on the output of an EEG classifier tuned to distinguish motor-imagery commands. This processing pipeline includes: (1) pre-filtering of EEG signal in some relatively broad frequency range, e.g., 0.5-45 Hz, (2) followed by the band-pass filtering of data in several overlapping bands (e.g., 6-10 Hz, 8-12 Hz, ..., 20-24 Hz), (3) spatial decomposition, e.g., CSP aimed at contrasting the data recorded within different motor-imagery states; (4) estimation of the instantaneous power in each of the spatial components isolated at step 3 and smoothing with a moving average (MA) with a time constant of 0.5 s; (5) application of a motor-imagery classifier (e.g., a classifier for delineating right and left hand imagery and resting conditions). Stages 3 and 5 include parameter adjustments for individual users; this adjustment requires an EEG sample recorded during the imagery states of interest.

5. MAIN NFBLAB OBJECTS AND THEIR DESCRIPTORS

5.1. Derived Signal

According to the NFBLab terminology, derived signal (DS) is the envelope of the narrow-band filtered virtual lead signal. In most cases, the DSs are scalars varying over time. The number of DSs that a user can define is limited only by the performance characteristics of the computer running the NFBLab. The DSs are described in a specific XML structure. Figure 5A shows the dialog for editing DS settings. This dialog can be called from the "Experiment design" module before launching an experiment.

FIGURE 4 | Correlation coefficient for the envelope obtained causally by each of classic (black curves) and novel (colored curves) filters with the ground-truth sequence vs. filters group delay for real EEG P4-channel data (A) and white noise signal (B). Example of T-cFIR envelope reconstruction for 0, 100, and 200 ms delays (C).

Additionally, DSs can be set and/or changed based on the data from the functional tasks. For this purpose, an interactive Data-driven filter designer module can be called during an experiment. This tool allows to choose the desired frequency band, perform spatial decomposition of the data and choose appropriate spatial filters. In addition to the interactive dialogs, advanced users can edit the XML file directly. The DSs can be linked to the NFB presentation module so that NFB output represents one of the DSs or a composite signal formed from several DSs.

In the example shown in Figure 5 we have defined a Derived Signal named Alpha to be obtained from the raw EEG data by spatial filter coefficients stored in OccipitAlphaFilter.txt and with 4-th order Butterworth signal with a passband of 911 Hz. Envelope smoothing is performed with Savitsky-Golay polynomial filtering.

5.2. Composite Signal

NFB often incorporates several signals, such as different EEG rhythms. For example, alpha-to-theta ratio feedback is based on the instantaneous power of the alpha band divided by the power of the beta band. For such computations, NFBLab implements the composite signal class, which is defined as an arbitrary mathematical function of two DSs. As shown in Figure 5B, this function is specified by the user. For example, to define a training protocol which requires calculation of the ratio of the occipital alpha

rhythm power to the frontal theta rhythm power, the DSs representing the corresponding bands are used, and the division function is specified. NFBLab also allows for more complex computations, such as assessment of correlation between virtual lead signals using such metrics as coherence or envelope correlation.

In addition to forming NFB, composite signals can be used to track the presence of artifacts in neural recordings, such as eye blinks, muscular artifacts, and mechanical artifacts caused by movements. A derived or composite signal representing artifact can be converted into a warning stimulus instructing the subject to correct the undesired behavior.

5.3. Experimental Block

NFB experiments usually consist of several blocks. In NFBLab, a block is defined as the task part where the stimulus settings and signal processing procedures remain unchanged. After a block is finished running, the recorded data is stored in the HDF5 file. Next, one or several actions from the following list are executed:

• the mean and the standard deviation are updated, which are then used to standardize the feedback signal, see section 4.3;

• the interactive Data-driven filter design module is launched to modify or create the DSs based on the data from functional tests;

• a sound signal is issued after the block completion.

FIGURE 5 | Example of the derived (A) and composite (B) signals settings along with the XML script lines the GUI based settings are translated to. Derived Signals description is stored in the vSignals section of the XML file.

• the experiment is suspended;

The following block types are currently supported:

• Baseline. Baseline blocks are useful for recording background activity or running functional tests. In this block type, either a cross is displayed in the user window or a text message with instructions is shown. The data collected during these blocks are utilized for calculation of the standardization parameters and as an input to the Data-driven filter design module for deriving spatial and temporal filters.

• Feedback. This is the main module during which the software displays NFB in the user-selected form. In addition to presenting real NFB, this block supports mock feedback generation. To activate this mode, the user has to specify the name of a prerecorded data file whose data will be used for calculation of the feedback signal. Additionally, to implement mock feedback condition, one can configure the software to use the data file recorded during one of the previous blocks.

• Video. This block type allows a video (specified in the block settings) to be played on the screen.

• Trials. This block type is used in the experiments recording stimulus evoked potentials. On each trial, a stimulus appears on the screen, and the brain response to the stimulus is measured. Different stimuli can be arranged in a predetermined order or randomized. Stimulus and EEG data synchronization is an important issue. One option for syncing the stimuli and the other signals being recorded is to have the stimulus accompanied by a change of brightness within a limited area in the screen corner that is then detected by

a photo-sensor attached to the screen. The sensor signal is then fed into one of the auxiliary channels (see also section 7). Stimuli identifiers and timestamps are saved to the data file. • Center-Out. This block type implements the center-out task (Georgopoulos et al., 1982). The block consists of a sequence of events, such as appearances of screen targets at the central and peripheral positions. The subject moves a cursor (using a joystick or mouse, or through BCI control) from one target to another. The stimuli are synced with the other signals in the data stream. The identifiers and timestamps of all the events are saved to the data file.

Based on these examples, new block types can be added to NFBLab, including the ones that require external visualization programs.

In the example presented in Figure 6, we define an experimental block named Real of type Feedback that will display feedback in the form of a circle with non-harmonic (random) boundary shape. The protocol will not result in an update of the statistics (mean and standard deviation) used for normalization of the feedback signal. The block will last for 120 s.

5.4. Experimental Protocol

The overall experimental protocol is defined by the sequence of blocks entered in the configuration settings (Figure 7). Upon the completion of each block, one or several events are executed, as explained in the previous section.

It is also possible to define groups of experimental blocks that consist of several blocks to be executed in a specific order.

FIGURE 6 | Experimental block settings can be adjusted using the GUI. The parameters specified In the GUI are translated Into the corresponding portion of the XML file. Description of each block is stored in the vProtocols section of the XML file. Here we defined a block named Real of type Feedback with feedback display in the form of a circle with uneven boundary shape. After the protocol is finished the statistics (mean and standard deviation) will not recalculated. Block duration is set to 120 s.

FIGURE 7 | Experimental protocol example. A group of experimental blocks called FBTraining Is defined to consist from of the actual feedback (FB) and rest (Rest) blocks. The repetition count for this unshuffled pair is set to 5. The actual sequence of blocks to be performed is then specified in the vPSequence section of the XML file. As shown in the right panel the experiment itself comprises blocks (Baseline,Open,Close). XML file reflects this construction within vPSequence section.

It is also possible to define randomization within each of the groups. A group of experimental blocks can be set up from the corresponding GUI and will be reflected in the vPGroups section of the XML file. In Figure 7, we define an experimental blocks group called FBTraining that consists of the actual feedback (FB) and rest (Rest) blocks. This sequence of two blocks will be repeated 5 times in this particular order (first FB, then Rest). The actual sequence of blocks to be performed is then specified in the vPSequence section of the XML file. In this example, the experiment comprises blocks (Baseline,Open,Close) and the above defined block group, as illustrated in the right panel of the figure. XML file reflects this construction within vPSequence section.

It is also possible to define a group of blocks with random shuffle requirement. In this setting, block order will

be randomized and the rate of appearance for block will be proportional to the specified value.

All experimental protocol settings can be specified using NFBLab's Experiment protocol editor, implemented as GUI (described next).

6. NFBLAB MODULES

6.1. Experiment Protocol Editor

The experimental protocol can be configured using the graphical user interface (GUI) (see Figure 8). This interface provides access to all the GUI components described in section 6 and allows to fully configure an experiment. The GUI allows to define and modify derived and composite signals, set experimental blocks parameters, organize blocks into groups, and specify the

experimental sequence. These settings are stored in an XML file and can be loaded later if a similar experiment needs to be conducted. After the experimental protocol is defined through the GUI or loaded from a file, the Experimental module GUI (described next) can be launched by pressing the Start button in the bottom of the Experiment protocol editor panel.

6.2. Experimental Module

The experimental module receives raw data from the EEG/ECoG/MEG recording devices, updates all signals, visualizes NFB according to the specifications in each block, controls the sequence of blocks, triggers inter-block events, and saves the experimental data. The experimental module interface (see Figure 9) consists of two main windows: the experimenter's window and the subject's window. In the experimenter's window (Figure 9A) processed signals are displayed in the upper part of the screen and raw signals in the lower part. Textual information about the experiment including the frame rate is displayed at the bottom of the screen. Additionally, the experimenter's window contains start/pause and restart buttons.

During an experiment, the subject faces the monitor that displays the subject's window. Depending on the current block, this window may show different visual stimuli or instructions. The sequence usually starts from a baseline block with a cross displayed in the screen (see Figure 9B). During the actual feedback session (i.e., block of type feedback), a visual indicator of NFB is presented. In the example shown in Figure 9C, the shape of a circle carries NFB information: subjects are instructed to make the circle's boundary smoother by modulating their brain activity.

Raw and processed signals, experimental parameters, signal and block properties are saved as an HDF5 (Hierarchical Data Format) file. This data storage format is widely supported by various libraries, including the ones implemented in Python (h5py package), Matlab (built-in function hdf5read) and R (package h5). The detailed structure of the experiment file can be found in an on-line file5.

6.3. Data-Driven Filters Designer

After the baseline data is collected, NFBLab offers to adjust the signal processing pipeline based on the functional tests. The temporal and spatial filters can be adjusted interactively using the Data-driven filters designer module. The adjustment can start, for example, with an interactive spectral analysis (see Figures 10A,D) that allows to fine-tune EEG frequency bands and redesign the temporal filter to make it, for example, a better match to subject's individual alpha range. Next, spatial decomposition techniques (ICA,SSD or CSP) (see Figures 10A,B) can be invoked to decompose the data into spatial components that represent both functionally relevant neural modulations and artifacts. Using a selector sub-menu (see Figure 10C), it is possible to specify spatial filters that correspond to the features of interest (e.g., sensorimotor rhythm, alpha-rhythm, etc.), and pin-point the spatial components

5 https://github.com/nikolaims/nfb/wiki/Experiment-file-structure

representing the artifacts to be removed. Different methods of spatial decomposition are available within this module.

Independent Component Analysis (ICA) decomposes the multichannel signal into the components representing the functionally relevant modulations and the artifacts (Bell and Sejnowski, 1995). Common Spatial Pattern (CSP) selects the components that contrast signal modulations between the pair of conditions of a functional test. The algorithm is based on solving the generalized eigenvalue problem (Koles et al., 1990) formed by a pair of autocorrelation matrices for the two different conditions being contrasted. Spatio-Spectral Decomposition (SSD) method also solves the generalized eigenvalue problem. This method maximizes the ratio of power in two frequency bands (the central band and the adjacent flanker sub-bands), and allows to find spatial filters that emphasize narrow-band oscillatory activity (Nikulin etal., 2011).

When the filter configuration module is called, the experiment pauses and the target signal manager opens (Figure 10A). At this point, parameters of spatial and temporal filters are shown and can be edited for each DS. As explained above, a spatial filter combines a notch matrix and a filter vector, and converts a multidimensional signal into a scalar. The filter parameters can be edited manually or determined with one of the methods described above (ICA, CSP, or SSD). For each type of analysis, the user can select the components to remove (the notch matrix maps input data into the space orthogonal to the components containing artifacts) and the components to contribute to NFB. The components are selected based on their scalp topography and patterns in the temporal and frequency domains.

7. FEEDBACK PRESENTATION DELAY CONTROL

Rapid presentation of NFB is the distinctive feature of the NFBLab software. To allow for continuous control of visual feedback presentation latency the NFBLab software supports the photo-sensor based loop. We use small square in the upper-right screen corner whose intensity is modulated by the actual feedback signal presented on the main portion of the same screen. The observed within this square intensity fluctuations are registered by the photo-sensor and fed to one of the auxiliary channels of the EEG device to be registered as p(t). When an appropriate photo-sensor has used the shape of p(t) closely matches the corresponding derived signal y(t) used as a feedback signal. Obtained in real-time conditions under causal restrictions signal y(t) is a delayed version of the ideal derived signal yo(t) that can be calculated offline by using zero-phase batch filtering. Thus, signal p(t) appears to be a slightly distorted version of y(t) and is delayed by some specific value with respect to y0(t). This lag can be determined by computing the estimate of cross-correlation

function R(t ) = E{p(t)y°(t and finding the lag Tmax VE{p2(f))E{y;(f))

corresponding to the maximum value Rmax of R(t ). Figure 11A shows the two signalsp(t) andyo(t) and Figure 11B illustrates the cross-correlation function. The lag estimated this way represents the true average total delay between the neuronal event and the moment it is reflected in the feedback presented to the subject.

FIGURE 8 | Graphical user interface of the experiment protocol editor. This interface provides access to all the GUI components described in section 6 and allows to fully configure an experiment.

FIGURE 9 | Graphical user interface of the experiment module: experimenter window (A), subject window for baseline (B), and feedback (C) blocks.

Based on the above considerations it is possible to compute a single pair (Rmax, Tmax) from the data recorded from the entire duration of the experiment in order to assess the average accuracy and latency of the feedback signal. However, given non-stationary nature of the EEG, a more informative approach is to use a sliding window based approach, so that for each data window one can obtain a pair of values characterizing the accuracy and the delay of the feedback that pertains to it and plot the joint distribution of this pair of values over accuracy x latency plane for many possibly overlapping windows. Example of a joint distribution of (accuracy,latency) pairs for the three different envelope estimation methods implemented within NFBLab is shown in Figure 11C. We have chosen the parameters of the three methods to equalize mean envelope estimation accuracy (y-axis). As we can see the Butterworth filter based approach yields the worst accuracy-delay combination as well as the largest variation in the envelope estimation accuracy as compared to the cFIR and the modified windowed FFT techniques. It is also possible to do such computation in a lagged online manner in order to continuously monitor the delay and accuracy parameters of the feedback presented.

8. EXAMPLES OF NFBLAB PROTOCOLS

All participants of experiments bellow provided a written consent approved by The Higher School of Economics Committee on Interuniversity Surveys and Ethical Assessment of Empirical Research in accordance with the Declaration of Helsinki.

8.1. NFB Paradigm Based on the Occipital Alpha Rhythm

Perhaps the simplest NFB implementation is the one that utilizes the prominent occipital alpha rhythm. Commonly, this kind of neurofeedback is performed on the basis of activity in a single EEG channel. Rhythm power is calculated using Fourier transform of 500-1,000 ms data-chunk preceding the feedback presentation time instance. This generally adopted procedure lacks spatial and temporal specificity and could be sub-optimal in terms of training efficiency induced. NFBLab implementation of this simple paradigm addresses both spatial and temporal specificity issues raised above.

To improve spatial specificity we use functional localization approach under which we record and analyze a segment of data collected within the eyes-closed condition and contrast it against the EEG data from the eyes-open condition. Aided with NFBLab interactive tool, we perform CSP spatial decomposition (ICA and SSD are also available) and choose the component with topography corresponding to the occipital generator of the alpha rhythm. This also allows focussing on the alpha frequency band specific to the individual subject. Additionally, the ICA can be performed with an interactive tool (as described in sections 4.1 and 6.3) that chooses spatial components corresponding to the artifact sources to be rejected.

In this example, the band-pass filter is set to the 911 Hz frequency range. The envelope detector comprises complex demodulation followed by smoothing with the 2nd order Savitzky-Golay filter using the 151 samples long window. The resulting delay of envelope estimation (end-to-end)

FIGURE 11 | Estimated p(t) and ideal non-causally estimated envelope yo(t) (A) are separated by a lag that can be assessed using the cross-correlation between p(t) and yo(t) (B). Joint distribution of envelope estimation accuracy (maximum value of the cross-correlation sequence) and total feedback presentation delay (lag of the cross-correlation sequence) for the three different methods (C) as estimated from a large number overlapping windows.

appears to be 131 ms achieved at sampling frequency 250 Hz. The reconstructed envelope matches the ideal Hilbert transform based envelope with an accuracy of 0.7, measured as correlation coefficient. In contrast, FFT and Butterworth filter based approaches, see Figures 3A,B operating at the same accuracy would result into 250-350 ms delay. To further reduce feedback presentation latency and improve temporal specificity of neurofeedback, the novel adaptive complex demodulation approach could also be used (see Figure 3C).

The experiment is divided into 19 blocks (see Figure 12A). During blocks 1 through 6, subjects open and close their eyes ("Open" and "Close" blocks). Within Block 7 ("CSP") of the pipeline, the Data-driven filter designer module is called to derive spatial filter to accommodate individual patterns of the occipital alpha rhythm generator and get rid of artifacts. The CSP analysis contrasts the EEG activity during the "Open" and "Close" blocks. This analysis assumes that the alpha rhythm originates from the same cortical source in both cases, and that the rhythm power increases (i.e., alpha synchronization) in the eyes-closed condition and decreases (i.e., alpha desynchronization) when the eyes are opened. Figures 12C-F shows the characteristics of the selected CSP component, including the segment of its time series containing an alpha episode (C), the spectra of this component for the "Open" and "Close" states (D), the spatial filter (E), and their topography (F).

Following the adjustment of the spatial filter, a resting state (i.e., Baseline block) is recorded to obtain EEG parameters (mean and standard deviation) needed for conversion of the signals into standardized quantities by subtracting the mean and dividing by the standard deviation. The standardization is needed to correctly display the NFB (Figure 12B). For these standardized quantities, the sample mean is 0 and the standard deviation is 1.

The NFB is visualized as a circle with rough (uneven) boundary. The stronger is the NFB signal, the smoother is the boundary. NFB training paradigm consists of five experimental sessions separated by 2-min breaks. Next, five sessions are run with mock NFB. The details of the experimental design are available on-line6. Figure 12H illustrates the processing diagram. The diagram consists of the standard blocks for this kind of experiments and, in addition, specifies unique features for each block, such as settings of spatial and temporal requirements.

Here, as an example, we describe results obtained with this paradigm from a single subject. EEG recordings were carried out using a wireless amplifier SmartBCI (Mitsar Co., Ltd). We collected EEG with 21 channels according to the standard 1020 scheme at 250 Hz sampling rate. The digital average ear (A1 + A2)/2 was used as the reference. As shown in the diagram in Figure 12, the overall duration of the experiment was less than 15 min.

To assess NFB efficiency, we analyzed changes in the NFB power over time within the data segments corresponding to the real and mock NFB training conditions. NFB changes were quantified as the slope of the data-points cloud representing average NFB power within a 1-min time window; 500 windows were randomly allocated to conduct a linear regression analysis. The obtained regression coefficients were significantly different for the mock and real NFB (p < 0.001, Wilcoxon's test). The average slope was significantly positive for real NFB (p < 0.001, Student's t-test, H0-coefficient = 0) and significantly negative for mock NFB (p < 0.001, t-test of the Student, H0 - the coefficient is 0). Thus, only during the real NFB sessions, alpha power experienced statistically significant increase.

6https://github.com/nikolaims/nfb/blob/master/tests/designs/alpha_nfb_settings. xml

FIGURE 12 | Neurofeedback training experiment example: (A) experimental protocol, (B) feedback stimulus, (C) CSP-component time series segment, (D) CSP-component power spectral density (PSD) for closed and open eyes probes, (E) CSP-component spatial filter, (F) CSP-component topography, (G) learning rate in two conditions (Real and Mock feedback), (H) signal processing diagram.

8.2. Motor-Imagery BCI

Brain-computer interfaces are another popular real-time EEG paradigm. The described software contains modules that allow for conducting experiments within P300 and motor-imagery BCI paradigms. Making motor-imagery work reliably in patients with neurological and motor deficits is always a challenge that requires individually designed signal processing tracts. In this example, we consider NFBLab implementation of a lower limb motor-imagery BCI (Figure 13) based on the beta components of the patient's sensorimotor rhythm (SMR) power. This signal is derived from multichannel EEG recordings using a combination of individually designed filters. The first processing stage ("ICA Data") consists of 15 motor-imagery blocks and 15 resting blocks; the duration of each block is 4 s. During a motor-imagery block, the subject imagined limb movements (left hand, right hand, or the legs). Following the collection of this data, "SMR ICA-component extraction" was conducted. The experimenter then used the "Data-driven filter designer" module to perform ICA decomposition of the recorded data and define a spatial filter for extracting SMR. At this stage, it was also possible to adjust the bandwidth that defined SMR and choose the components representing the artifacts to be spatially filtered out from the data. Finally, during the last part of the experiment ("BCI Session"), a BCI loop was implemented for controlling a video game: the output signal was streamed through the LSL interface to a BCI game controlled by two commands: command 1 corresponding

to the resting state and command 2 corresponding to the lower limb motor-imagery state.

We tested this experimental protocol at the Neurotlon competition that took place in St. Petersburg, Russia in 2017 (supported by the interdisciplinary union "Neuronet" and the Ministry of Industry and Trade). BCI game consisted of controlling a virtual submarine that had to avoid obstacles by adjusting its position in the vertical dimension with discrete up or down moves while progressing toward the finish line (Figure 13B). EEG recordings were conducted with NVX52 amplifier (Medical Computer Systems Ltd). We sampled 26 EEG channels at 500 Hz with the reference to A1-A2 (see the schematics of channels in Figure 13E). One paraplegic subject imagined leg movements during the motor-imagery state. The parameters of the spatial filter are shown in Figures 13C-F. Figure 13C shows a representative trace of the ICA component representing SMR. Figure 13D shows the spectrum of this signal, where two distinct peaks are present in the alpha and beta bands for the resting state data (blue). The spectrum shows a clear SMR desynchronization in the motor-imagery state (orange). Subject-specific frequency band (17-21 Hz), marked by vertical lines in Figure 13D, was selected. The spatial-temporal filter assured a good separation between the motor-imagery and resting states (Figure 13E). As depicted in Figure 13G, NFBLab based implementation of the foot motor-imagery BCI allowed the subject to achieve a very low command latency (< 200 ms)

FIGURE 13 | BCI experiment example: (A) experimental protocol, (B) BCI game, (C) ICA-component time series segment, (D) ICA-component power spectral density (PSD) for resting state (Rest) and motor task (Legs), (E) ICA-component spatial filter, (F) ICA-component topography, (G) cross-trials averaged envelope in two states (Rest and Legs).

which resulted in nearly instantaneous gradual control of the game. The details of this experiment design are available online within the NFBLab GitHub resource7 Figure 13H shows the signal processing diagram implemented by the NFBLab in this example.

Moreover, this NFBLab based BCI implementation has been recently used to control the lower-limb exoskeleton in another paraplegic patient with complete leg paralysis. Several challenges had to be addressed in this exoskeleton paradigm, including eliminating mechanical, electrical and EMG artifacts. All these have been successfully handled by the NFBLab software.

9. DISCUSSION

BCI and NFB paradigms are real-time EEG/ECoG/MEG paradigms that are relevant to both basic neuroscience and clinical applications. While the idea of recording neural activity and feeding it back to the subject dates back to the 1960s and 1970s, the BCI and NFB fields have experienced a renaissance during the last decade, with experts from multiple disciplines developing systems that extract and process neural information and convert it into commands to external devices or NFB. Research and development of the modern BCI and NFB systems require advanced software tools. Although the existing software for real-time EEG processing has been

7https://github.com/nikolaims/nfb/master/tests/designs/bci_2states_settings.xml

useful for many applications, we need a next-generation of tools to move forward and address problems related to improving specificity of the feedback signal which will lead to better efficacy of both NFB and BCI technologies. These new software tools should be sufficiently flexible to allow for development of new experimental protocols and their modifications. At the same time, the software should provide a core functionality, with essential EEG processing methods incorporated, such as removal of artifacts, detection and allocation of rhythmic components through spatio-temporal filtering, and evaluation of instantaneous power and phase of oscillatory EEG components. Furthemore, this software should be available to the broad cohort of researchers and should support a wide variety of encephalographic devices, including MEG equipment with unique combination of spatial and temporal resolution properties.

The proposed NFBLab platform offers a solution that meets all these requirements. It supports flexible configuration of the experimental protocol and contains a number of components for EEG data acquisition and processing. These components include functionality that allows to design signal processing chains having high spatial and temporal fidelity. The software inherently supports scripting which in turn allows for natural reproducibility of the paradigms and promotes sharing specific experimental designs by different researchers. Moreover, NFBLab is an open-source project that can incorporate new signal processing procedures and custom modules. The architecture

of the software is modular and can be further extended by the community of users.

NFBLab uses the LSL interface that operates in two directions: multi-channel data acquisition, and streaming results of EEG data processing. This allows, on the one hand, for ensuring compatibility of NFBLab with the equipment of nearly all manufacturers, and on the other hand, implementing a variety of NFB stimuli and controlling an unrestricted variety of external devices.

NFBLab platform employs XML language to define signal processing pipelines, specify characteristics of the experimental blocks and determine the sequence of these blocks to be administered during the experiment. This allows for a fully automatic administration of the experiments with preset or randomized sequences of experimental blocks. Also, mock condition mode can be chosen blindly to the experimenter. This facilitates double-blind designs that must be used to establish the efficacy of NFB interventions. Additionally, XML scripts define the rules for combining different EEG-derived signals and for scaling the feedback based on the signal statistics. We believe that such scripting based approach facilitates reproducible neurofeedback and BCI research. XML scripts can be modified if needed both via NFBLab interface or via a standard text editor.

Although NFBLab offers a significant advantage over existing real-time EEG tools it still lacks certain features. In particular, a graphical programming interface once implemented would allow to create signal processing chains and experimental designs by dragging and dropping icons corresponding to specific signal processing operations and/or experimental protocols. LSL provides universal access to EEG data obtained by a broad variety amplifiers from various vendors. This universality comes at the cost of a delay introduced by the LSL protocol. Although not very significant in most of the applications, it may become a problem in certain cases. Therefore, we believe that implementing direct access to the DAC registers of EEG acquisition devices from some selected vendors. Then, in combination with the advanced signal processing algorithms, we will be able to explore the realm of zero-latency or negative latency feedback for the first time. Currently, only visual feedback is implemented. Adding auditory and tactile feedback options would further reduce efficient feedback delivery time because of the inherent properties of the corresponding sensory systems.

We have been successfully using this software in a wide range of our own projects related to real-time processing of electromyographic, electroencephalographic and electrocorticographic signals. Here we described in detail two implementations of EEG-based NFB and BCI paradigms, and mentioned the third one. The constantly growing functionality

of the software goes beyond these examples. The compatibility of NFBLab platform with external hardware and software and its flexibility in adjusting experimental parameters make it useful for numerous research fields and clinical applications. In conclusion, NFBLab is a versatile software package that eases development of NFB systems and employs advanced algorithms to improve NFB efficiency. As such, it advances a variety of NFB and BCI applications, including NFB-based treatment of neurological conditions.

10. SOFTWARE AVAILABILITY

The source code and installation instructions for NFBLab can be found at https://github.com/nikolaims/nfb. The manual is available at https://nfb-lab.readthedocs.io/en/latest/index.html. In order to facilitate learning how to use our software, we have developed a tutorial and included it in the documentation. In the provided example, the EEG device is emulated by an LSL stream created within the NFBLab and streaming the prerecorded data read from a file. After examining this example, a prospective user can follow the tutorial to get familiar with the basic concepts of NFBLab. The XML files implementing the NFB and BCI experiments described in section 8 are available at https://github.com/nikolaims/nfb/blob/master/tests/designs/ alphanfbsettings.xml and https://github.com/nikolaims/nfb/ master/tests/designs/bci2statessettings.xml.

AUTHOR CONTRIBUTIONS

NS developed the software and wrote the first version of the manuscript. KV tested the software and added application examples. SZ developed low-level routines for equipment support and maintained the equipment within application examples. ML provided general guidance, wrote and revised the paper. AO conceived the study, designed the software, developed the earlier version of the software, wrote and revised the paper.

FUNDING

The study has been funded by the Center for Bioelectric Interfaces NRU Higher School of Economics, RF Government grant, ag. No 14.641.31.0003.

ACKNOWLEDGMENTS

The authors would like to thank Dr. Mark E. Pflieger for the useful discussions related to the architecture of the initial version of the presented software.

REFERENCES

Ang, K. K., Chua, K. S., Phua, K., Wang, C., Chin, Z. Y., Kuah, C. W. et al. (2015). Randomized controlled trial of EEG-based motor imagery brain-computer in-terface robotic rehabilitation for stroke. Clin. EEG Neurosci. 46, 310-320. doi: 10.1177/1550059414522229

Bell, A., and Sejnowski, T. (1995). An information-maximization approach to blind separation and blind deconvolution. Neural Comput. 35, 83-105.

Coben, R., Linden, M., and Myers, T. (2010). Neurofeedback for autistic spectrum disorder: a review of the literature. Appl Psychophysiol. Biofeedback. 7, 11291159. doi: 10.1007/s10484-009-9117-y

Congedo, M., Lubar, J. F., and Joffe, D. (2004). Low-resolution electromagnetic tomography neurofeedback. IEEE Trans. Neural Syst. Rehabil. Eng. 12,387-397. doi: 10.1109/TNSRE.2004.840492 Darvas, F., Pantazis, D., Kucukaltun-Yildirim, E., and Leahy, R. (2004). Mapping human brain function with MEG and EEG: methods and validation. Neuroimage. 23(Suppl. 1):S289-S299. doi: 10.1016/j.neuroimage.2004.07.014 Enriquez-Geppert, S., Huster, R. J., and Herrmann, C. S. (2017). EEG-neurofeedback as a tool to modulate cognition and behavior: a review tutorial. Front. Hum. Neurosci. 11:51. doi: 10.3389/fnhum.2017.00051 Frolov, A., Kozlovskaya, I., Biryukova, E., and Bobrov, P. (2018). Use of robotic devices in post-stroke rehabilitation. Neurosci. Behav. Physiol. 48, 1053-1066. doi: 10.1007/s11055-018-0668-3 Georgopoulos, A. P., Kalaska, J. F., Caminiti, R., and Massey, J. T. (1982). On the relations between the direction of two-dimensional arm movements and cell discharge in primate motor cortex. J. Neurosci. 2, 1527-1537. doi: 10.1523/JNEUROSCI.02-11-01527.1982 Greenblatt, R., Ossadtchi, A., and Pflieger, M. (2005). Local linear estimators for the bioelectromagnetic inverse problem. IEEE Trans. Signal Proc. 53, 3403-3412. doi: 10.1109/TSP.2005.853201 Hamalainen, M., and Ilmoniemi, R. (1994). Interpreting magnetic fields of the

brain: minimum norm estimates. Biol. Eng. 32, 35-42. Kamiya, J. (1969). "Operant control of the EEG alpha rhythm and some of its reported effects on consciousness," in Altered States of Consciousness, ed C. Tart (New York, NY: John Wiley and Sons), 489-501. Koles, Z. J., Lazaret, M. S., and Zhou, S. Z. (1990). Spatial patterns underlying population differences in the background EEG. Brain Topogr. 2, 275-284. doi: 10.1007/BF01129656 Larsen, S., and Sherlin, L. (2013). Neurofeedback: an emerging technology for treating central nervous system dysregulation. Psychiatr. Clin. North Am. 36, 163-168. doi: 10.1016/j.psc.2013.01.005 Lee, W. R., Caccetta, L., Teo, K. L., and Rehbock, V. (2006). Optimal design of complex fir filters with arbitrary magnitude and group delay responses. IEEE Trans. Signal Process. 54, 1617-1628. doi: 10.1109/TSP.2006. 872542

Lofthouse, N., Hendren, R., Hurt, E., Arnold, L. E., and Butter, E. (2012). A review of complementary and alternative treatments for autism spectrum disorders. Autism Res. Treat. 2012:870391. doi: 10.1155/2012/870391 Nikulin, V. V., Nolte, G., and Curio, G. (2011). A novel method for reliable and fast extraction of neuronal EEG/MEG oscillations on the basis of spatio-spectral decomposition. Neuroimage 55, 1528-1535. doi: 10.1016/j.neuroimage.2011.01.057 Oblak, E. F., Lewis-Peacock, J. A., and Sulzer, J. S. (2017). Self-regulation strategy, feedback timing and hemodynamic properties modulate learning in a simulated fmri neurofeedback environment. PLoS Comput. Biol. 13:e1005681. doi: 10.1371/journal.pcbi.1005681 Ossadtchi, A., Shamaeva, T., Okorokova, E., Moiseeva, V., and Lebedev, M. (2017). Neurofeedback learning modifies the incidence rate of alpha spindles, but not their duration and amplitude. Sci. Rep. 7:3772. doi: 10.1038/s41598-017-04012-0

Pascual-Marqui, R. (2002). Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. Methods Find Exp. Clin. Pharmacol. 24(Suppl. D):5-12.

Ponomarev, V., Mueller, A., Candrian, G., Grin-Yatsenko, V., and Kropotov, J. (2012). Group independent component analysis (gICA) and current source density (CSD) in the study of EEG in ADHD adults. Clin. Neurophysiol. 125, 83-97. doi: 10.1016/j.clinph.2013.06.015 Renard, Y., Lotte, F., Gibert, G., Congedo, M., Maby, E., Delannoy, V., et al. (2010). Openvibe: an open-source software platform to design, test and use brain-computer interfaces in real and virtual environments. Presence 19, 35-53. doi: 10.1162/pres.19.1.35 Schalk, G., McFarland, D., Hinterberger, T., Birbaumer, N., and Wolpaw, J. (2004). BCI2000: a general-purpose brain-computer interface (BCI) system. IEEE Trans. Biomed. Eng. 51, 1034-1043. doi: 10.1109/TBME.2004.827072 Shin, H., Law, R., Tsutsui, S., Moore, C., and Jones, S. (2017). The rate of transient beta frequency events predicts behavior across tasks and species. eLife 6:e29086. doi: 10.7554/eLife.29086 Sitaram, R., Ros, T., Stoeckel, L., Haller, S., Scharnowski, F., Lewis-Peacock, J., et al. (2016). Closed-loop brain training: the science of neurofeedback. Nat. Rev. Neurosci. 18, 86-100. doi: 10.1038/nrn.2016.164 Smetanin, N., and Ossadtchi, A. (2017). "Express estimation of brain rhythm power for low-latency neurofeedback," in The First Biannual Neuroadaptive Technology Conference (Berlin). Sterman, M. B., MacDonald, L. R., and Stone, R. K. (1974). Biofeedback training of the sensorimotor electroencephalogram rhythm in man: effects on epilepsy. Epilepsia 15, 395-416. doi: 10.1111/j.1528-1157.1974.tb04016.x Thatcher, R., and Lubar, J. (2008). Z Score Neurofeedback. Clinical Applications. San

Diego, CA: Academic Press. White, D. J., Congedo, M., and Ciorciari, J. (2014). Source-based neurofeedback methods using EEG recordings: training altered brain activity in a functional brain source derived from blind source separation. Front. Behav. Neurosci. 8:373. doi: 10.3389/fnbeh.2014.00373 Wolpaw, J. (2012). Brain-Computer Interfaces: Principles and Practice. Oxford:

Oxford University Press. Zoefel, B., Huster, R. J., and Herrmann, C. S. (2011). Neurofeedback training of the upper alpha frequency band in EEG improves cognitive performance. Neuroimage 54, 1427-1431. doi: 10.1016/j.neuroimage.2010.08.078

Conflict of Interest Statement: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Smetanin, Volkova, Zabodaev, Lebedev and Ossadtchi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

B Статья "Digital filters for low-latency quantification of brain rhythms in real time"

Authors. Nikolai Smetanin, Anastasia Belinskaya, Mikhail Lebedev and Alexei Ossadtchi

Abstract. Objective. The rapidly developing paradigm of closed-loop neuroscience has extensively employed brain rhythms as the signal forming real-time neurofeedback, triggering brain stimulation, or governing stimulus selection. However, the efficacy of brain rhythm contingent paradigms suffers from significant delays related to the process of extraction of oscillatory parameters from broad-band neural signals with conventional methods. To this end, real-time algorithms are needed that would shorten the delay while maintaining an acceptable speed-accuracy trade-off. Approach. Here we evaluated a family of techniques based on the application of the least-squares complex-valued filter (LSCF) design to real-time quantification of brain rhythms. These techniques allow for explicit optimization of the speed-accuracy trade-off when quantifying oscillatory patterns. We used EEG data collected from 10 human participants to systematically compare LSCF approach to the other commonly used algorithms. Each method being evaluated was optimized by scanning through the grid of its hyperparameters using independent data samples. Main results. When applied to the task of estimating oscillatory envelope and phase, the LSCF techniques outperformed in speed and accuracy both conventional Fourier transform and rectification based methods as well as more advanced techniques such as those that exploit autoregressive extrapolation of narrow-band filtered signals. When operating at zero latency, the weighted LSCF approach yielded 75% accuracy when detecting alpha-activity episodes, as defined by the amplitude crossing of the 95th-percentile threshold. Significance. The LSCF approaches are easily applicable to low-delay quantification of brain rhythms. As such, these methods are useful in a variety of neurofeedback, brain-computer-interface and other experimental paradigms that require rapid monitoring of brain rhythms.

Journal of Neural Engineering

<D

Cross Mark

RECEIVED

29 November 2019

REVISED

23 March 2020

ACCEPTED FOR PUBLICATION

14 April 2020

PUBLISHED

29 July 2020

PAPER

Digital filters for low-latency quantification of brain rhythms in real time

Nikolai Smetanin®, Anastasia Belinskaya, Mikhail Lebedev and Alexei Ossadtchi©

Center for Bioelectric Interfaces, Higher School of Economics, Moscow, 101000, Russia

E-mail: [e-mail], [e-mail], [e-mail] and [e-mail]

Keywords: brain-state contingent paradigm, closed-loop neuroscience, neurofeedback, brain rhythm, envelope, instantaneous phase, feedback latency, state-dependent stimulation, optimal filtering, Hilbert transform, complex-valued filter, prediction

Abstract

Objective. The rapidly developing paradigm of closed-loop neuroscience has extensively employed brain rhythms as the signal forming real-time neurofeedback, triggering brain stimulation, or governing stimulus selection. However, the efficacy of brain rhythm contingent paradigms suffers from significant delays related to the process of extraction of oscillatory parameters from broad-band neural signals with conventional methods. To this end, real-time algorithms are needed that would shorten the delay while maintaining an acceptable speed-accuracy trade-off. Approach. Here we evaluated a family of techniques based on the application of the least-squares complex-valued filter (LSCF) design to real-time quantification of brain rhythms. These techniques allow for explicit optimization of the speed-accuracy trade-off when quantifying oscillatory patterns. We used EEG data collected from 10 human participants to systematically compare LSCF approach to the other commonly used algorithms. Each method being evaluated was optimized by scanning through the grid of its hyperparameters using independent data samples. Main results. When applied to the task of estimating oscillatory envelope and phase, the LSCF techniques outperformed in speed and accuracy both conventional Fourier transform and rectification based methods as well as more advanced techniques such as those that exploit autoregressive extrapolation of narrow-band filtered signals. When operating at zero latency, the weighted LSCF approach yielded 75% accuracy when detecting alpha-activity episodes, as defined by the amplitude crossing of the 95th-percentile threshold. Significance. The LSCF approaches are easily applicable to low-delay quantification of brain rhythms. As such, these methods are useful in a variety of neurofeedback, brain-computer-interface and other experimental paradigms that require rapid monitoring of brain rhythms.

1. Introduction

Investigations of neural oscillations have been and continue to be an area of intensive research, particularly with the advancement of neuroimaging techniques, such as noninvasive electroencephalography (EEG) and magnetoencephalography (MEG), invasive electrocorticography (ECoG) and stereo EEG (sEEG). Many experimental paradigms and methods have been developed in order to deal with specific types of neuronal oscillations particularly, experimental and clinical approaches for either enhancing or suppressing various brain rhythms [1, 2]. Existing paradigms for exploration of oscillatory activity

fall in one of the two categories: (1) studies where stimuli are presented regardless of the specifics of the ongoing neural activity and (2) studies implementing a closed-loop design where the stimuli formation depends on the characteristics of the ongoing brain activity [3]. In this paper, we evaluate signal processing methods that are useful for the second category where the closed-loop is formed on the basis of the oscillatory activity quantified in real-time. We refer to this experimental approach as brain-rhythm contingent paradigm (BRCP) .

As shown in figure 1(A), BRCP operates via three key steps: (1) data acquisition, (2) target signal extraction, and (3) stimulus generation. During the data

Figure 1. Schematics of BRCP. (A) Signal flow diagram in the closed-loop BRCP system (B) Signal processing pipeline (C) The sources of delays contributing to the total latency of BRCP. Delays with technical and fundamental sources are marked in blue and red, respectively.

acquisition step, brain activity is measured, usually with multiple spatially distributed electromagnetic sensors, and streamed to a computer. Within the target signal extraction step, a computer routine extracts the parameters of oscillatory activity from the multichannel data in real-time, typically amplitude and phase. At the stimulus generation step, these parameters are converted into a feedback delivered to the brain either directly, using stimulation applied to the nervous tissue (also called neuromodulation), or using natural senses: vision, hearing or touch. The feedback can be presented continuously and modulates the parameters of an ongoing stimulation, or discretely and contributes to selecting a stimulus from a predefined set.

Many implementations of BRCP have been developed allowing us to explore a variety of goal-directed behaviors dependent on a closed-loop design [3]. The most distinct versions of BRCP are: neurofeedback [4-7], brain computer interface (BCI) [8-10], closed-loop brain stimulation [11-13], and brain state-contingent stimulus delivery [14-20]. As such, BRCPs are applicable to the treatment of a range of neurological disorders [21], [22], cognitive enhancement therapy [23], and post-stroke rehabilitation [9], [10], [24]. Additionally, brain-state contingent stimulus presentation can reduce the overall session duration and reduce participants' fatigue, particularly when neural patterns of interest consist of short-lived episodes of activity [25, 26]. Despite the conceptual differences between the

implementations of BRCP reported in the literature, their signal processing pipelines have many common features, as depicted in figure 1(B). The generic pipeline incorporates the spatial and temporal filtering steps in order to emphasize specific neural sources and frequency components, respectively. These filtering steps are often followed by quantification of the oscillatory envelope and phase.

Here we focus primarily on the multichannel recording methods with millisecond-scale temporal resolution, such as EEG and MEG. These techniques allow for exploration of the fine rhythmic structure of brain activity and, in principle, can enable closed-loop, instantaneous interaction with brain circuits. Yet, as explained in detail below, time-lags of both technical and fundamental nature occur during the online extraction of oscillatory parameters from the ongoing brain activity. These delays vary, depending on the concrete implementation and can be significantly reduced with the use of optimized signal processing methods.

Temporal specificity of BRCP is characterized by the overall delay between the onset of the neural event of interest and the time the corresponding feedback is issued (figure 1(C)). This overall BRCP latency incorporates time-lags related to different factors, some of them technical (i.e. software and hardware delays), others, fundamental (i.e. required to collect a sufficient amount of neural data to isolate a rhythmic component). The technical time-lags typically do not exceed 100 ms; this time is needed for data

transmission by the hardware and low-level software processing. These delays can be reduced by specific hardware and software solutions. From our experience, it is feasible to reduce the technical delay to 10-30 ms depending on the feedback delivery specifics. The fundamental lag cannot be reduced as easily because it is composed of the time needed to collect a snapshot of neural oscillatory data that could be then quantified with an appropriate algorithm like bandpass filtering and envelope/phase analysis. As such, this lag is fundamental: a longer observation time is required to improve the resolution over spectral frequencies. Suboptimal approaches used for extracting the envelope and/or phase of rhythmic neural components may mount to several hundreds of milliseconds of fundamental delay. This is undesirable in most closed-loop paradigms [27].

The adverse effect of feedback delays was reported in many neurofeedback applications and the BCI studies where participants were supposed to learn to modulate their own brain activity - with or without a mental strategy recommended by the experimenter [28]. A similar problem with long feedback delays was highlighted in several studies on evidence-based learning. Thus, back in 1948, Grice showed that learning to discriminate complex visual patterns drastically depends on the feedback signal latency [29]. Impaired performance caused by delayed feedback was also demonstrated in the studies on motor learning, e.g in the prism adaptation task [30]. According to Rah-mandad etal [20], behavioral learning is impaired when feedback delay is misperceived by a participant. Moreover, an elongated feedback delay impairs the sense of agency during the BCI control [31], the finding also corroborated by the simulation study [32] showing that feedback delay and temporal blur adversely influence automatic (strategy free) learning in BCI tasks.

Temporal specificity is also an important consideration for the experimental settings with closed-loop brain stimulation and brain state triggered stimulus delivery [13, 33, 34]. Indeed, a key requirement for these methods is an accurate estimation of instantaneous oscillatory features and the timely delivery of stimuli to efficiently interfere with oscillatory neural patterns.

In the present study, we evaluate several approaches aimed at reducing the delay between neuronal events and the corresponding feedback in BRCP. We pay special attention to a family of methods based on the least-squares complex-valued filter (LSCF) design methodology to identify finite impulse response (FIR) filter weights and explicitly control the processing latency. The proposed approach is capable of extracting the instantaneous amplitude and phase of neural oscillations with a shorter latency and higher accuracy as compared to other

techniques. With LSCF, the user can explicitly specify the desired delay which facilitates a flexible control of the speed-accuracy trade-off in the closed-loop neuroscientific experiments..

2. Methods

Our basic assumption is that the measured neural activity x[n] is a sum of the narrow-band signal s[n] (neural activity targeted by a BRCP) and colored broad-band noise n[n].

x[n] = s[n] + n[n] (1)

The narrow band neural activity s[n] utilized in a BRC paradigm can be represented as the real part of analytic signaly[n] [35]:

y[n] = s[n] + jsh [n] = a[njn] (2)

where Sh[n] is the imaginary part of the analytic signal (often called second quadrature of the original signal), a[n] is the instantaneous amplitude of the narrowband process at time n, and ^>[n] is its instantaneous phase at time n. Sequence a[n], often called envelope, is a relatively smooth real-valued sequence that skirts the local peaks of a narrow-band process. Importantly, envelope a[n] and phase ^>[n] can be instantaneously computed from the complex-valued analytic signal y[n] as

a[n] = V^{y[n]}2 + 3{y[n]}2,

*] = iH) (3)

As illustrated in figure 2(A), analytic signal y[n] can be computed from x[n] using narrow-band Hilbert transform [35] through the following steps. First, the discrete Fourier transform (DFT) is applied to the data segment of x[n]. Then, the obtained DFT coefficients are masked so that the coefficients corresponding to the frequency band of interest (designated on the positive-frequency semi-axis) are doubled and the others are set to zero. Lastly, the inverse Fourier transform (iDFT) is applied to obtain y[n]. This approach will inevitably deliver inaccurate values of y(n) for the beginning and the end of the data segment it is applied to. To avoid this distortion, the segment needs to be augmented with flanker data—a procedure that is impossible to conduct in real time where the data segment ends at the last sampled point. This is because obtaining the flanker in this case would require an extension into the unknown future. Unlike the real-time operations, y[n] is estimated accurately when the data segment is extended with true flanker data in the offline analysis. Then, using equations (3)

Figure 2. Methods for narrow-band signal envelope estimation. (A) Ideal non-causal system for ground truth signal extraction. (B) Envelope detector based on rectification of the band-filtered signal. (C) Sliding window narrow-band Hilbert transform-based method. (D) Complex-valued FIR family filters that perform FIR causal filter based approximation of the ideal non-causal system. (D') Complex-valued FIR filter design.

one can compute the ideal envelope a[n] and phase 0[n] sequences. The rest of figure 2 shows a summary of the available causal techniques that can be applied in real-time to compute the envelope and instantaneous phase of a narrow-band component extracted from a broad-band signal.

2.1. Conventional methods

For didactic purposes, we start with the most basic approaches to extraction and quantification of the narrow-band processes.

2.1.1. Rectification and smoothing of the band-filtered signal (rect)

This is conceptually the most straightforward way to estimate the envelope a[n]. With this method, low-pass filtering is applied to the rectified narrowband filtered signal, as illustrated in figure 2(B). This approach can be mathematically expressed as:

where * is the convolution operator, hbp is the impulse response of the band-pass filter, | ■ | denotes the absolute value (i.e. the rectification step), and hp is the impulse response of the low-pass filter that performs smoothing of the rectified signal. Without loss of generality, we can assume that both hbp and hip are linear phase FIRs designed by a Hanning window method with the number of taps Nbp and Nip, correspondingly. The cutoff frequency for the low-pass filter fc is set to correspond to one-half of the expected bandwidth of the narrow-band process, i.e. fc = (f — f1 )/2. FIR filters have a linear phase and therefore the total delay D and the number of taps in the individual symmetric filters Nbp and Nip are related as:

D = NP - 1 + NbP - 1

(5)

ar[n] = hip[n] * hbp[n] *x[n]

(4)

It follows from equation (5) that a given value of D can be achieved for multiple selections of Nip and Nbp. In order to ensure the best envelope reconstruction accuracy for a given group delay value D, we

used grid search over variable Nbp in our comparative analysis. The parameter Nip then follows directly from (5). Since Nbp and Nip are positive, this method can estimate the envelope values only for a positive delay.

2.1.2. Sliding window narrow-band Hilbert transform (hilb)

This is a second most commonly used method that is based on the use of the analytic signal y[n] computed using Hilbert transform [35]. There are various implementations of this approach. In this paper, we resorted to the use of the windowed discrete Fourier transform (DFT). The DFT is calculated on each window of length Nt which is zero-padded to the length of Nf samples. We then proceed as in the ideal approach but after computing the inverse DFT we retain only the (Nt —D)-th element of the resultant complex-valued sequence as an estimate of the analytic signal with delay D. This filter performs two operations simultaneously: band-pass filtering and extracting the analytic signal. The analytic signal is then used to instantaneously estimate the oscillatory envelope and phase. This algorithm is illustrated in figure 2(C).

In matrix representation and using temporal embedding to form vector x[n], this can be represented as:

yh [n] = Nf ' wH<-DWAfx[n] >

(6)

f

Ne = Na corresponding to the central point of the extrapolated segment of length 2Ne. In the original work, this value was used to determine the current phase that was then extrapolated for the next cycle to determine the time for stimulation. This approach was originally designed for D = 0, so we used it as a benchmark when we comparing the performance of different methods operating at zero latency.

Using matrix notation, we can formalize this method as follows.

2

yp [n] = Nff ' wNa W+PAR(p){ x[n] }

(7)

where vector x[n] contains the last Nt samples of x[n], i.e. x[n] = [x[n],x[n — 1],...x[n — Nt + 1]] , WAf is the Nf -by-Nt modified DFT matrix with zeros in the k-th row for k outside of the range [f1Nf/fs, f2Nf/fs], which corresponds to the band of interest Af= [f1, f2] and wH(_

_D is the Hermitian transposed (Nt —D)-th column of the DFT matrix.

The overall delay of this method is explicitly determined by parameter D. The parameters to be optimized in this method are window length Nt and zero-padded length of the signal Nf.

2.1.3. Sliding window Hilbert transform with autoregressive prediction of the narrow-band filtered signal (ffiltar)

This approach was proposed by Chen et al [33] and applied in brain rhythm contingent transcranial magnetic stimulation (TMS) [13]. In this method, the sliding window vector x[n] containing Na last samples is forward-backward bandpass filtered. Then, flanker Ne samples are truncated to eliminate edge artifacts, and 2Ne samples are forward predicted by using an autoregressive (AR) model fitted to the immediately preceding Na — Ne samples. Finally, Hilbert transformation of the prediction is used to estimate the analytic signal with zero latency at n = Na — Ne +

where x[n] contains forward-backward filtered last Na samples of the x[n], PAR(p) denotes AR model based prediction operation that augments the truncated data segment and adds 2Ne predicted samples by using a p-th order AR model, W+ is the Nf -by-(Na + Ne) modified DFT matrix with zeroed rows corresponding to the range of negative digital frequency values [—n, 0] and w$ is the Hermitian transposed Na-th row of the DFT matrix. We used a Butterworth filter of the order k as the narrow-band filter for the forward-backward filtering part as suggested in [33].

One of the disadvantages of this approach is that it has multiple parameters that need to be adjusted to achieve the optimal performance. To attain the best performance in this study, we searched over the parameter grid composed of the following variables: AR order p, number of flanker samples Ne, and Butterworth filter's order k.

2.2. Least squares complex-valued filter (LSCF) design approach

In order to build the analytic signal y[n] that corresponds to an estimate of narrow-band signal s[n] extracted from the noisy broad-band measurements x[n], one can apply the ideal narrow-band Hilbert transform filter [36, 37]. The complex-valued frequency response of this combined filter can be defined as:

HD(ew)=

2e-jwD, w e [wc — Sw, wc + Sw] C [0,n] 0, otherwise

(8)

where Sw = 2nSf is half of the pass band width, and D is the group delay measured in samples. This filter is non-causal for any finite delay D and therefore it cannot be applied in real time. To reconstruct the analytic signal causally, one can find a causal complex-valued finite impulse response (FIR) filter b = {b[n]} of length Nt that approximates the ideal complex-valued frequency response HD(ejw) [38]. This filter can be then applied in real-time asyc[n] = x[n]*b[n] using the convolution of b and x[n], the procedure

that by design incurs a fixed processing delay of D samples.

Causal complex-valued FIR filter b[n] can be established by solving the least squares optimization problem. Various definitions of the cost functions lead to different filters.

2.2.1. Frequency domain least squares (cfir) This is the most straightforward approach to quantifying narrow-band components present in the data. The least squares filter design strategy consists of finding the complex-valued vector of the FIR filter weights b of length Nt that minimizes the L2 norm between the FIR filter frequency response obtained by the DFT and the discrete appropriately sampled version hD of an ideal response HD(eju) in the frequency domain. To increase the frequency resolution, we use a truncated DFT-matrix W with the dimensions Nf x Nt (Nf > Nt). The use of transform matrix-based formulation of the DFT in this case is equivalent to the DFT of Nt samples long vector zero-padded to length Nf. Since WHW = NfI, the solution bcfir of the normal equations for the optimization problem (2.2) can be found by a simple inverse DFT of the desired complex-valued characteristics of the narrow-band Hilbert transformer:

}cfir :

argmin ||Wb — hD||L2 b

1

bcfir = N WHhD

(9)

where Sx is the diagonal matrix formed from the samples of the power spectral density of x[n] estimated over the entire training data segment. Note that only WHW = NfI remains true while WWH = NtI. Therefore, at the optimum ||Wb — hD|| = 0 and thus b is the least squares approximation of the ideal filter that can be computed even for negative D.

Panels D and D' of figure 2 illustrate cfir and wcfir approaches. The delay D corresponds to the slope of the phase response within the pass-band and theoretically can be set to an arbitrary value. The analytic optimization procedure then aims at finding such complex-valued vector of the FIR filter coefficients b for which the desired transfer function hD is approximated optimally in the least squares or weighted least squares sense.

Conceptually, having in mind the two tasks of optimal envelope and instantaneous phase estimation we could have formulated two separate optimization problems and used two different sets of weights implementing two different band-pass complex-valued filters delivering optimal accuracy in estimation of envelope and phase approximation with the specified delay. In this case, however, due to non-linearity of the target functional, we would have to perform an iterative optimization in order to find the optimal vector of weights b for the FIR filter.

2.2.3. Time domain least squares (tcfir) This is the last approach from cfir family. It is based on minimization of the squared distance in the time domain between the complex delayed ground truth signaly[n — D] and the filtered signal x[n]*b[n]:

Equation (9) is similar to equation (6) but can be used with negative delays, D < 0. This simple method, however, does not take into account the second order frequency domain statistics of the target signal and could be further improved. Note that the cfir approach with parameters Nt, Nf and D > 0 matches the sliding window narrow-band Hilbert transform approach with the same parameters Nt, Nf and D.

2.2.2. Frequency domain weighted least squares (wcfir) This is the method that optimizes the filter by constructing weights from the power spectral density of the input signal x[n]. We thus formulate the weighted frequency domain least squares design technique via optimization problem

bwcfir = argmin ||S1/2(Wb — hD )||l2 (10) b

whose solution can be found by solving the normal equations:

btcfir = argmin ||x[n] * b[n] — y[n — D]||l2 b

(12)

^wcfir

(WHSxW) WHSxh

XhD

(11)

The ground truth signal y[n — D] is obtained non-causally from the training samples via an ideal zero-phase Hilbert transformer (8). According to Par-seval's theorem, this approach is equivalent to the wcfir approach. However in contrast to the frequency domain formulation, it allows for implementation of recursive schemes (12) and, therefore, may potentially account for non-stationarity in the data. One of the most straightforward solutions is to use recursive least squares [39] in order to update filter weights online.

2.3. Methods comparison

2.3.1. Data and preprocessing

We compared the performance of the methods described above using resting state EEG data recorded from ten participants engaged in neurofeedback experiments. Resting state data was collected prior to neurofeedback runs. Prior to the experiment, all participants provided written consent approved by The

Higher School of Economics Committee on Inter-University Surveys and Ethical Assessment of Empirical Research in accordance with the Declaration of Helsinki. EEG data was obtained using 32 AgCl electrodes placed according to the 10-20-system with the ground electrode at AFz position and reference electrodes on both ears. The impedance for all electrodes was kept below 10 KOhm. The signal was sampled at 500 Hz using an NVX-136 amplifier (Medical Computer Systems Ltd) and bandpass-filtered in the 0.5—70 Hz band. These preprocessing filters incurred an overall delay of no more than 10 milliseconds for the bandwidth of interest (8-12 Hz).

We used 2 minutes of resting state recordings for each participant. The first minute of the data was used for implementing a parametric grid search and the second minute - for performance testing. These processing steps corresponded to the settings of our neurofeedback experiments where, after recording a subject's resting state activity for one minute, we adjusted the signal processing parameters and then proceeded to the main body of the experiment. To eliminate eye artifacts, we performed independent component analysis (ICA) on the training data, identified eye movement-related components by means of the mutual information spectrum [40], and removed the artifactual data characterized by the highest mutual information with the two frontal channels Fp1 and Fp2. Only the parietal P4 channel of the cleaned data was used for the analysis and it played the role ofx[n] in (1).

2.3.2. Individual alpha band

We determined individual alpha ranges with the following procedure: estimate magnitude spectrum using Welch method with a 4 second 90%-overlap boxcar window, find the frequency f0 with maximal signal-to-noise ratio (SNR) in the 8-12 Hz range, define individual band as [ f1,f2] interval, where f1 = f0 — 2 Hz and f2 = f0 + 2 Hz. SNR was defined as the difference between mean magnitude in the band Aband = [ f1,f2] and mean magnitude in the flankers band Aflankers = [f — 2,f1]U[f2,f2 + 2] divided by the mean magnitude in the flankers band Aflankers:

SNR -

m(Aband) — m(Aflankers)

m(Aflankers)

where m(A) is the mean magnitude in band A.

(13)

2.3.4. Performance metrics

To measure the accuracy of the envelope and instantaneous phase estimation obtained with the described techniques, we used the following metrics. All metrics were based on the ground truth envelope and instantaneous phase information extracted from the ground truth signal s[n] obtained non-causally from the real EEG data as described above and then shifted to match the specific delay value D of causal processing. We will denote the shifted ground truth envelope and phase sequences as a[n — D] and ¿[n — D] correspondingly.

To assess the performance of different methods for estimating the envelope, we calculated the Pearson correlation coefficient between the estimated envelope a[n], obtained causally with each of the methods, and the appropriately shifted ground truth envelope sequence a[n — D] over the 1-minute long test data segment:

J2 (a[n — D] — ma)(a[n] — ma)

neNa

(a[n — D] — ma)2

neNa

J2 (a[n] — ma)2

neNa

(14)

where Na = D... N — 1 is the set of time indices with N = 30000, ma and m-a are sample mean values of a[n] and a[n] over set Na

To assess the performance of instantaneous phase estimation, we used bias b^, absolute bias |b^| and the standard deviation a^ with respect to the delayed ground truth phase ¿[n — D] at the time when predicted 4>[n\ phase crosses 0. These metrics reflect the bias, absolute bias, and the variance in determining zero crossings(negative-to-positive direction) of the delayed signal s[n — D].

b* =

IN I

E ¿[n — D]

neNfi

(15)

a^ =

E —D]—^>2 (16)

2.3.3. Ground truth signal

As the ground truth signal we used the non-causally computed analytic signal y[n] obtained using the ideal procedure described in section 2 and illustrated in figure 2(A) and corresponding to the individual alpha band [f1,f2] of each subject. After the analytic signal y[n] was computed, we estimated both envelope and instantaneous phase without any additional delay according to equation (3).

whereA^ = {n : n e Na, sign(<p[n]) > sign(4>[n — 1])} is the set of time instances when 4>[n] crosses 0.

2.3.5. Grid search procedure

To ensure that the methods being compared operated optimally, we defined a grid search space for each of them as described in table 1 and looked for the combination ofparameters that ensured the best performance for each technique. Here we use the following short method names: rect for envelope detector based on rectification of the band-filtered signal, cfir for the frequency domain LSCF, wcfir for the frequency

r

a

1