[EEG] Signal Filtering in Python
Fast Fourier Tranform 및 filter 설명
Python을 사용하여 신호처리에서 사용되는 filtering을 알아보고자 합니다.
Filter
신호 처리에서 필터(Filter)란 특정한 신호에서 원하지 않는 신호를 차단하거나 원하는 신호만 통과시키는 과정을 의미합니다.
필터의 종류는 다양하며, 그 중에서 band-pass filter 및 band-stop filter을 구현하고자 합니다.
Types of filter(Source: Wikipedia).
Fast Fouier Transform
시간 영역(time domain)의 신호를 주파수 영역(frequency domain)으로 변환하기 위해 고속 푸리에 변환(fast fourier transform, FFT)를 사용합니다.
View of a signal in the time and frequency domain(Source: NTi Audio)).
푸리에 변환(Fourier transform)의 수식적 이해를 얻고 싶다면 Dark programmer님의 Fourier transform의 이해와 활용을 참고하시길 바랍니다.
예제로는 시계열 데이터인 EEG data를 사용하겠습니다.
An eeg signal in the time domain
Single-channel EEG 데이터는 4초간 측정되었으며, 250 Hz의 sampling frequency가 사용되었습니다.
Python에서는 주로 scipy 모듈을 활용하여 푸리에 변환을 하며, tensorflow를 사용하여 tensor를 조작할 때에는 tensor의 값을 변경하지 못한다는 제약이 있기 때문에 tensorflow 내장 함수를 사용하고 있습니다.
따라서 scipy 및 tensorflow를 활용하여 FFT를 계산하고 시각화하고자 합니다.
from scipy import fft
def fft_scipy(x, sfreq): # sfreq = sampling frequency
# FFT
X = fft.fft(x) / x.shape[-1]
# FFT magnitude
fft_magnitude = abs(X)
# FFT shift
fft_shift = fft.fftshift(fft_magnitude)
freqs = fft.fftshift(fft.fftfreq(x.shape[-1], d=1 / sfreq))
return fft_shift, freqs
import tensorflow as tf
def fft_tensorflow(x, sfreq):
x = tf.cast(x, dtype=tf.complex64)
# FFT
X = tf.signal.fft(x) / x.shape[-1]
# FFT magnitude
fft_magnitude = abs(X)
# FFT shift
fft_shift = tf.signal.fftshift(fft_magnitude)
freqs = tf.signal.fftshift(fft.fftfreq(x.shape[-1], d=1 / sfreq))
return fft_shift, freqs
1. FFT
tensorflow의 FFT는 데이터 타입이 float일 경우 에러가 납니다.
“TypeError: Value passed to parameter ‘input’ has DataType float64 not in list of allowed values: complex64, complex128”
먼저 데이터 타입을 복소수(complex number)로 변환해줍니다.
다음으로 scipy.fft.fft 또는 tf.signal.fft를 사용하여 고속 푸리에 변환(FFT)을 계산할 수 있으며, 하단의 그림과 같은 결과를 확인하실 수 있습니다.
Double Sided FFT - with normalization
또한, FFT 결과값을 신호의 길이인 sample 개수로 나누어주어야 normalization 되면서 제대로 된 값을 얻을 수 있습니다. Normalize 하지 않을 경우 다음과 같이 y 값이 높은 것을 확인하실 수 있습니다.
Double Sided FFT - without normalization
2. FFT magnitude
푸리에 변환을 통해 두 종류의 그래프를 얻을 수 있습니다.
- 푸리에 스펙트럼 그래프(Fourier spectrum graph): 주파수에 따른 magnitude 그래프
- 위상 그래프(Phase graph): 주파수에 따른 angle 그래프
푸리에 스펙트럼 그래프를 얻기 위해서는 FFT의 output인 복소수값에 절대값을 씌워야 합니다. python의 내장함수 abs를 통해 쉽게 절대값을 구할 수 있습니다.
Double Sided FFT - without FFTShift
3. FFT shift
Nyquist Theorem에 따라 원래 실수(real) 신호의 최대 주파수 f_max = sfreq / 2가 되며, 푸리에 변환을 하면 -f_max에서 +f_max까지의 주파수 성분이 나옵니다.
scipy.fft.fftfreq 함수를 통해 주파수 구간을 확인할 수 있습니다.
sfreq = 250
freqs = scipy.fft.fftfreq(signal.shape[-1], d=1/sfreq)
print(freqs)
# output
[ 0. 0.24975025 0.4995005 0.74925075 0.999001
1.24875125 1.4985015 1.74825175 1.998002 2.24775225
...
122.37762238 122.62737263 122.87712288 123.12687313 123.37662338
123.62637363 123.87612388 124.12587413 124.37562438 124.62537463
124.87512488 -124.87512488 -124.62537463 -124.37562438 -124.12587413
-123.87612388 -123.62637363 -123.37662338 -123.12687313 -122.87712288
-122.62737263 -122.37762238 -122.12787213 -121.87812188 -121.62837163
...
-2.4975025 -2.24775225 -1.998002 -1.74825175 -1.4985015
-1.24875125 -0.999001 -0.74925075 -0.4995005 -0.24975025]
sfreq는 sampling frequency로써, 시간 축에서 몇 초에 한번 data를 획득했냐를 의미합니다. 예제 데이터의 경우 0.004초에 한번 샘플링 했으므로 sfreq=250이며, 최대 주파수는 125Hz가 됩니다.
또한 output을 보면 0부터 등장하는데요, FFT를 하게 되면 양의 주파수 축이 먼저 나오고 다음으로 음의 최대 주파수부터 나오기 때문입니다.
우리에게 익숙한 x축의 형태로 변환하기 위해서 scipy.fft.fftshift 또는 tf.signal.fftshift를 통해 shift 해야 합니다.
sfreq = 250
freqs = scipy.fft.fftfreq(signal.shape[-1], d=1/sfreq)
print(scipy.fft.fftshift(freqs))
# output
[-124.87512488 -124.62537463 -124.37562438 -124.12587413 -123.87612388
-123.62637363 -123.37662338 -123.12687313 -122.87712288 -122.62737263
...
-2.4975025 -2.24775225 -1.998002 -1.74825175 -1.4985015
-1.24875125 -0.999001 -0.74925075 -0.4995005 -0.24975025
0. 0.24975025 0.4995005 0.74925075 0.999001
1.24875125 1.4985015 1.74825175 1.998002 2.24775225
2.4975025 2.74725275 2.997003 3.24675325 3.4965035
...
122.62737263 122.87712288 123.12687313 123.37662338 123.62637363
123.87612388 124.12587413 124.37562438 124.62537463 124.87512488]
FFT의 성질에 의해, 주파수 성분에서의 -f_max와 +f_max(혹은 -sfreq/2 ~ +sfreq/2) 신호 이외에는 동일한 신호 성분이 반복됩니다.
Double Sided FFT with FFTShift
음의 주파수 부분을 생략하면 최종적으로 다음과 같이 나타낼 수 있습니다.
y, x = fft_scipy(signal, sfreq=250) # signal shape: [1001, ]
# y, x = fft_tensorflow(signal, sfreq=250)
onesided = signal.shape[-1] // 2
plt.plot(x[onesided:], y[onesided:])
plt.title("Fast Fourier Transform", fontsize=15)
plt.xlabel('Frequency (Hz)', fontsize=15)
plt.ylabel('Amplitude', fontsize=15)
plt.ylim([0, 1.5])
plt.tight_layout()
plt.show()
One Sided FFT
Band-pass filter
대역 통과 필터(band-pass filter)는 특정 주파수 대역만을 통과하고 다른 주파수 대역은 차단하는 역할을 합니다.
Band-pass filter frequency response(Source: Analog Devices).
from scipy.signal import butter, filtfilt, lfilter
def butter_bandpass_filter(signal, sfreq, lowcut, highcut, order=5):
nyq = 0.5 * sfreq
low = lowcut / nyq
high = highcut / nyq
b, a = butter(order, [low, high], btype='band')
y = filtfilt(b, a, signal)
return y
예제로 사용한 데이터에 20 - 30 Hz까지 band-pass filter를 걸어보면 다음과 같은 그래프가 나타납니다.
FFT Band-pass filter
Band-pass filterining으로 20 - 30 Hz의 주파수만 남아있는 것을 확인하실 수 있습니다.
Filter 함수로는 scipy.signal.lfilter 또는 scipy.signal.filtfilt가 사용됩니다.
두 개의 차이는 필터링 된 신호가 지연되느냐에 있습니다.
어떤 필터든 간에 신호가 통과하는데 아주 작게라도 시간이 걸릴 것입니다. 결국 그런 물리적인 크기가 전기적 파장에 영향을 주어 시간축에서는 신호위상의 지연이 발생하게 됩니다. 이러한 현상을 군지연(Group Dealy)라고 합니다. scipy.signal.lfilter는 필터링에 따른 신호의 지연을 반영하지만 scipy.signal.filtfilt는 원래 신호에서 발생하는 시점과 동일한 시간으로 유지합니다.
- lfilter: 실제 전자 필터와 유사한 순방향 필터링 입니다.
- filtfilt: 필터링할 때 신호를 이동시키지 않는 영위상 필터링(zero-phase filtering) 입니다. 실시간(online)으로는 사용할 수 없으며, 신호 획득 이후의 후처리(offline)에서만 사용 가능합니다.
Comparison of the Scipy function filtfilt and lfilter(Source: Endolith)
Band-stop filter
대역 저지 필터(band-stop filter)는 다른 말로 notch filter 또는 band-reject filter라고 하며, 특정한 주파수 대역만을 차단합니다. Band-pass filter와 반대의 역할을 하는 것이죠.
Band-stop filter frequency response(Source: Analog Devices).
def butter_bandstop_filter(signal, sfreq, lowcut, highcut, order=5):
nyq = 0.5 * sfreq
low = lowcut / nyq
high = highcut / nyq
i, u = butter(order, [low, high], btype='bandstop')
y = filtfilt(i, u, signal)
return y
예제로 사용한 데이터에 20 - 30 Hz까지 band-stop filter를 걸어보면 다음과 같은 그래프가 나타납니다.
FFT Band-stop filter
Band-pass filter와 정반대로 작동하며, 20 - 30 Hz를 제외하고는 신호가 남아있는 것을 알 수 있습니다.
이로써, scipy 및 tensorflow를 활용하여 각 내장함수 안에서 FFT가 어떻게 이루어지는지 살펴보았으며 scipy로 band-pass filter 및 band-stop filter를 구현해보았습니다.
References
[1] Fast Fourier Transform Code
[2] Fourier Transform Code
[3] lfilter vs filtfilt
[4] Group delay
[5] Band-stop filter Code