Извлечение компонентов ЭЭГ из сигнала в MATLAB

У меня есть простой сигнал ЭЭГ в MATLAB, такой как показан на рисунке ниже. И я хотел извлечь компоненты ЭЭГ в соответствии со следующей таблицей.

  • Дельта - до 4 Гц;
  • Тета - 4 -> 8 Гц
  • Альфа - 8 -> 13 Гц
  • Бета - 13 -> 30 Гц
  • Гамма - 30 -> 100 Гц

В первой попытке решить эту проблему я попытался создать полосовой фильтр с помощью 'fdatool' из MATLAB для извлечения сигнала 'theta' компонента, но безуспешно.

Вдобавок приложил код для фильтра, полученный с помощью 'fdatool'.

function Hd = filt_teta
%FILTROPARA TETA Returns a discrete-time filter object.

%
% M-File generated by MATLAB(R) 7.9 and the Signal Processing Toolbox 6.12.
%
% Generated on: 05-May-2011 16:41:40
%

% Butterworth Bandpass filter designed using FDESIGN.BANDPASS.

% All frequency values are in Hz.
Fs = 48000;  % Sampling Frequency

Fstop1 = 3;           % First Stopband Frequency
Fpass1 = 4;           % First Passband Frequency
Fpass2 = 7;           % Second Passband Frequency
Fstop2 = 8;           % Second Stopband Frequency
Astop1 = 80;          % First Stopband Attenuation (dB)
Apass  = 1;           % Passband Ripple (dB)
Astop2 = 80;          % Second Stopband Attenuation (dB)
match  = 'stopband';  % Band to match exactly

% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.bandpass(Fstop1, Fpass1, Fpass2, Fstop2, Astop1, Apass, ...
                      Astop2, Fs);
Hd = design(h, 'butter', 'MatchExactly', match);

Есть предложения, как мне решить проблему?

Спасибо всем

4 голоса | спросил Ricardo Frederico Leote Mota 5 Maypm11 2011, 21:39:55

2 ответа


0

Один более простой подход может состоять в том, чтобы просто взять БПФ и обнулить частотные компоненты, отличные от конкретного диапазона, который вас может заинтересовать, а затем инвертировать БПФ, чтобы вернуться во временную область.

Имейте в виду, что вам придется обнулять положительную частоту и отрицательную частоту, чтобы убедитесь, что сигнал в частотной области сопряжен симметрично относительно частоты 0. Если вы этого не сделаете, вы получите сложный сигнал после вычисления обратного БПФ.

EDIT: Например, следующий код создает две синусоиды во временной области, соответствующий ДПФ (вычисленный с помощью БПФ), а затем один из удаленных пиков.

t = 0:0.01:0.999;
x = sin(t*2*pi*4) + cos(t*2*pi*8);
subplot(2,2,1);
plot(x)
title('time domain')
subplot(2,2,2);
xf = fft(x);
plot(abs(xf))
title('frequency domain');
subplot(2,2,3);
xf(9) = 0; xf(93) = 0;  % manual removal of the higher frequency
plot(abs(xf));
title('freq. domain (higher frequency removed)');
subplot(2,2,4);
plot(ifft(xf));
title('Time domain (with one frequency removed)')

Несколько вещей на заметку. Частотная область в DFT имеет несколько различных диапазонов: смещение постоянного тока (постоянное значение), которое является частотой 0; положительный частотный диапазон, в котором (для длины N исходного сигнала) находятся записи от 1 до N /2; отрицательный частотный диапазон, который представляет собой записи от N /2 до N-1; Обратите внимание, что это не опечатка, самая высокая частота (та, что в N /2) дублируется и является одинаковым значением как для положительных, так и для отрицательных частот. (Некоторые люди используют fftshift, чтобы показать, как это может нарисовать человек, но на самом деле это просто для взгляда /понимания.)

Что касается частоты, которую нужно удалить, вам придется выяснить это самостоятельно, но я могу дать вам подсказку. Самая высокая частота (в частотном положении N /2) будет соответствовать самой высокой частоте, представляемой вашей системой, то есть fs /2, где fs - ваша частота дискретизации. Вы можете масштабировать соответственно, чтобы выяснить, какие из них следует отрицать.

Если вы не отрицаете соответствующие отрицательные частоты правильно, вы получите воображаемый сигнал при обратном fft.

Последний комментарий. Этот метод будет работать только в том случае, если вы можете позволить себе весь свой сигнал раньше времени, поскольку вам нужно использовать ДПФ для всего сигнала. Если вы хотите сделать это в реальном времени, вам нужно будет создать какой-то фильтр, как вы делали ранее.

ответил Chris A. 6 Mayam11 2011, 01:39:36
0

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

ответил Hephaestus 5 Maypm11 2011, 21:59:23

Похожие вопросы

Популярные теги

security × 330linux × 316macos × 2827 × 268performance × 244command-line × 241sql-server × 235joomla-3.x × 222java × 189c++ × 186windows × 180cisco × 168bash × 158c# × 142gmail × 139arduino-uno × 139javascript × 134ssh × 133seo × 132mysql × 132