Power Spectrum in IDL

오늘은 IDL에서 파워 스펙트럼(Power Spectrum)을 구하는 간단한 방법을 소개하고 관련 예제를 살펴보기로 하겠습니다. IDL에서는 원래 파워 스펙트럼을 구하는 전통적인 방법이 있습니다. 바로 FFT라는 내장함수를 사용하여 FFT(Fast Fourier Transform)를 수행하고 이 결과를 제곱하여 파워 스펙트럼을 산출하는 방식입니다. 하지만 오늘은 이 방법 대신 좀 더 간편한 방법을 소개하고자 하는데요. FFT_POWERSPECTRUM이라는 다소 긴 이름의 내장함수를 사용하는 방법입니다. 단, 이 방법을 사용하는데 있어서는 약간의 제한이 있습니다.

1) 대상 자료가 1차원 배열이어야 한다.

2) IDL 8.4부터 사용 가능하다.

이 점들을 염두에 두고 일단 시작해보겠습니다. 먼저 예제로 사용할 1차원 데이터를 만들어야 하는데요. 시간에 따라 어느 정도의 주기성을 보이는 데이터를 사용해야 합니다. 그래서 이러한 가상의 데이터를 하나 만들어보겠습니다. 그 과정은 다음과 같습니다.

dt = 1/1000.

t = [0:4:dt]

freq1 = 2.1

data = COS(2*!pi*t*freq1)

data += RANDOMN(seed, t.length)/3

여기서는 1/1000초 단위로 측정된 4000개의 데이터 포인트들이 있는데, 시간에 따른 변동성은 대략 2.1 정도의 주파수를 갖고 있다고 가정하였습니다. 즉, 1초마다 2.1 주기의 변동을 보이는 데이터라고 가정하였습니다. 그리고 RANDOMN 함수를 사용하여 약간의 노이즈에 해당되는 난수값들을 더해줌으로써 잡음에 해당되는 고주파수의 변동성 패턴도 추가하였습니다. 이 원본 데이터를 다음과 같이 표출해보면 그 모습은 다음 그림과 같습니다. 이 그림을 보면 이러한 변동성의 패턴이 잘 보입니다.

win1 = WINDOW(DIMENSIONS=[800, 400])

p = PLOT(t, data, XTITLE='Time', YTITLE='Signal', COLOR='dodger blue', /CURRENT)

이제 이러한 데이터에 대하여 파워 스펙트럼을 구하고자 합니다. 이를 위하여 FFT_POWERSPECTRUM 함수를 다음과 같이 사용하였습니다.

f = FFT_POWERSPECTRUM(data, dt, FREQ=freq)

여기서 첫번째 인자는 원본 데이터 배열이고, 두번째 인자는 원본 데이터의 시간 간격에 해당되는 값입니다. 이 값을 명시하지 않을 경우에는 무조건 1로 간주가 되는데, 여기서는 1/1000이라고 가정했었기 때문에 이 값을 부여하였습니다. 그리고 FREQ라는 키워드를 사용하여 freq라는 배열을 돌려받게 되는데, 이것은 나중에 파워 스펙트럼을 그래프로 나타낼 때 X축에 위치할 주파수 값들입니다. 일단 이 결과를 다음과 같은 방법으로 표출하면 그 모습은 다음 그림과 같습니다.

win2 = WINDOW(DIMENSIONS=[800, 400])

pf = PLOT(freq, f, /YLOG, XTITLE='Frequency (Hz)', YTITLE='Power', $

  COLOR='orange red', /CURRENT)

이 그림에서 X축은 헤르쯔(Hz) 단위의 주파수이고 Y축은 파워(Power)에 해당됩니다. 이 그림을 보면 원본 데이터의 시간적 변동 패턴에 대한 주파수를 찾아낼 수 있습니다. X축상에서 앞부분에서 잘 보이는 것 같으므로, 이 부분만 자세히 들여다보기 위하여 다음과 같은 방식으로 X축의 범위를 제한한 그림을 또 그려보았습니다.

win3 = WINDOW(DIMENSIONS=[800, 400])

pf = PLOT(freq, f, /YLOG, XTITLE='Frequency (Hz)', YTITLE='Power', $

  COLOR='orange red', XRANGE=[0, 100], /CURRENT)

이 그림을 보면 파워가 큰 주파수의 값이 2.1Hz 근처임을 쉽게 확인할 수 있습니다. 물론 작은 파워에 해당되는 높은 주파수 값들도 상당히 많은데 바로 노이즈(Noise)에 해당됩니다. 이와 같이 FFT_POWERSPECTRUM 함수를 사용하면 변동 패턴에 존재하는 주파수들을 비교적 쉽게 확인할 수 있습니다.

앞서 사용했던 예제 데이터는 대세가 되는 변동 패턴의 주파수가 한 개인 경우인데, 이러한 주파수가 더 많은 경우에도 파워 스펙트럼을 산출하면 쉽게 파악이 가능합니다. 앞선 예제에서 예제 데이터를 생성하는 부분을 다음과 같이 수정해 봅시다.

dt = 1/1000.

t = [0:4:dt]

freq1 = 2.1

freq2 = 8.4

data = COS(2*!pi*t*freq1) + sin(2*!pi*t*freq2)

data += RANDOMN(seed, t.length)/3

여기서는 기존의 2.1외에 8.4라는 주파수에 해당되는 패턴까지 추가하여 새로운 데이터를 만들어 보았습니다. 나머지 그림을 표출하는 방법은 위와 동일합니다. 이 원본 데이터의 모습은 다음 그림과 같습니다.

그림으로 봐도 대세가 되는 주파수가 하나만 있지는 않다는 것이 보입니다. 이 데이터에 대하여 똑같은 방식으로 파워 스펙트럼을 구해보면 그 결과는 다음 그림들과 같습니다.

그림에서 볼 수 있듯이, 2.1과 8.4에 해당되는 주파수가 가장 지배적인 주파수들로서 쉽게 확인이 됩니다. 따라서 이와 같이 파워 스펙트럼을 잘 활용하면, 다양한 주기성을 갖고 변화하는 데이터들의 변동 패턴을 쉽게 알아내는 것이 가능합니다. 다음 회에서는 이러한 방법을 실전 데이터에 응용한 사례를 소개해보기로 하겠습니다.