파이선 Matplotlib 라이브러리 사용 그래프 작성 I

※파이선 코딩 초보자를 위한 텐서플로우∙OpenCV 머신 러닝 2차 개정판 발행

http://blog.daum.net/ejleep1/1175

파이선 코딩 초보자를 위한 텐서플로우∙OpenCV 머신 러닝 2차 개정판 (하이퍼링크) 목차 pdf 파일

 

아나콘다 편집기에서 파이선 코드 작업에서 그래프 작성을 위한 matplotlib 라이브러리 사용법을 사례를 통해 알아보기로 하자. 아나콘다를 설치한 직 후에 아나콘다 네비게이터를 실행하게되면 아래와 같이 base(root) 환경이 나타나게 된다.

base(root)는 Environments에서 특정한 가상환경을 설정하기 이전에 기본적으로 제공되는 파이선 편집 환경으로서 아주 기본적인 파이선 라이브러리들은 이미 설치가 되어 있다. 예를 들자면 numpy 또는 matplotlib 들일 것이며 다음과 같이 Environments에 들어가서 해당 라이브러리가 설치되어 있는지 쉽게 확인이 가능하다. 따라서 별도의 설치가 불필요하다.

시간에 따라 변하는 함수 그래프 작성

matplotlib를 활용한 사례로서 다음과 같은 함수를 대상으로 그래프를 작도해 보자. 여기서 t는 시간 변수를 뜻하며 그 범위는 0.5초로 둔다.

f(t) = 0.5*cos(2π*60*t+π/2) + cos(2π*120*t)

상기 함수의 그래프 작성을 위해서 코드 초반에 라비브러리를 불러 들이자.

 

import matplotlib.pyplot as plt

import numpy as np

import math

 

시간 t에 따라 변화하는 함수 f(t)의 그래프 작성을 위해서 샘플링 주파수 Fs를 1 kHz 로 두자. 즉 시간 변수 t 매 1초 당 1000개의 점을 찍어 연결하여 그래프를 작성하기로 한다.

 

Fs = 1000 # 샘플링 주파수

dt = 1/Fs # 샘플링 시간 구간

te= 0.5 # 그래프 작성을 위한 가로 축 시간 구간 법위

t = np.arange(0, te, dt) # Time vector

 

한 점당 출력 시간은 0.001초이며 0.0 ~ 0.5초에 이르는 구간에서 함수를 작도한다.

하나의 함수 그래프를 작도하기 위해서는 다음과 같은 코드작업이 기본적으로 필요하게 된다.

⓵ 그래프 그룹의 번호(num), 해상도(dpi)와 배경 색상(facecolor)을 지정할 필요가 있다.

⓶ 변수(t), 함수 y=f(t), 선의 색상지정(‘r’: 빨간색)

⓷ 변수 구간의 최대상한 범위 지정(0, 0.05)

⓸ 가로축 즉 xlabel 변수명(‘time(sec))

⓹ 세로축 즉 ylabel 변수명(‘y’)

⓺ 출력된 그래프를 현재의 폴더에 저장 (./그래프이미지파일명, dpi)

 

plt.figure(num=1,dpi=300, facecolor='red') #그래프 그룹의 번호(num),해상도(dpi)

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

plt.xlim( 0, 0.05)

plt.xlabel('time(sec)')

plt.ylabel('y')

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

 

출력된 그래프는 다음과 같다.

 

시간에 따라 변하는 함수의 스펙트럼 그래프 작성 사례

시간에 따라 변하는 함수 f(t) = 0.5*cos(2π*60*t+π/2) + cos(2π*120*t) 의 특성을 살펴보자. 앞 항의 0.5는 주파수가 60Hz이며 위상 값이 (π/2)인 cosine 함수의 진폭을 나타낸다. 아울러 뒷항은 진폭이 1.0이며 주파수가 120Hz를 의미한다.

즉 앞서 작도한 함수f(t)는 푸리에 급수(Fourier Seiries)로 표현하면 60Hz 와 120Hz 주파수 항들을 가지므로 해 당 주파수의 진폭을 사용하여 표현된다. 이와 같이 고유 주파수와 이들의 배수를 사용하여

