GSR (Galvanic Skin Response) 전기피부반응 신호를 분석하여 교감신경계의 신호를 해석하기 위해 사용되는 방법입니다.
해당 방법은 비음수 희소 역합성을 활용하여 GSR 신호를 저주파 성분과 고주파 성분으로 분리합니다.
분리된 신호에서 특정 신호 처리 및 분석 작업을 수행하며 이를 수행하는 코드를 분석합니다.
if (length(signalIn)/sr < 80)
display('Signal not enought large. May be longer than 80 senconds')
end
if (sum(isnan(signalIn))>0)
display('Signal contains NaN')
end
입력 신호의 길이가 최소한 80초 이상인지를 체크합니다.
신호에 NAN 값이 포함되어 있는지 또한 함께 체크합니다.
signalAdd = [signalIn(1)*ones(20*sr,1); signalIn; signalIn(end)*ones(60*sr,1)];
signalIn 에 추가적인 데이터를 덧붙이는 작업을 수행합니다.
해당 코드는 경계 효과를 줄이기 위해 사용이 됩니다.
signaIn 은 GSR 신호로, 길이가 N인 신호입니다.
이 신호의 시작과 끝에 각각 20초, 60초 분량의 데이터를 붙입니다.
signalIn(1) 은 입력 신호의 첫 번째 값인 상수로 ones(20*sr, 1)은 길이가 20*sr이며 모든 값이 1인 벡터입니다.
벡터에 signaIn(1) 을 곱하게 되면, 첫 번째 값이 20초동안 지속해서 나타난 꼴을 하게 됩니다.
마찬가지로 signalIn(end) 인 끝 값을 이용하여 60초 분량의 벡터를 생성합니다.
signalIn = [3; 5; 7];
sr = 2; % 샘플레이트(Hz)
signalAdd = [3; 3; ... (40개); 3; 5; 7; 7; 7; ... (120개)]
예시데이터를 이용한다면 위와 같은 모양이 될 것입니다.
if (sr > 8)
[P,Q] = rat(8/sr);
signalAdd = resample(signalAdd,P,Q);
sr = 8; % 샘플 레이트를 8Hz로 고정
end
입력 신호의 샘플 레이트가 8Hz 보다 높다면, 이를 8Hz으로 리샘플링 합니다.
8과 sr의 비율을 가장 근사한 분수 형태로 표현합니다.
그리고 P와 Q는 이 분수의 분자와 분모가 됩니다.
예를들어 sr이 16이라면, P는 1 Q는 2가 됩니다.
y = resample(x, P, Q)
기존 샘플레이트에 (P/Q)를 곱한 값을 새로운 샘플레이트로 신호를 구성합니다.
P > Q 경우, 고주파 신호가 추가되며 데이터 포인트가 증가합니다.
P < Q 경우, 저주파 필터링 후 일부 샘플이 제거됩니다.
다음의 코드는 Toeplitz 행렬을 사용하여 고주파 성분(SCR)을 모델링 합니다.
여러 주파수에서 신호를 분석하기 위한 다중 분석 기법을 구성합니다.
Rzeros = zeros(N+L, Lreg);
Rzeros = zeros(N+L,Lreg);
srF = sr * [0.5 0.75 1 1.25 1.5]; % Multianalysis
for j = 1:length(srF)
t_rf = 0:(1/srF(j)):10; % 10초 동안의 반응 모델링
rf_biexp = exp(-t_rf/taus(2)) - exp(-t_rf/taus(1));
rf_est = taus(3) * rf_biexp;
rf_est = rf_est / sqrt((sum(rf_est.^2))); % 표준화
Rzeros(1:N, 1 + (j-1)*Lreg : j*Lreg) = ...
toeplitz([rf_est zeros(1, N-length(rf_est))], zeros(1, Lreg));
end
다음에서는 Toeplitz 행렬을 저장할 공간을 초기화합니다.
N은 70초 동안의 신호 길이.
L은 추가적인 신호 길이 (10초)
Lreg 는 모델링할 60초의 신호 길이를 뜻합니다.
srF = sr * [0.5 0.75 1 1.25 1.5]; % Multianalysis
기본 샘플 레이트의 절반, 75%, 100% … 를 고려하여 분석하겠다는 정보입니다.
신호를 다양한 샘플 레이트로 분석하기 위함입니다.
for j = 1:length(srF)
t_rf = 0:(1/srF(j)):10; % 10초 동안의 반응 모델링
t_rf 는 10초 동안의 시간 축을 정의합니다.
시간의 간격은 샘플링 주파수(srF(j))에 따라 결정이 됩니다.
10초 동안의 신호 반응을 모델링하기 위한 시간 벡터 생성이 목적인 코드입니다.
rf_biexp = exp(-t_rf/taus(2)) - exp(-t_rf/taus(1));
쌍지수 감쇠 함수입니다.
GSR 신호의 고주파 성분은 자극에 대한 급격한 방응과 점진적인 감소로 나타납니다.
자극이 발생하면 피부의 전도도가 빠르게 증가하고, 이후 교감 신경의 활성도가 낮아지며 전도도는 서서히 감소합니다.
이러한 특성을 수학적으로 잘 나타내는 모델이 해당하는 쌍지수 감쇠 함수입니다.
rf_est = taus(3) * rf_biexp;
함수에 특정한 스케일링 계수를 곱하여 신호를 특정한 크기로 조정합니다.
rf_est = rf_est / sqrt((sum(rf_est.^2))); % 표준화
sum(rf_est.^2) 는 신호의 에너지를 계산하는 식입니다.
sqrt 를 통해 에너지의 제곱근을 계산하게 됩니다.
신호를 에너지로 나누어, 정규화된 신호를 생성합니다.
결과적으로 정규화된 신호는 에너지가 1로 유지됩니다.
Rzeros(1:N, 1 + (j-1)*Lreg : j*Lreg) = ...
toeplitz([rf_est zeros(1, N-length(rf_est))], zeros(1, Lreg));
입력 신호의 시간적 관계를 모델링 합니다.
행렬의 각 행은 신호 벡터의 오른쪽으로 이동된 버전을 포함합니다.
행렬의 열이 1 + (j-1)*Lreg : j*Lreg 범위에 해당합니다.
즉, 다중 분석(j)의 각 srF 에 대해 서로 다른 Toeplitz 행렬이 생성됩니다.
toeplitz([a, b, c], [0, 0, 0])
a 0 0
b a 0
c b a
0 c b
0 0 c
데이터의 첫 번째 열은 고주파 성분의 값과 그 뒤에 패딩이 붙어지는 형태입니다.
이어서 열을 보면 하나씩 점점 밀리는데, 시간의 흐름에 따라서 신호의 고주파 성분이 위치하는 곳을 나타냅니다.
[beta,~, activationHist,~,~,~] = lasso(R, signalCutIn, sr, Kmax, epsilon);
signalEst = R * beta + b0;
위 코드에서 R이 앞서 toeplitz 통해 변환한 행렬입니다.
signalCutln 은 실제 입력 신호입니다.
Kmax 는 Lasso 알고리즘의 최대 반복 횟수.
epsilon 은 최적화 종료 기준입니다.
다음은 위 코드에서 사용되는 lasso 함수의 수식입니다.
def lasso_regression(X, y, lambda_val, max_iter=1000, tol=1e-6):
"""
Simple Lasso Regression implementation using coordinate descent.
X: Input feature matrix
y: Target vector
lambda_val: Regularization parameter
max_iter: Maximum number of iterations
tol: Convergence tolerance
"""
n_samples, n_features = X.shape
beta = np.zeros(n_features) # Initialize coefficients to 0
for iteration in range(max_iter):
beta_old = beta.copy()
for j in range(n_features):
# Compute the residual excluding feature j
residual = y - (X @ beta - X[:, j] * beta[j])
# Update beta[j] using the soft-thresholding rule
rho = X[:, j].T @ residual
if rho < -lambda_val:
beta[j] = (rho + lambda_val) / (X[:, j] @ X[:, j])
elif rho > lambda_val:
beta[j] = (rho - lambda_val) / (X[:, j] @ X[:, j])
else:
beta[j] = 0 # Set beta[j] to 0 if |rho| <= lambda_val
# Check for convergence
if np.sum(np.abs(beta - beta_old)) < tol:
break
return beta
다음은 residual 을 구하는 방식을 수식을 나타낸 것입니다.
y는 실제 데이터의 값이죠.
그리고 이 값에서 X@beta 가 모델이 예측한 결과값입니다.
이때 X@beta 에서 j번째에 해당하는 feature을 제외하고 실제 신호에서 감산합니다.
residual 은 현재 모델에서 j번째 신호를 제외한 예측값을 나타냅니다.
잔여신호를 뜻하게 됩니다.
실제 신호에서 모델이 예측한 부분을 빼고 남는 값이기 때문입니다.
이어서 코드 내에서 rho를 구하는 수식입니다.
내적을 사용합니다. 내적을 통해 j번째 신호와 잔여 신호 간의 상관성을 측정합니다.
두 벡터가 잘 정렬되어 상관성이 높을수록 내적 값이 커집니다.
#j 번째 신호
X_j = [1, 2, 3]
#잔여신호
residual = [4, 5, 6]
#내적 계산 output
output = (1*4) + (2*5) + (3*6)
#output 은 32
위와 같이 두 신호가 비슷하게 정렬되어 있다면, output 은 크게 나타납니다.
#j 번째 신호
X_j = [1, −1, 1]
#잔여신호
residual = [4, 5, 6]
#내적 계산 output
output = (1*4) + (-1*5) + (1*6)
#output 은 5
그런데 한쪽 신호를 바꿔서 두 신호의 상관성을 낮춘다면
결과는 5로 기존에 상관성이 컸을 때의 결과보다 작아집니다.
#j 번째 신호
X_j = [-1, -2, -3]
#잔여신호
residual = [4, 5, 6]
#내적 계산 output
output = (-1*4) + (-2*5) + (-3*6)
#output 은 -32
이때 다음과 같이 음의 상관성이 큰 경우도 있을 것입니다.
이러한 경우들을 코드에서는 lambda_val 에 따라서 beta[j]를 업데이트 합니다.
for j in range(n_features):
# Compute the residual excluding feature j
residual = y - (X @ beta - X[:, j] * beta[j])
# Update beta[j] using the soft-thresholding rule
rho = X[:, j].T @ residual
if rho < -lambda_val:
beta[j] = (rho + lambda_val) / (X[:, j] @ X[:, j])
elif rho > lambda_val:
beta[j] = (rho - lambda_val) / (X[:, j] @ X[:, j])
else:
beta[j] = 0 # Set beta[j] to 0 if |rho| <= lambda_val
이러한 경우들을 코드에서는 Soft Thresholding 규칙에 따라 처리합니다.
양의 상관이 큰 경우에는 과적합 방지를 위해서 신호의 기여를 줄여줍니다.
정의 방향으로 진행되는 여기서는 -통해 신호의 기여가 줍니다.
반대로, 음의 상관이 큰 경우는 -방향으로 설명이 되니 이곳읠 신호를 축소하기 위해서는 +lambda를 사용합니다.
lamda 갑이 엄청난게 커버리면, 대부분으 beat[j] 는 0으로 업데이트가 이뤄집니다.
실제로 lamda 값이 커지면 커질수록 신호가 살아남는 Feature의 개수는 매우 줄어듭니다.
signalEst = R * beta + b0;
R은 앞서 Toeplitz 행렬로 생성한 시간적 특성과 고주파 성분을 반영한 행렬입니다.
beta는 Lasso 알고리즘으로 각 열에 대한 기여도를 나타내는 가중치 벡터입니다.
두 행렬과 벡터를 곱하면 각 시간 이동된 고주파 성분과 가중치를 조합하여 전체 신호를 복원 가능합니다.
SCL = R(:,1:6)*beta(1:6) + b0;
저주파 성분을 계산하는 데에 쓰이는 방식입니다.
앞선 SignalEst 를 만드는 수학적 공식이 그대로며
쓰이는 데이터만 다릅니다.
Toeplitz 행렬의 첫 6개의 열을 사용합니다.
가중치 벡터의 6번째 값까지 사용합니다.
SCRline = beta(7:end);
SCRaux = zeros(Lreg, 5); % 초기화
SCRaux(:) = SCRline;
driver = sum(SCRaux', 1); % 고주파 성분 합산
SCRline: [ 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26]
SCRaux:
[[ 7. 8. 9. 10. 11.]
[12. 13. 14. 15. 16.]
[17. 18. 19. 20. 21.]
[22. 23. 24. 25. 26.]]
SCRaux`:
[[ 7. 12. 17. 22.]
[ 8. 13. 18. 23.]
[ 9. 14. 19. 24.]
[10. 15. 20. 25.]
[11. 16. 21. 26.]]
Driver: [58. 62. 66. 70. 74.]
고주파 성분은 Toepliz 행렬의 나머지 열과 베타 값으로 계산합니다.
Lreg 는 60초 동안의 샘플 수를 뜻합니다.
결과적으로..
SCL : 저주파 성분
SCR: 고주파 성분
이 됩니다.
ind = find(SCRaux > 0); % SCR 값이 양수인 인덱스만 선택
scr_temp = SCRaux(ind); % 양수 값 추출
[scr_ord, ind2] = sort(scr_temp, 'descend'); % 크기 순으로 정렬
scr_fin = scr_ord(1); % 가장 큰 값 선택
ind_fin = ind(ind2(1)); % 해당 인덱스 저장
% 최소 거리 조건을 만족하는 값만 선택
for i = 2:length(ind2)
if all(abs(ind(ind2(i)) - ind_fin) >= dmin)
scr_fin = [scr_fin scr_ord(i)];
ind_fin = [ind_fin ind(ind2(i))];
end
end
driver = zeros(size(SCRaux)); % 최종 SCR 벡터 초기화
driver(ind_fin) = scr_fin; % 선택된 값만 유지
SCR값 중 0보다 큰 값을 추출합니다.
반응의 강도 순으로 내림차순을 수행합니다.
자극 간의 최소 거리를 만족하지 못하는 항목을 제거하고,
선택된 항목만을 최종적으로 driver 에 저장합니다.
scr_max = scr_fin(1); % SCR 최대값
threshold = rho * scr_max; % 임계값 설정 (rho 비율로 계산)
driver(driver < threshold) = 0; % 임계값 이하인 항목 제거
scr_fin 에서 가장 큰 값을 기준으로 임계값을 계산합니다.
임계값은 최대값에 특정 비율을 곱하여 설정합니다.
임계값 이하인 반응을 노이즈로 간주하여 제거합니다.
여기까지 GSR 데이터를 처리하고 분석하는 Sparse EDA 알고리즘을 분석하였습니다.