Python으로 FFT처리하여 그래프를 그려보자

[사진사랑]

Matplotlib과 그 종속 항목 (Dependencies)은 Windows, macOS, Linux에 대해 wheel 패키지의 형태로 배포됩니다.
아래의 명령어로 공식 배포판을 설치하세요.

python -m pip install -U pip

python -m pip install -U matplotlib

이렇게 하면 matlab의 라이브러리가 올라간다.

그런다음 

다음과 같이 프로그램을 한다.

60Hz에 진폭 1인 시그널과 120Hz에 진폭 2인 시그널이 합해진 신호를 분석하기 위한 신호의 예제 수식을 아래와 같이 정하자.

y=1⋅cos(2π⋅60⋅t+π/2)+2cos(2π⋅120⋅t)

진폭 0.6과 1.0인 두 cosine 신호의 합성으로 만들어져 있음을 알수 있다. 우선 1초동안 2kHz로  샘플링 하여 결과를 보기위한 소스코드는 다음과 같다.

import matplotlib.pyplot as plt

import numpy as np

import math Fs = 2000 # Sampling frequency

T = 1/Fs # Sample interval time

te= 0.5 # End of time

t = np.arange(0, te, T) # Time vector # Sum of a 50 Hz sinusoid and a 120 Hz sinusoid

noise = np.random.normal(0,0.05,len(t))

x = 1*np.cos(2*np.pi*60*t+np.pi/2) + 2*np.cos(2*np.pi*120*t)

y = x # + noise # Sinusoids plus noise

먼저 라이브러리를 import하고, t시간 배열을 0초부터 te초까지 arange()함수를 사용하여 만들어 준다. 예제에서는 2kHz로 샘플링하여 0.5초 까지 데이터를 만들어 준다. 샘플링을 많이해주면 FFT 후에 더 깨끗한 그래프 결과를 볼 수 있지만 의도적은 적은 샘플링을 해주어서 그 의미를 추측해보자. 만들어진 데이터에 노이즈를 적당히 실어주기 위해서는 예제에서 리마크한 #을 제거해 주어 y=x+noise와 같이 해주면 실제 신호같은 효과를 줄 수 도 있다. 만들어진 신호는 아래와 같이 간략히 그래프로 그려주자.

plt.figure(num=1,dpi=100,facecolor='white')

plt.plot(t,y,'r')

plt.xlim(0, 0.05)

plt.xlabel('time($sec$)')

plt.ylabel('y') plt.savefig("./test_figure1.png",dpi=300)

를 프로그램하면 아래와 같은 그래프가 그려진다.

기본적인 Signal이 생성되었으니 FFT를 처리하면 된다. 분석할 신호의 개수를 알아내기 위해서 len(y)하여 FFT의 샘플개수를 결정한다. 이때 matlab등에서는 일반적으로 nextpow2()등의 함수를 사용하여 취한 샘플개수보다 많은 2n개의 샘플개수를 NFFT로 결정하지만 파이썬에서는 그렇게 해주지 않아도 알아서 내부에서 처리한다. 다음에 f0=np.arange(NFFT)*Fs/NFFT 식으로 더블사이드 주파수 배열을 구한다. 나머지 반대쪽은 folding하여 싱글사이드와 동일하므로 샘플링 주파수의 절반은 f0=f0[range(math.trunc(NFFT/2))] 방법으로 잘라낸다. 이때 trunc()는 정수화하여 에러를 피하기위해 사용한다.

이제 Y=np.fft.fft(y)/NFFT 를 사용하여 FFT를 수행하고 복소수 값을 Y에 리턴하여 대입해 준다. 마찬가지고 싱글사이드 주파수 영역의 값을 구한다. 해석하여 구하여진 주파수별 복소수값은 2*abs()함수를 사용하여 amplitude를 구할 수 있으며, 이때 2를 곱하는 것은 더블사이드에 동일한 크기의 amplitude값이 있으므로 2를 곱하여 표현해 주어야 원래 신호의 진폭 값에 거의 근사해 진다. (전자공학에서는 신호의 크기에 대하여 상대적 비교만 하기때문에 좀더 간편하게 FFT 구현에 접근하지만 기계공학에서는 분석치의 물리적인 값의 표현에 좀더 민감하다. 그래서 전자공학쪽 FFT amplitude 예제에는 가끔 2의 배수가 표현되어 있지 않다.)  주파수 스팩트럼은 신호의 주파수 대역별 에너지 개념이라는 것을 생각하여 이해하자. phase역시 간략히 np.angle()함수를 사용하여 구하면 된다.