함수 f(t)를 해당 주파에서 진폭을 표현한 것이 스펙트럼 그래프에 해당한다. 물론 함수 자체가 예에서와 같이 cosine 급수형태처럼 표현되는 것이 아니라 지정된 점의 수만큼 숫자 값으로 주어진 경우 이들을 푸리에 급수로 변환하려면 FFT 알고리듬이 필요하며 파이선 numpy 라이브러리가 제공하는 FFT루틴(fft.fft(∙))을 사용하기로 한다. 루틴명이 fft.fft(∙) 로소 fft 가 두 번 중복되는 이유는 fft라는 변수명이 사용자에 의해 흔하게 사용될 수도 있기 때문에 이로 인한 에러를 피할 수 있도록 의도적으로 작명한 것으로 보인다. fft 알고리듬은 다소 복잡하므로 사용 결과에만 주목하기로 하자.

⓵ len(y)를 실행하면 데이터 점들의 수 n 값이 얻어진다.

⓶ 이 n 값을 FFT 계산을 위한 NFFT 값 설정에 사용하자. 일반적으로 2N 값이 사용되지만 반드시 2N 이 아니라도 파이선 내부에서 알아서 맞춰주게 된다. 가급적이면 맞춰주는 것이 좋을 듯하다.

⓷ fft.fft(∙)루틴을 사용하여 주파수를 계산한다.

⓸ 진폭 amplitude_Hz를 계산하자.

⓹ 아울러 위상도 계산한다.

 

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

n=len(y) # Length of signal

NFFT=n

k=np.arange(NFFT)

f0=k*Fs/NFFT # double sides frequency range

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

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

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

amplitude_Hz = 2*abs(Y)

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

 

FFT 계산이 완료되었으면 묶음 2(num=2)에서 3개의 그래프를 subplot 명령을 사용하여 311, 312, 313 형식으로 순차적으로 그래프를 작도하자. 311이란 3개의 rows 1 column 형식 중에서 첫 번째로 시그널 자체를 출력하고 312에서 두 번째로 진폭을, 313에서 세 번째 위상을 출력하자.

subplots_adjust 명령에서 hspace 와 wspace는 패딩의 높이와 폭 비율을 나타낸다.

 

# figure 1 ..................................

plt.figure(num=2,dpi=100,facecolor='green')

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')

 

두 번째 그림에서 주파수 스펙트럼 그림인 312를 작도해 보자. 앞 그래프 fig1의 311 그래프 가로축 상의 틱 간격이 plt.xticks(np.arange(0,500,20)) 명령을 사용하여 20으로 바뀐 점이다. Default 처리하면 간격 틱을 볼 수 없을 정도이므로 적절한 값을 사용하자. xlim 및 ylim 명령을 사용하여 그래프의 가로축 및 세로축 구간 범위를 설정한다. 마지막으로 fig1과는 달리 plt.grid(∙) 명령을 사용하여 그래프 상에 격막을 형성한다.

 

# Plot single-sided amplitude spectrum.

plt.subplot(3,1,2)

plt.plot(f0, amplitude_Hz,'blue')

plt.xticks(np.arange(0,500,20))

plt.xlim( 0, 200)

plt.ylim( 0, 1.2)

plt.xlabel('frequency($Hz$)')

plt.ylabel('amplitude')

plt.grid()

 

다음의 311,312,313 출력 결과를 살펴보자.

 

특히 fft 루틴에 의해 계산된 주파수 계산 결과를 보고자 한다면 아래와 같이 shell에서 amplitude_Hz 를 출력하도로 하자.  주파수 60HZ 와 120Hz에서 각각 주어진 진폭 0.6 과 1.0 을 확인할 수 있다.

 

 

#fft_01.py

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


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

# Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
noise = np.random.normal(0,0.05,len(t))
x = 0.6*np.cos(2*np.pi*60*t+np.pi/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=300,facecolor='white')
plt.plot(t,y,'r')
plt.xlim( 0, 0.05)
plt.xlabel('time($sec$)')
#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)
m=k
f0=k*Fs/NFFT    # double sides frequency range
l=f0
f0=f0[range(math.trunc(NFFT/2))] # single sides frequency range
Y=np.fft.fft(y)/NFFT             # fft computing and normaliation
Y=Y[range(math.trunc(NFFT/2))]   # single sides frequency range
amplitude_Hz = 2*abs(Y)
phase_ang = np.angle(Y)*180/np.pi

# figure 1 ..................................
plt.figure(num=2,dpi=100,facecolor='green')
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')

# Plot single-sided amplitude spectrum.
plt.subplot(3,1,2)
plt.plot(f0,amplitude_Hz,'blue')   #  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,'black')   #  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()