[Python] 스펙트럼 분석 (spectral analysis, spectrum analysis)
Python 분석과 프로그래밍/Python 통계분석 2021. 10. 3. 18:44이번 포스팅에서는 주기성을 가지는 센서 시계열 데이터, 음향 데이터, 신호 처리 데이터 등에서 많이 사용되는 스펙트럼 분석 (spectral analysis) 에 대해서 소개하겠습니다.
많은 시계열 데이터는 주기적인 패턴 (periodic behavior) 을 보입니다. 스펙트럼 분석 은 시계열 데이터에 내재되어 있는 주기성(periodicity) 을 찾아내는 분석기법으로서, spectral analysis, spectrum analysis, frequency domain analysis 등의 이름으로 불립니다.
스펙트럼 분석을 할 때 우리는 먼저 시계열 데이터를 시간 영역(time domain)에서 주파수 영역(frequency domain)로 변환을 합니다. 시계열의 공분산은 스펙트럼 확률밀도함수 (spectral density function) 로 알려진 함수로 재표현될 수 있습니다. 스펙트럼 밀도함수는 시계열과 다른 주파수의 사인/코사인 파동 간의 제곱 상관계수로 구하는 주기도(periodogram) 를 사용하여 추정할 수 있습니다. (Venables & Ripley, 2002)
아래의 그림은 시계열(time series)과 스펙트럼(spectrum, frequency) 간의 관계를 직관적으로 이해할 수 있게 도와줍니다.
[ Time Series vs. Spectrum/ Frequency, Fourier Transform(FT) vs. Inverse FT ]
참고로, 주파수(frequency)는 단위 시간 당 반복되는 이벤트가 발생하는 회수를 의미합니다. 주파수를 측정하는데 사용되는 단위인 Hz (Hertz, 헤르쯔) 는 전자기파의 존재를 처음으로 증명한 Heinrich Rudolf Hertz 의 이름에서 따온 것인데요, 가령 2 Hz 는 1초에 2번의 반복 사이클이 발생하는 시계열의 경우, 4 Hz 는 1초에 4번의 반복 사이클이 발행하는 시계열의 주파수에 해당합니다. (아래 그림 예시 참조)
주기성을 띠는 시계열 데이터나 공간 데이터로 부터 주파수 (spectrum, frequency) 를 계산하거나, 반대로 주파수로 부터 원래의 데이터로 변환할 때 고속 푸리에 변환 (Fast Fourier Transform, 줄여서 FFT) 알고리즘을 많이 사용합니다. 시계열을 주파수로 변환하는 과정 자체를 Fourier Transform 이라고 하며, 이산 푸리에 변환 (Discrete Fourier Transform, DFT) 은 시계열 값을 다른 주파수의 구성요소들로 분해함으로써 구할 수 있습니다. 지금까지 가장 널리 사용되는 고속 푸리에 변환 (FFT) 알고리즘은 Cooley-Tukey 알고리즘 입니다. 이것은 분할-정복 알고리즘 (divide-and-conqer algorithm) 으로서, 반복적으로(recursively) 원래의 시계열 데이터를 훨씬 적은 개수의 주파수로 분할합니다.
[ Fast Fourier Transform(FFT) vs. Inverse FFT]
이제 Python 의 spectrum 모듈을 사용해서 스펙트럼 분석을 해보겠습니다. 터미널에서 $ pip install spectrum 으로 먼저 spectrum 모듈을 설치해줍니다.
## installment of spectrum python library at a terminal
$ pip install spectrum
spectrum 모듈에서 Periodogramm data_cosine 메소드를 불러옵니다.
예제로 사용할 샘플데이터로서, data_cosine() 메소드를 사용해서 200 Hz 주파수 (freq=200) 의 주기성을 띠는 코사인 신호(cosine signal)에 백색 잡음 (amplitude 0.1) 도 포함되어 있는 이산형 샘플 데이터 1,024개 (N=1024, sampling=1024) 를 생성해보겠습니다.
## -- importing modules
from spectrum import Periodogram, data_cosine
## generating a toy data sets
## data contains a cosine signal with a frequency of 200Hz
## buried in white noise (amplitude 0.1).
## The data has a length of N=1024 and the sampling is 1024Hz.
data = data_cosine(N=1024, A=0.1, sampling=1024, freq=200)
## plotting the generated data
plt.rcParams['figure.figsize'] = [18, 8]
plt.plot(data)
plt.show()
위의 시계열 그래프가 1,024개 데이터 포인트가 모두 포함되다보니 코사인 신호의 주기성이 눈에 잘 안보이는데요, 앞 부분의 50개만 가져다가 시도표로 그려보면 아래처럼 코사인 파동 모양의 주기성을 띠는 데이터임을 눈으로 확인할 수 있습니다.
plt.rcParams['figure.figsize'] = [18, 8]
plt.plot(data[:50])
plt.show()
spectrum 모듈의 Periodogram() 메소드에 위에서 생성한 샘플 데이터를 적용해서 periodogram 객체를 만든 후에 run() 함수를 실행시켜서 스펙트럼을 추정할 수 있고, p.plot() 으로 간단하게 추정한 스펙트럼을 시각화할 수 있습니다.
## the Power Spectrum Estimation method provided in spectrum.
## creating an object Periodogram.
p = Periodogram(data, sampling=1024)
## running the spectrum estimation
p.run() # or p()
p.plot(marker='o')
스펙트럼 분석 결과 중 개수(N), 주파수의 세기 (power of spectrum dendity) 속성값을 확인할 수 있습니다.
## some information are stored and can be retrieved later on.
p.N
[Out] 1024
p.psd
# [Out] array([1.06569152e-04, 1.74100855e-03, 3.93335837e-03, 1.10485152e-03,
# 8.28499278e-03, 7.64848301e-04, 2.94332240e-03, 7.35029918e-03,
# 1.21273520e-03, 3.72342384e-05, 1.87058542e-03, 4.35817479e-03,
# 6.60803332e-04, 3.01094523e-03, 3.14914528e-03, 3.40040266e-03,
# 1.70112710e-04, 4.04579597e-04, 1.73102721e-03, 5.17048610e-03,
# 5.89537834e-03, 3.43004928e-03, 2.60894164e-03, 2.29556960e-03,
# 6.23836725e-03, 7.35588289e-03, 3.88793797e-03, 2.25691304e-03,
# ... 중간 생략함 ...
# 2.57852156e-03, 1.23483225e-03, 1.49139300e-03, 9.09082081e-04,
# 5.24999741e-04])
우리가 위에서 data_cosine(N=1024, A=0.1, sampling=1024, freq=200) 메소드를 사용해서 200 Hz 의 데이터를 생성하였기 때문에, 이 샘플 데이터로부터 주파수를 분해하면 200 Hz 가 되어야 할텐데요, 스펙트럼 분석 결과도 역시 200 Hz 에서 스파이크가 엄청 높게 튀는군요.
plt.rcParams['figure.figsize'] = [10, 6]
plt.plot(p.psd)
plt.show()
* R 을 이용한 스펙트럼 분석은 https://rfriend.tistory.com/64 를 참고하시기 바랍니다.
이번 포스팅이 많은 도움이 되었기를 바랍니다.
행복한 데이터과학자 되세요! :-)