# Calculate FFT ....................

n=len(y) # Length of signal

NFFT=n # ?? NFFT=2^nextpow2(length(y)) ??

k=np.arange(NFFT) f0=k*Fs/NFFT # double sides frequency range

f0=f0[range(math.trunc(NFFT/2))] # single sied frequency range

Y=np.fft.fft(y)/NFFT # fft computing and normaliation

Y=Y[range(math.trunc(NFFT/2))] # single sied frequency range

amplitude_Hz = 2*abs(Y)

phase_ang = np.angle(Y)*180/np.pi

로 코딩된다.

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

#St=0.0005

#Fs = 1/St                    # Sampling frequency
Fs = 2000                    # Sampling frequency
T = 1/Fs                     # Sample interval time
te= 0.5                     # End of time
t = np.arange(0, te, T)   # Time vector

# Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
noise = np.random.normal(0,0.05,len(t))
x = 1*np.cos(2*np.pi*60*t+np.pi/2) + 2*np.cos(2*np.pi*120*t)
y = x # + noise     # Sinusoids plus noise

# if energy was collapsed by forced code like belows...
#BB=range(math.trunc(len(t)/2), len(t))
#y[BB]=0

plt.figure(num=1,dpi=100,facecolor='white')
plt.plot(t,y,'r')
plt.xlim( 0, 0.05)
plt.xlabel('time($sec$)')
plt.ylabel('y')

plt.savefig("./test_figure1.png",dpi=300)


# Calculate FFT ....................
n=len(y)        # Length of signal
NFFT=n      # ?? NFFT=2^nextpow2(length(y))  ??
k=np.arange(NFFT)
f0=k*Fs/NFFT    # double sides frequency range
f0=f0[range(math.trunc(NFFT/2))]        # single sied frequency range

Y=np.fft.fft(y)/NFFT        # fft computing and normaliation
Y=Y[range(math.trunc(NFFT/2))]          # single sied frequency range
amplitude_Hz = 2*abs(Y)
phase_ang = np.angle(Y)*180/np.pi

# figure 1 ..................................
plt.figure(num=2,dpi=100,facecolor='white')
plt.subplots_adjust(hspace = 0.6, wspace = 0.3)
plt.subplot(3,1,1)

plt.plot(t,y,'r')
plt.title('Signal FFT analysis')
plt.xlabel('time($sec$)')
plt.ylabel('y')
#plt.xlim( 0, 0.1)

# Amplitude ....
#plt.figure(num=2,dpi=100,facecolor='white')
plt.subplot(3,1,2)

# Plot single-sided amplitude spectrum.

plt.plot(f0,amplitude_Hz,'r')   #  2* ???
plt.xticks(np.arange(0,500,20))
plt.xlim( 0, 200)
plt.ylim( 0, 1.2)
#plt.title('Single-Sided Amplitude Spectrum of y(t)')
plt.xlabel('frequency($Hz$)')
plt.ylabel('amplitude')
plt.grid()

# Phase ....
#plt.figure(num=2,dpi=100,facecolor='white')
plt.subplot(3,1,3)
plt.plot(f0,phase_ang,'r')   #  2* ???
plt.xlim( 0, 200)
plt.ylim( -180, 180)
#plt.title('Single-Sided Phase Spectrum of y(t)')
plt.xlabel('frequency($Hz$)')
plt.ylabel('phase($deg.$)')
plt.xticks([0, 60, 120, 200])
plt.yticks([-180, -90, 0, 90, 180])
plt.grid()

plt.savefig("./test_figure2.png",dpi=300)
plt.show()

실행을 하면 아래와 같이 그려진다.

60Hz는 1, 120Hz는 2라는 크기를 가진 COS함수를 합성된 파형이 위 파형이고 위 파형을 FFT처리하면 중간에 있는 파형이다.

세번째 그래프는 위상으로 추적한 그래프이다.

단 이프로그램은 2개의 COS신호를 수식에 의해서 합성하고 합성된 것을 FFT변환한것이니 시뮬레이션 하는데 유용할 것이다.