목차
Python을 사용한 지 약 2년이 좀 지난 거 같다. 기계공학을 전공한 나로서는 아직도 최적의 프로그램 코드 작성이 아직 버겁다. 최근에는 현장에서의 dataset을 이용한 데이터 분석 및 이상감지 알고리즘에 대한 연구개발을 진행 중이다. 이상감지를 위해서는 데이터의 주파수 분석을 해야 하고, 현장에서의 데이터 양이 많기 때문에 한 개의 주기만 추출하여 관리할 필요가 있다. 구글링을 통해 해당내용을 아무리 검색해 보아도 주기 추출은 찾아볼 수 없었고, FFT변환에 대해서도 코드 설명 및 사용방법에 대해서는 잘 다뤄지고 있었으나, 코드 안의 기능에 대해서는 쉽게 나와있는 자료는 없었다. 지금까지 여러 자료조사 경험을 통해 얻는 기능들에 대해 정리해 보기로 했다.
1. 주기 데이터란?
일반적으로 주기를 갖는 데이터의 양상은 아래 그림과 같다. 단순해 보이긴 하지만 FFT 변환을 통해 해당 데이터가 얼마나 많은 주파수를 갖고 있는지 확인할 수 있다.
데이터를 분석해 보니 모든 주파수가 실제 데이터를 표기하는 것으로 확인됐다. 그림 2에서 거의 모든 frequency에서 amplitude 값이 존재하는 것처럼 보인다. 신호공학에서 이러한 부분에 대해 잘 설명되어 있지만, 여기서는 실제 주파수인 1X, 2X들의 값만을 다루기로 한다.
주기 특성을 갖는 데이터는 크기와 주기가 있는 여러 개의 신호들의 총합으로 나타낼 수 있다. 공업수학 및 신호처리에 대한 자세한 설명은 생략한다. 단순하게 아래 수식만 참고하자. 최근에는 python과 같은 프로그램 라이브러리에서 많은 함수들을 제공하고 있어, 대략적인 개요만 이해하고 있다면 해당 함수를 활용해서 데이터 분석을 진행할 수 있다. 아래 수식을 단순하게 신호 생성 함수라 하자.
2. 진동데이터 생성하기
Python에서 Numpy 라이브러리의 FFT를 이해하기 위해 위의 수식을 이용해서 진동데이터를 생성해 보겠다. 우선 필요한 수식 함수를 갖고 있는 라이브러리를 위해 numpy와 그래프 라이브러리인 matplotlib을 코드 첫 줄에서 추가했다.
import numpy as np
import matplotlib.pyplot as plt
현실성 있는 데이터를 만들기 위해 저주파에서 가장 큰 진폭(amplitude)을 대입하였고, 고주파 대역으로 갈수록 점점 작게 설정하겠다. 실제 데이터들을 다뤄보면 가장 큰 진폭을 갖는 주파수가 1X 주파수임을 확인할 수 있다. (1X 주파수 : 회전 주파수로 한번 회전 시 특성을 갖는 주파수다.) 주파수를 이용하여 한 주기를 추출하는 내용은 맨 아래쪽에 기술되어 있다.
신호생성 코드에서 총 4개의 복합 신호를 이용하였다. amplitude는 2, 1, 0.5, 0.2로 점차 작아지게, 주기는 주파수는 10, 20, 30, 40으로 커지게 설정하였다. Fs는 sampling frequency로 1초에 1000개의 데이터가 생성된다고 생각하면 되겠다. 만일 측정되는 데이터라면 1초에 1000개의 데이터를 습득한다고 생각하자. T는 주기로 Fs의 역수다. 수식 1을 풀어서 코드로 작성하면 아래와 같이 쓸 수 있다.
Fs = 1000
T = 1/Fs
end_time = 1
time = np.linspace(0, end_time, Fs)
amp = [2, 1, 0.5, 0.2]
freq = [10, 20, 30, 40]
signal_1 = amp[0]*np.sin(freq[0]*2*np.pi*time)
signal_2 = amp[1]*np.sin(freq[1]*2*np.pi*time)
signal_3 = amp[2]*np.sin(freq[2]*2*np.pi*time)
signal_4 = amp[3]*np.sin(freq[3]*2*np.pi*time)
signal = signal_1 + signal_2 + signal_3 + signal_4
plt.plot(time, signal)
plt.show()
신호생성 코드 작성 시 for문을 사용하면 amp와 freq의 개수에 따라 좀 더 많은 신호 조합을 이루는 결과물을 얻을 수 있다.
signal_list = []
for i in range(len(amp)):
signal_list.append(amp[i]*np.sin(freq[i]*2*np.pi*time))
signal = sum(signal_list)
3. FFT 변환
numpy의 FFT 함수를 이용하여 만들어진 신호를 주파수 도메인으로 변환시켜 보자(FFT 변환)
s_fft = np.fft.fft(signal) # 추후 IFFT를 위해 abs를 취하지 않은 값을 저장한다.
amplitude = abs(s_fft)*(2/len(s_fft)) # 2/len(s)을 곱해줘서 원래의 amp를 구한다.
frequency = np.fft.fftfreq(len(s_fft), T)
plt.xlim(0, 50)
plt.stem(frequency, amplitude)
plt.grid(True)
plt.show()
단 두줄로 FFT변환이 완료되었다. np.fft.fft(signal) 함수로 데이터를 FFT로 변환한 후 abs(s_fft) 함수로 변환된 값의 절대치를 만들었다. fft() 함수를 이용해 주파수 도메인으로 변환하면 실수와 복소수로 구성된 값으로 나오기 때문에 abs를 통해 절대치로 바꿔주면 된다.
abs(s_fft) x (2/len(s_fft)) 와 같이 뒤쪽에 2/len(s_fft)를 곱해줬는데, 실제 amplitude를 구하기 위해서 scaling을 해준 것이다. np.fft.fftfreq(len(s_fft), T) 함수는 x축을 만들어 주기 위한 것으로 데이터의 개수는 신호의 개수(len(s_fft))와 일치시켜줘야 하고, 정확한 주파수 값을 얻기 위해서는 T(주기)를 입력해 주어야 한다.
변환된 그래프에서 amplitude와 frequency가 초기 입력한 값과 동일함을 알 수 있다. 10Hz의 amplitude의 값이 2를 나타내고 있으며, 10Hz 바로 양 옆을 보면 0이 아님(그림 4에서 빨간색 점선 동그라미)을 확인할 수 있다. 그림 2에서 X1, X2에 대한 설명을 잠시 했던 것과 같이 그림 4의 그래프에서 X1은 10Hz가 되고 2X는 20Hz가 될 것이다. 즉, 1X, 2X 등의 값들 이외에는 변환과정에서 부수적으로 발생되는 값이라 생각하면 된다. IFFT(inverse fft)를 수행할 때 해당 값들을 제외한 모든 값을 0으로 만드는 이유다.
4. 주기 추출 (1X)
이제 그림 3에서의 데이터에서 한 주기만 추출해 보도록 하자. 위쪽에서 설명한 것과 1X는 회전 주파수라 했다. 회전 주파수는 한 바퀴의 회전 성분을 의미하기 때문에 1X의 주파수를 이용하면 된다. 일반적인 데이터에서는 1X는 가장 큰 진폭을 갖는 주파수를 구하면 된다.
fft_freq = frequency.copy()
peak_index = amplitude[:int(len(amplitude)/2)].argsort()[-1]
peak_freq = fft_freq[peak_index]
주파수 값을 변경하기 위해 fft_freq로 우선 복사를 했다. FFT 변환을 하면 사실상 양의 구역과 음의 구역 2개의 구역으로 데이터가 반환되기 때문에 amplitude[:int(len(amplitude)/2]와 같이 양의 값만을 이용하여 진폭에 대한 각각의 index를 최대값으로부터 argsort()를 이용하여 정렬했다. [-1]을 사용하면 최대값이 반환된다. peak_index에는 최대 진폭의 인덱스가 입력되었으니 이 인텍스로 최대 진폭을 갖는 주파수를 확인한다
fft_1x = s_fft.copy()
fft_1x[fft_freq!=peak_freq] = 0
filtered_data = 2*np.fft.ifft(fft_1x)
cycle = round(Fs/peak_freq)
ifft를 이용하여 fft값을 인버스 하는 과정이다. 주의할 점은 fft변환 후 절대값(abs)을 취한 값이 아닌 fft변환 직후의 값을 인버스 해야 한다. s_fft 변수를 따로 만든 이유다. fft_1x[fft_freq != peak_freq] = 0 에서와 같이 peak_freq를 제외한 모든 주파수를 0으로 만든 후 ifft를 통해 인버스를 한다. 이때, scaling을 위해 2를 곱해 주었다. 한 주기의 데이터 개수를 구하기 위해 sampling frequency에서 peak frequency를 나눠 주었다.
한 cycle을 수집하기 위한 데이터 개수 =
sample frequency / 회전 주파수(fft에서 peak frequency) =
[단위] 데이터개수 / 1 cycle = (데이터 개수 / 1 sec) / (1 cycle / 1 sec)
plt.subplot(2, 1, 1)
plt.title('1X sin wave')
plt.plot(filtered_data[:400], marker='o', color='darkgreen', alpha=0.3)
plt.subplot(2, 1, 2)
plt.title('1-period graph')
plt.plot(signal[:400], marker='o', color='lightgrey')
plt.plot(signal[:cycle], color='indigo')
plt.subplots_adjust(left=0.125, bottom=0.1, right=0.9, top=0.9, wspace=0.2, hspace=0.35)
plt.show()
최종 확인을 위해 가시화 코드를 구성했다. 전체 데이터에서 400개까지만 도표화한다. subplot을 이용해서 첫 번째는 filtering 된 sin파 그래프를 확인하고, 두 번째에서는 실제 데이터에서 지금까지 구한 cycle이 진짜 한 주기만 나타내는지 체크한다.
5. 마치며
관련된 업무를 하는 도중에 한 주기를 추출하기 위해 필요한 방법은 모터의 회전속도, 생산속도등의 데이터를 추가로 얻어야만 했다. 또는, 회전 이론에서는 keyphasor event를 이용하라고 한다. 미리 준비되어 있다면 이들을 이용하는 것이 바람직하다. 만일 모든 조건이 다 구성되어 있다면 해당 로직으로 속도 확인을 하는 것도 좋을 듯하다.
댓글