Python에서 수행해 본 간단한 FFT 코드

일요일 아침(이 글은 평일에 예약 발행되겠지만)이네요.. 오늘 아침은 꽤 상쾌하고 약간 몽롱한.. 뭐 아무튼 기분이 좋아지는 아침이네요^^. 요즘은 뭔가를 마무리하는 단계에서 오는 급급하게 일을 처리하는 몇몇 모습들을 제 주변 동료들이 보여주고 있는데요. 다들 그러지 말고 R&D에서 오는 error와 그걸 debug하는 그냥 일상적인 과정이라고 생각했으면 하는데요.. 뭐 정작 저도 잘 되지 않으니 문제입니다만.. 뭐 아무튼 역시 웃으면 다 좋아진다는 건 확실합니다^^.

[바로가기]에서 저는 원래 그 글의 목적이 함수의 입력을 키워드를 이용한다는 것이 주제였는데... 거기서 다룬 예제인 sine 함수 그리기를 이용해서 [바로가기]에서 클래스까지 설명을 했네요. 이번에는 이왕 삼각함수를 다룬김에 FFT(겁나게 빠른 푸리에 변환이라는 뜻의 Fast Fourier Transform [바로가기])을 수행하는 코드를 예제로 다뤄볼까 합니다.

항상 이렇게 예제를 다루다 보면 전혀 관계없는 부분이 좀 걸리는데요. 이번의 경우는 그래프로 예쁘게 표현하기 였네요..ㅠㅠ. 그래서 검색하다가 괜찮은 예제코드가 보이길래 그걸 참조합니다.[바로라기] 근데 제가 참조한 이 사이트는 원래 각종 plot을 보여주는 곳인데 소스코드 단을 이용해서 그래프 결과가 나오는 것 같은데요. 꽤 ... 재미있더라구요. 언제 한번 들어가서 놀아야겠어요^^

import numpy as np
import matplotlib.pyplot as plt
import sinClass

test1 = sinClass.sinWaveForm(amp = 1, freq=1, endTime = 5)
test2 = sinClass.sinWaveForm(amp = 2, freq=5, endTime = 5)
test3 = sinClass.sinWaveForm(amp = 4, freq=10, endTime = 5)

t = test1.calcDomain()
resultTest1 = test1.calcSinValue(t)
resultTest2 = test2.calcSinValue(t)
resultTest3 = test3.calcSinValue(t)

Ts = test1.sampleTime 					# sampling interval
Fs = 1/Ts 						# sampling rate
t = test1.calcDomain()					# time vector

12번행까지는 [바로가기]에서 다룬 부분인데요. 테스트를 위한 시험 신호를 만들기 위해 test1, test2, test3까지를 생성했구요. 각각 주파수를 1, 5, 10Hz로 설정해서 그 결과를 resultTest1,2,3에 받았습니다. 흠... 9번과 16번에 t를 계산하는 코드가 두번 들어갔네요.(ㅠㅠ) 뭐 아무튼.. 시간을 t에 저장하구요. samplTime은 0.01초로 설정되어 있는 것인데 그게 Ts입니다. 그 역수를 샘플 주파수 Fs로 둡니다.

y = resultTest1 + resultTest2 + resultTest3

n = len(y) 					# length of the signal
k = np.arange(n)
T = n/Fs
freq = k/T 					# two sides frequency range
freq = freq[range(int(n/2))] 			# one side frequency range

Y = np.fft.fft(y)/n 				# fft computing and normalization
Y = Y[range(int(n/2))]

이제 FFT 예제에 적용해 볼 시험 신호를 만들는게 y입니다. y는 resultTest1,2,3을 모두 더했습니다. 결국 FFT한 결과에서 1, 5, 10Hz가 결과로 나와야만 되는 거죠^^. 위 코드에서는 y의 길이를 잡고(n) 거기서 주파수영역을 정하고 오늘의 결론인 FFT는 간편히 numpy(np)의 fft를 이용한 것 뿐입니다. 아.. 그걸 다시 n으로 나눠준 것은 그렇지 않으면 몇 백... 혹은 아주 큰 값이 y축에 나타나서 그래프를 보기 좀 그래서 간단히 nomalize를 해준 겁니다. 이제 이걸 잘 표현하면 되죠^^

fig, ax = plt.subplots(2, 1)
ax[0].plot(t, y)
ax[0].set_xlabel('Time')
ax[0].set_ylabel('Amplitude')
ax[0].grid(True)
ax[1].plot(freq, abs(Y), 'r', linestyle=' ', marker='^') 
ax[1].set_xlabel('Freq (Hz)')
ax[1].set_ylabel('|Y(freq)|')
ax[1].vlines(freq, [0], abs(Y))
ax[1].grid(True)
plt.show()

위 코드는.. 음.. Python 초보인 저한테는 꽤 재미있는 코드였는데... 뭐 다른 분들은 아닐 수도 있구요^^. 아무튼 한 plot에 두 그림을 나란히 놓고 비교하기 위해 subplot을 사용하는 겁니다. 첫 번째는 ax[0] 두번때는 ax[1]입니다. 각각 xlabel과 ylabel을 정했구요. 6번행에서 linestyle을 공란으로 두고 marker를 위를 향한 삼각형을 했습니다. 그리고, vlines라는 함수를 이용해서 라인을 바로 그런거죠. 흠.... 이게 MATLAB에서는 명령 한 줄인데... 그게 있을 거라고 생각하고 검색했는데 없어서 저렇게 아이디어를 내고 찾아서 구현했습니다.ㅠㅠ. 좀 더 간편한 방법이 있으면 알려주세요^^

아무튼 위 그림이 그 결과 입니다.위의 시간 영역에서의 그래프를 FFT를 수행해보니 1, 5, 10Hz의 성분이 0.5, 1, 2의 상대적 크기로 나타났다고 되어 있네요. 제가 각각 시험신호를 만들때 그 크기를 1, 2, 4로 amp를 잡았는데 그것도 잘 반영되어 나타났습니다.^^. 물론 저렇게 시뮬레이션 시간이 5초로 좀 충분하면 사실 어느 정도는 눈으로도 주파수 성분을 알 수 있습니다만....

이렇게 짧은 시간으로 바꾸면 저 주파수 성분 세 개가 모두 보이지는 않습니다.^^. 뭐 아무튼 이렇게 FFT를 수행하는 그것도 numpy를 이용해서 간편히 수행하는 것과, 그걸 잘 그래프로 표현한 예제를 다루었네요^^. 전체 코드는 GitHUB에 올려 두었습니다. [바로가기]