AMUSE noise preprocessing 을 이용한 EEG frequency analysis

EEG signal processing with AMUSE 알고리즘

AMUSE라는 알고리즘을 사용하면,

신호에 있는 noise를 제거하여 원신호를 복구할 수 있다고 합니다. 

 

Matlab을 사용하였으며, 특별히 Matlab에 라이브 스크립트를 사용했습니다.

 

AMUSE 함수는 ICALAB을 다운받으신 후 압축을 풀고, 환경 변수 설정에서 압축 푼 폴더를 저장하면 amuse 함수를 사용할 수 있습니다.

 

또한, 신호를 그래프로 확인하기 위하여 다음 매트랩 제공 FFT 예제를 사용했습니다.

https://kr.mathworks.com/help/matlab/ref/fft.html

고속 푸리에 변환(Fast Fourier Transform) - MATLAB fft - MathWorks 한국

코드는 다음과 같이 작성하였습니다.

load('AD.mat')
eeg = c_outeeg.data;
Fs = 1000;          %샘플링 개수(Fs)는 1초에 1000개 입니다.
L = 4000;           %총 샘플링 수(L)는 4초로 4000개 입니다.
T = 1/Fs;           %주기(T) = 1sec/1000 = 0.001, 즉 0.001초 당 1개를 샘플링합니다.
t = (0:L-1)*T;      %time domain 그래프를 보기 위해서 t의 범위를 0.001초 단위로,
                    %매트랩에서 변수 = (a:b)의 의미는 크기 a~b인 배열을 말합니다.
                    %예를 들어 지금 우리의 수식에선, t = (0:3999)*0.001으로 설정되어 있습니다.
                    %따라서 배열(행렬) t는 범위로 0 ~ 3.999를 가지고, 0.001씩 증가하는 요소를 가지게 됩니다.
                    

sum_time_signal = zeros(32, 4000);                  %zeros(32,4000)은 32x4000 영행렬을 만들어주는 함수입니다.
for people = (1:70)                                 %사람마다 AMUSE의 결과가 다르므로, '사람'을 뽑아주는 for문을 만듭니다.
    data = eeg(:, :, people);                       %처음 eeg배열의 크기는 32(채널)x4000(데이터)x70(사람)이기 때문에,
                                                    %32(채널)x4000(데이터)만 불러올 수 있게 됩니다.
    system = amuse(data);                           %data에 AMUSE 알고리즘을 적용하여 linear system을 예측합니다.
    noisy = system * data;                          %선형예측도로 정렬시킨(그러나 노이즈가 포함된) data를 얻어냅니다.
    free = [noisy(1:10, :); zeros(22, 4000)];       %data에서 선형예측도가 좋은 10개를 남기고, 나머지를 영행렬로 바꾸어 노이즈를 제거합니다.
    time_signal = inv(system) * free;               %다시 계측신호로 돌려 놓아, 노이즈가 제거된 time domain singal을 얻습니다.
    
    sum_time_signal = sum_time_signal + time_signal;%이 신호들을 계속 누적시켜 더해갑니다.
end                                                 %그래서 위에서 sum_time_signal의 크기를 32x4000으로 해놓은 것입니다.
                                                    %기존에 있는 행렬에 누적시키려면, size가더하는 행렬의 size와 같아야하기 때문입니다.
                                                    
avg_time_signal = sum_time_signal/70;               %우리는 사람70명을 1명으로 평균 내야하기 때문에, 70으로 나누어줍니다.

avg_freq_signal = [];                               %time domain signal을 frequency domain signal로 저장할 행렬을 만들어 줍니다.
for channel = (1:32)                                %지금 avg_time_signal의 크기는 32(채널)x4000(데이터)이므로,
                                                    %각 채널 별로 나누어 푸리에 변환을
                                                    %시키기 위한 for문을 만들어 줍니다.
    Y = fft(avg_time_signal(channel, :));           %푸리에 변환은 fft, fast fourier transform 함수를 사용하면 바로 바뀝니다.
    P2 = abs(Y/L);                                  %변환된 신호 Y의 요소 전체는 (a + jb)인 복소수입니다.
                                                    %복소수 (a + jb)를 absolute value(절댓값)하면 root(a제곱 + b제곱)이 됩니다. 
                                                    % 즉, abs를 취하면 magnitude가 됩니다.
                                                 
    P1 = P2(1:L/2 + 1);                             %P2는 양방향 스펙트럼이기 때문에 주파수분석을 위한 단방향 스펙트럼P1을 구합니다.
    P1(2:end - 1) = 2*P1(2:end - 1);
    avg_freq_signal = [avg_freq_signal; P1];        %배열(행렬)avg_freq_signal 에 차곡차곡 쌓습니다.
end
f = Fs*(0:(L/2))/L;
%이렇게 time domain에서 70명을 평균냈고, frequency domain으로 변환하여 단방향 진폭 스펙트럼을 구했습니다.
f1 = figure;
subplot(2,1,1);
plot(t, avg_time_signal(1, :));title('Time domain'); xlabel('t (ms)'); ylabel('X(t)');
subplot(2,1,2);
plot(f, avg_freq_signal(1, :)); title('Frequency domain'); xlabel('f (Hz)'); ylabel('magnitude');

그래프는 다음과 같습니다.

 

AD 치료 전

 

AD 치료 후