파이선 Matplotlib 오실로스코프 듀티 파형 FFT 애니메이션 예제 V

 다음의 예에서 처럼 정현파 기본 주파수의 정수배에 해당하는 주파수 파형의 선형 합성으로 이루어지는 사각형 파형을 포함하는 듀티 파형을 FFT 처리한 스펙트럼을 애니메이션 해 보기로 한다.

애니메이션을 위한 코드는 앞서 작성했던 사인파형 FFT 에니메이션 코드와 거의 구조가 같으며 차이점은 scikit 라이브러리의 signal.square 명령을 사용하여 데이터를 구성하면 된다. 아래 왼쪽 그림이 time domain에서의 노이즈가 섞여 있는 듀티파형이며 오른쪽이 주파수 스펙트럼으로서 30Hz가 main frequency임을 알 수 있다.

실제 코드에서 time domain 에서의 파형을 보려면 line.set_data(z, y) 명령을 살리고 스펙트럼을 보려면 line.set_data(z, amplitude_Hz) 명령을 살려 쓰면 된다. 아울러 반드시 아나콘다 프롬프트 창에서 첨부된 코드를 command line 방식으로 실행시켜야 한다.

 

#spectrum_analyzer

from scipy import signal
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation
from matplotlib.animation import FuncAnimation, PillowWriter
import math

# initialization function: plot the background of each frame
def init():
    line.set_data([], [])
    return line,

# animation function.  This is called sequentially
def animate(i):
    Fs = 1000
    t= 640
    x = np.linspace(0, t, Fs)
    z = x
    noise = 2.0 * np.random.normal(0,0.05,len(x))
    y = signal.square(2 * np.pi * 0.05 * x) + noise
    #y = np.sin(2 * np.pi * 60 * (x - 0.01 * i))  + np.sin(2 * np.pi * 120 * (x - 0.01 * i))  + 10.0 * noise
    #+ 0.5 * np.sin(20 * np.pi * 120 * (x - 0.01 * i))
# 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))]
    z = x[range(math.trunc(NFFT/2))]
    Y=Y[range(math.trunc(NFFT/2))]          # single sied frequency range
    amplitude_Hz = 2*abs(Y)
    phase_ang = np.angle(Y)*180/np.pi
    #line.set_data(z, y)
    line.set_data(z, amplitude_Hz)
    return line,

fig = plt.figure(num=1, dpi=100, facecolor='white')
ax = plt.axes(xlim=(0, 320), ylim=(-2, 2))
line, = ax.plot([], [], lw=2,color='red')
plt.xlabel('Hz')
plt.ylabel('y')
# call the animator.  blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=200, interval=20, blit=True)
# save the animation video file 
#anim.save('dutyspectrum.gif',writer='matplotlib.animation.PillowWriter') 
plt.show()