R, Python 분석과 프로그래밍의 친구 (by R Friend)

데이터 변환 방법으로서

(1) 표준화

(2) 정규분포화

(3) 범주화

     - 이산형화

     - 이항변수화

(4) 개수 축소

(5) 차원 축소

     - 주성분분석

     - 요인분석 

(6) 시그널 데이터 변환

   - 푸리에 변환 (FFT: Fast Fourier Transform)

    - 웨이블릿 변환(Wavelet Transform)

의 6개 구분 중에서

전파, 진동, 소리, 파도, 빛 등 시간(time domain)에 따라 주기성(periodicity)을 띠면서 파형을 형성하는 데이터를 주파수(spectrum domain) 대역별로 세기로 변환하는 푸리에 변환(Fourier Transform)에 대해서 알아보겠습니다.  참고로, 푸리에(Jean-Baptiste Joseph Fourier, 1768~1830)는 프랑스의 수학자이자 물리학자로서, 푸리에 변환은 바로 이분의 이름을 딴 것이랍니다.  저 아래에 애니메이션의 오른쪽 하단에 사진이 푸리에가 되겠습니다.

푸리에 변환은 물리나 공업 분야에서 폭넓게 사용되고 있기에, 특히 시그널 데이터, 기계 데이터, 센서 데이터 등의 주기성을 띤 데이터를 분석하는 분이라면 반드시 알아야 할 변환이라고 하겠습니다.

푸리에 변환을 이해하려면 삼각함수(특히 사인함수, 코사인함수), 미적분, 함수의 사칙연산, 함수의 직교성, 푸리에 급수와 계수 등에 대해서 알아야 하는데요, 어려운 수학 공식은 다른 교재를 참고하시기 바라며, 이번 포스팅에서는 직관적으로 이해할 수 있는 그림과 예시를 들어서 가급적 쉽게 설명을 해보겠습니다.

아래의 그림처럼 시간에 따른 진폭 데이터를 => 주파수별 세기의 데이터로 변환하는 작업을 푸리에변환이라고 하며, 그 반대를 푸리에역변환이라고 합니다.

[ 푸리에변환과 푸리에역변환 ]

먼저, 일정한 괘도를 회전하는 운동을 하는 시간함수로 나타내는 법에 대해서 알아보겠습니다.  아래 애니메이션의 왼쪽이 일정한 속도로 반지름이 일정한 괘도를 회전하는 운정이 되겠구요, 오른쪽이 이 회전운동을 시간축에 옮겨놓았을 때의 모양입니다.  전형적인 코사인(cosine) 형태를 띠고 있습니다.

* 출처: http://www.di.fc.ul.pt/~jpn/r/fourier/fourier.html

아래는 왼쪽의 반지름 거리와 시간 주기가 다른 4개의 회전운동을 오른쪽에 시간 축에 진폭을 나타낸 그래프가 되겠습니다.  아래 보는 것처럼 회전운동은 시간을 축으로 해서 진폭이 변화하는 값을 사인(sine) 또는 코사인(cosine) 함수로 나타낼 수 있습니다.

* 출처: http://www.di.fc.ul.pt/~jpn/r/fourier/fourier.html

이렇게 주기성을 띤 회전운동을 시간함수로 나타낼 수 있는데요, 이 시간함수는 사실 여러개의 주파수를 띤 시간함수들이 합해진 것입니다.  아래 예시로 든 그림에서는 주기(주파수)가 다른 3개의 시간함수가 함쳐져서 1개의 시간함수를 형성하고 있는데요, 주파수별로 필터링을 해서(주파수 성분을 구한다고 함) 세기가 큰 (peaks) 주파수를 헤아리면 되겠습니다. 시간함수를 알면 주파수 스펙트럼을 구할 수 있고, 주파수 스펙트럼을 알면 이들을 합쳐서 시간함수를 구할 수 있게 됩니다.

[ 시간함수와 주파수 스펙트럼의 관계 ]

* 그림 출처: aragec.com556

주파수는 1초에 파동 cycle이 몇 번 반복되느냐를 나타내는 말로서, 단위는 Hz(헤르츠)를 사용합니다.  아래의 3개의 파형을 예로 들면, 3개 파형 모두 진폭은 -1 ~ +1 로 동일한 반면에 주기는 모두 다릅니다. (즉, 주파수가 모두 다름)  첫번째 파형은 1초에 2회 주기이므로 주파수는 2Hz, 두번째 파형은 1초에 4회 주기이므로 4Hz, 세번째 파형은 1초에 6회 주기이므로 6Hz 주파수가 되겠습니다.

[ 주기와 진폭 ]

