BCI/1주 BCI 기초 다지기

[1주 BCI 기초 다지기] 6일차: 파워 스펙트럼 분석, PSD

_\oyo/_ 2025. 1. 28. 21:57

6일차: 특징 추출 및 분석 (25-01-28)

1. Power Spectral Density (PSD) 개념

PSD란?

신호의 주파수 성분별 에너지 분포를 나타내는 분석 방법.

신호가 각 주파수 대역에서 얼마나 많은 파워(에너지)를 가지는지를 정량적으로 보여줌.

  • 단위: 보통 μV^2/Hz 또는 V^2/Hz.

PSD 사용 이유

  • 주파수 대역별 신경 활동 분석
    • EEG의 Delta, Theta, Alpha, Beta, Gamma 대역에서의 각 대역의 에너지 크기를 분석.
  • 신호의 주파수 특성 시각화
    • 신호가 특정 주파수에서 강한 성분을 가지는지 확인 가능.
    • 노이즈나 특정 패턴 감지 가능.
  • 주파수 대역의 특징 추출
    • 머신러닝 학습 떄 사용되는 특징으로 활용 가능.

PSD 계산 방법

여러 가지 방법이 있지만, EEG에서는 보통Welch’s Method를 사용함.

이는 FFT 기반으로 신호를 여러 구간으로 나누고, 각 구간에서 평균을 계산해 파워 스펙트럼을 안정적으로 추정할 수 있는 방법.

2. MATLAB에서의 실습

오늘 실습에서는 “GigaDB(https://gigadb.org/dataset/view/id/100542/File_page/15)의 sess01_subj24_EEG_ERP.mat”를 사용하고자 한다.

데이터 불러오기 및 탐색

% 데이터 로드
load('sess01_subj24_EEG_ERP.mat');

% 파일에 포함된 변수 확인
whos;

해당 코드를 통해 데이터를 불러온 후 확인해 보면, 주요 변수들이 각각 포함된 것을 알 수 있다.

변수명으로 미루어 보아, EEG 신호 데이터는 EEG_ERP_test와 EEG_ERP_train 구조체로 저장되어 있으므로 해당 변수를 더블클릭 해, 이 구조체 안의 데이터 필드를 확인하고 탐색해 보자.

변수명과 형식을 기반으로 구조체 필드를 분석해 보면,

  • smt (800x1980x62)
    • 62개 채널의 EEG 데이터(800개의 trial x 1980개의 시간 포인트 x 62개의 채널).
    • 이 필드가 우리가 사용할 EEG 신호 데이터일 가능성이 높음.
  • x (969760x62)
    • 더 큰 크기의 데이터로, 아마도 원본 데이터 전체에 해당할 것.
  • t (1x1980)
    • 시간 벡터로 보이며, 1980개의 샘플이 어떤 시간에 해당하는지 나타냄.
  • fs (1000)
    • 샘플링 주파수. 1초에 1000개의 샘플을 수집했음을 의미.
  • chan (1x62 cell)
    • 각 채널의 이름(Fp1, Fp2 등)이 저장된 셀 배열.
  • y_dec, y_logic, y_class
    • 이벤트 정보 또는 자극 조건을 나타내는 필드로 추정.
  • EMG 및 관련 필드
    • 근전도(EMG) 관련 정보.

위와 같을 것이라 유추할 수 있다.

여기서 우리가 사용할 데이터는 EEG신호인 smt, 샘플링 주파수인 fs, 채널 이름인 chan, 시간 정보인 t이다.

% 주요 데이터 추출
eeg_data = EEG_ERP_test.smt;   % EEG 신호 데이터
fs = EEG_ERP_test.fs;          % 샘플링 주파수
channel_names = EEG_ERP_test.chan; % 채널 이름
time = EEG_ERP_test.t;  % 시간 벡터

데이터 구조를 익히기 위해, 특정 채널의 신호를 시각화하는 방법을 알아보자.

eeg_data는 3차원(샘플 x 시간 x 채널) 배열이므로, 예시로 Fp1 채널의 데이터를 아래와 같이 시각화 할 수 있다.

% 특정 채널의 신호 시각화
channel_idx = 1; % 시각화할 채널 번호
signal = squeeze(eeg_data(:, :, channel_idx)); % 채널 데이터 추출

figure;
plot(time, mean(signal, 1)); % 모든 샘플에 대한 평균 신호
title(['채널 ', channel_names{channel_idx}, '의 평균 신호']);
xlabel('Time (ms)');
ylabel('Amplitude (\\muV)');

Power Spectral Density 계산 및 시각화

이제 Welch's Method를 사용하여 PSD를 계산해 보자. pwelch 함수를 활용해 PSD를 계산하고, 이를 통해 주파수 대역별 에너지를 시각화 할 것이다.

  • [psd, f] = pwelch(분석할 신호, 각 구간에 적용할 윈도우 함수, 각 구간의 중첩 길이, FFT 길이, 샘플링 주파수);

이때, 각 구간에 적용할 윈도우 함수의 경우 구간의 경계에서 신호가 급격히 변화하면 FFT 결과에 누설(Leakage)이 생길 수 있다.

이를 방지하기 위해 Hamming Window를 사용해 신호의 경계를 부드럽게 만들어 사용하자.

해당 함수는 hamming(샘플 수)로 불러올 수 있다.

원래 PSD는 선형 스케일로 계산되는데, 해당 주파수 스펙트럼을 dB 스케일로 변환하여 표현하면 아래와 같은 장점이 있다.

  • 작은 값과 큰 값을 시각적으로 비교하기 쉬움.
  • 신호 대 잡음비(SNR)를 분석하거나 신호 강도를 로그 스케일로 보기 편함.

따라서 PSD는 dB 스케일로 변환하여 시각화 해보자.

해당 내용을 종합하여, 아래와 같이 코드를 짤 수 있다.

% 채널 1 데이터 추출
channel_idx = 1;
signal = squeeze(eeg_data(:, :, channel_idx)); % 채널 데이터 추출
signal_mean = mean(signal, 1); % 모든 샘플에 대한 평균 신호

% PSD 계산 (Welch's Method)
[psd, f] = pwelch(signal_mean, hamming(256), 128, 256, fs);

% PSD 시각화
figure;
plot(f, 10*log10(psd)); % dB 스케일로 표현
title(['채널 ', channel_names{channel_idx}, '의 Power Spectral Density']);
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');

 

그래프에서 x축은 주파수 대역이 Hz 단위로, y축은 각 주파수 대역의 에너지가 dB 단위로 나타난다.

낮은 주파수에서 높은 에너지가 나타나는 형태를 확인할 수 있는데, 이는 일반적인 EEG 신호의 특성으로 Delta~Alpha 대역에서 에너지가 많이 분포하는 경향이 있다.

저작자표시 비영리 변경금지