참고로, 악기 음 조율할 때 사용하는 소리굽쇠는 440Hz 의 '라'음을 낸답니다. 1초에 주기가 440회 진동이 발생한다는 뜻입니다. 서울/경기 지역의 '별이 빛나는 밤에' 라디오 주파수가 95.9kHz 인데요, 이는 전파가 1초에 9만 5천9백번 진동한다는 의미입니다.  음악에서는 주파수가 낮을 수록(현이 굷고 길수록) 저음이 나고, 주파수가 높을 수록(현이 가늘고 짧을 수록) 고음이 납니다.

이제 R을 가지고 1초에 2pi 만큼을 단위 시간 구간으로 해서 진폭과 주기를 달리한 4개의 사인함수 그래프도 그려보고, FFT (Fast Fourier Transform) 변환 실습을 해보도록 하겠습니다.

4개의 사인함수를 각각 순서대로 그려보면 아래와 같습니다. 첫번째와 두번째 그래프는 진폭(높이, y축)이 '2'로서 동일하고, 세번째와 네번째 그래프는 진폭이 '5'로서 동일합니다.  1초에 몇번의 주기가 있는지, 즉 주파수에 해당하는 주기는 순서대로 1, 2, 3, 7로서 뒤로 갈수록 증가하는 사인함수 그래프로서, 주파수가 커질수록 1초당 주기 갯수가 많아집니다 (진동 회수가 증가).

> # 사인함수 시간에 따른 그래프 > par(mfrow = c(1,4)) > plot(x, signal.1, type='l', ylim=c(-5,5)); abline(h=0, lty=3) # 진폭 2 & 주기 1인 사인함수 > plot(x, signal.2, type='l', ylim=c(-5,5)); abline(h=0, lty=3) # 진폭 2 & 주기 2인 사인함수 > plot(x, signal.3, type='l', ylim=c(-5,5)); abline(h=0, lty=3) # 진폭 5 & 주기 3인 사인함수 > plot(x, signal.4, type='l', ylim=c(-5,5)); abline(h=0, lty=3) # 진폭 5 & 주기 7인 사인함수

다음으로, 위 4개의 시간에 따른 사인함수를 합한 후에 그래프로 나타내보겠습니다.  위 4개의 개별 시간에 따른 사인함수가 주기성을 띠므로 아래의 1개로 합쳐진 시간함수도 일정한 주기성을 띠게 됩니다.  (푸리에 변환은 지금 하는 작업과는 거꾸로, 여러개의 시간함수들이 합쳐진 시간함수를 개별 시간함수들로 분해 해서 각 개별 시간함수들의 주파수 성분을 구하는 것입니다.  푸리에 변환할 때는 파형이 주기성을 띤다고 가정하고 변환을 진행합니다.) 

> # 사인함수 4개 합치기 (sine function summation) > signal.1234 <- signal.1 + signal.2 + signal.3 + signal.4 > head(signal.1234, n=30) [1] 0.000000 1.749660 3.442051 5.022470 6.441235 7.655888 8.633062 9.349913 9.795059 9.968964 [11] 9.883774 9.562593 9.038247 8.351592 7.549453 6.682284 5.801674 4.957812 4.197038 3.559593 [21] 3.077684 2.773944 2.660382 2.737852 2.996074 3.414214 3.961977 4.601179 5.287703 5.973757 > > # 사인함수 4개 합친 그래프 그리기 > par(mfrow = c(1,1)) > plot(x, signal.1234, type='l', main = "Sum of signal.1&2&3&4", + xlab = "Time", ylab = "Amplitude") > abline(h=0, lty=3)

R로 푸리에 변환할 때는 stats 패키지의 fft() 함수를 사용합니다.  fft()함수를 적용한 값을 관측치 개수의 반개((N-1)/2 = (201-1)/2 = 100 )로 나누어서 표준화를 시켜주고, abs()함수를 적용해 절대값을 취하게 됩니다.  plot()으로 그래프를 그려보면 좌우로 대칭인 그래프가 그려지는데요, 절반을 기점으로 해서 똑같은 대칭 그래프라서 왼쪽의 값(그래프)만을 사용하면 되겠습니다.  

이 예제에서는 주파수가 1Hz, 2Hz, 3Hz, 7Hz 짜리 4개의 사인함수를 더한 것이었으므로 1Hz, 2Hz,  3Hz, 7Hz 지점에서 피크(peak)를 치는 스펙트럼을 보여줄 겁니다.  따라서 1~20Hz 까지만 뽑아서 자세히 본 그래프가 두번째 그래프가 되겠습니다. 

> # 푸리에 변환 (FFT: Fast Fourier Transform) > library(stats) > N <- length(x) > > fft_x_abs <- abs(fft(signal.1234)/((N-1)/2)) # 표준화, 절대값 > plot(fft_x_abs, type="h") >   > # 앞의 주파수 20개만 그래프 그려보기 > plot(fft_x_abs[1:20], type="h")

참고로, Python 을 이용한 스펙트럼 분석은 https://rfriend.tistory.com/690 를 참고하시기 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감 ~♡' 단추를 꾸욱 눌러주세요.^^