[Monitoring System] 4. UDP 통신으로 실시간 진동 데이터 FFT / STFT시각화(numpy, tensorflow, pytorch)

결론부터 말하면 stft는 부하가 너무 심해서 실시간으로 시각화하는 것은 무리무리무리이다. 일단 matplotlib 자체에서 단순 시각화를 하기에도 데이터 양이 너무 많기 때문에 일종의 녹화기능을 넣어서 원하는 만큼의 데이터를 fft / stft 처리해서 따로 시각화하기로 결정하였다. 그래도 FFT 정도는 실시간으로 가능할 것 같아서 gpu accelerator를 사용해서 구현하였다.  그러나 CUDA 환경이 아닌 곳에서도 쓸 수 있도록 numpy로도 코드를 추가했다.

이전까지 코드는 요기에 있다. https://dnai-deny.tistory.com/65

[Monitoring System] 3. UDP 통신으로 실시간 진동 데이터 Plotting(feat. matplotlib)

- PyTorch로 구현한 버전

import torch
import matplotlib.pyplot as plt
import matplotlib.animation as animation

class RealTimeAnimation:
    """
    Animate original voltage signal graph and fft result.
    """
    def __init__(self):
        # define plot
        self.fig = plt.figure(figsize=(20, 10))
        self.fig.suptitle('Vibration Data', fontsize=20)
        plt.subplots_adjust(hspace=0.3)

        self.original = self.fig.add_subplot(211)
        self.original.title.set_text('Signal Graph')
        self.original.set(xlabel="samples", ylabel="mV")

        self.fft_graph = self.fig.add_subplot(212)
        self.fft_graph.title.set_text('FFT Graph')
        self.fft_graph.set(xlabel="frequency", ylabel="amplitude")

        # connect animation method
        self.animation = animation.FuncAnimation(self.fig, self.update, interval=50)

    def update():
    ...
        # update fft graph
        Fs = 25600.0
        Ts = 1/Fs
        n = len(graph_data)

        # pytorch fft
        graph_data = torch.Tensor(graph_data)
        graph_data.to("cuda:0")
        fft_data = torch.fft.rfft(graph_data)
        frequency = torch.arange(0.0, Fs/2.0, Fs/n)


        self.fft_graph.clear()
        self.fft_graph.plot(frequency[:len(fft_data)//2], np.abs(fft_data[:len(fft_data)//2]))
        self.fft_graph.title.set_text('FFT Graph')
        self.fft_graph.set(xlabel="frequency", ylabel="amplitude")

__init__ 함수에 fft graph를 그릴 부분을 새로 넣어주고, update 함수에서 fft를 처리하도록 추가하였다. FFT 자체는 신호데이터만 있으면 되기 때문에 rfft(음수 주파수를 날려주는 버전의 fft, 원하지 않으면 fft로 하면 된다)를 적용해 graph data에 대한 fft data를 얻었다.

그리고 해당하는 frequency 계산을 위해서 Fs(Sampling rate)의 절반으로 잡고 그래프 시각화를 위한 x좌표인 frequency를 계산하였다.

Sampling rate는 실제 측정하고 싶은 주파수 범위의 적어도 두 배를 잡아줘야한다. 이에 대한 내용은 나이퀴스트 주파수와 엘리어싱에 대해 찾아보면 더 자세하고 전문적인 이유를 알아볼 수 있다. 

FFT는 좌우가 대칭적으로 나오기 때문에 한쪽(반)만 보면 된다다. 따라서 plotting할 때 fft_data의 절반만 시각화하였다. 절대값 함수를 취해주는 이유는, 계산 결과가 복소수인 경우를 피하기 위해서이다.

짜잔! 실행해보면 잘 나온다. 아래 두 블럭은 각각 numpy버전, tensorflow 버전이다.

from numpy.fft import rfft, rfftfreq
frequency = rfftfreq(n, Ts)[:-1]
fft_data = (rfft(graph_data)/n)[:-1] * 2

import tensorflow as tf

with tf.device("/device:GPU:0"):
    fft_data = tf.signal.rfft(input_tensor=tf.cast(graph_data, tf.float32))
    frequency = tf.range(0.0, tf.divide(Fs,2.0), tf.divide(Fs, tf.cast(n, tf.float32)))

 

2. 녹화 및 일시정지 버튼 추가

class RealTimeAnimation:
    """
    Animate original voltage signal graph and fft result.
    """
    def __init__(self):
        ...

        # define buttons
        self.paused = False
        ax = plt.axes([0.8, 0.025, 0.1, 0.04])
        ax2 = plt.axes([0.68, 0.025, 0.1, 0.04])
        ax3 = plt.axes([0.56, 0.025, 0.1, 0.04])
        ax4 = plt.axes([0.44, 0.025, 0.1, 0.04])

        # plot control button
        self.pause_button = Button(ax, "Pause/Restart")
        self.pause_button.on_clicked(self.toggle_event)

        # start record button
        self.pos = 0
        self.record_button = Button(ax2, "Start Recording")
        self.record_button.on_clicked(self.start_recording)

        # data save button
        self.save_button = Button(ax3, "Save Data")
        self.save_button.on_clicked(self.save_data)

        # show entire data analysis button
        self.show_button = Button(ax4, "Show Analysis")
        self.show_button.on_clicked(show_result)

    def update(self, i) -> None:
        ...

    def toggle_event(self, *args, **kwargs) -> None:
        """
        On click event function of pause button.
        Stop / Restart plotting graphs.
        """

        if self.paused:
            self.animation.event_source.start()
        else:
            self.animation.event_source.stop()

        self.paused = not self.paused

    def start_recording(self, *args, **kwargs) -> None:
        data_lock.acquire()
        self.pos = len(entire_data) - 1
        data_lock.release()
        print("Recording started.")

    def save_data(self, *args, **kwargs) -> None:
        """
        Save entire voltage data in csv.
        During save data, graph plotting will be paused.
        Using lock.
        """

        self.animation.event_source.stop()
        print("Saving...")

        data_lock.acquire()
        f = open("data.csv", "w")
        target = csv.writer(f)
        target.writerow(entire_data[self.pos:])
        data_lock.release()

        print("Saved {0} Data. Restart plotting.".format(len(entire_data)))
        self.animation.event_source.start()

버튼 추가는 말 그대로 버튼을 추가하고 on click 이벤트를 연결하는 부분이다. recording start 버튼을 누르면 현재 data의 인덱스를 저장하고, save data 버튼을 누르면 저장된 인덱스부터 현재 인덱스까지를 csv 형태로 저장해준다. 부하를 막기위해 저장 중에는 plotting을 잠시 멈추도록 하였다.

pause/restart은 plotting을 잠깐 멈추고 재시작할 수 있도록 설정하는 버튼이다.

show result 버튼은 다음 단락에 구현한다.

 

3. 저장한 데이터를 FFT / STFT 분석하기

파일을 새로 만든 후 저장된 데이터를 원본 / FFT / STFT 처리 결과 세 가지로 plotting 한다. 모듈로 호출하도록 할 것이기 때문에 하나의 함수로 등록해두었다.

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import numpy as np
import torch
import csv

def show_result(e) -> None:
    """
    Show original signal graph, fft and stft result of saved data.

    Args:
        e: event
    """

    # read data
    data_file = open("data.csv", "r")
    reader = csv.reader(data_file, delimiter=',', quoting=csv.QUOTE_NONNUMERIC)
    entire_data = []
    for row in reader:
        entire_data = row
        print(entire_data[:10])
        break

    if len(entire_data) <= 1:
        print("No data saved. Please try after save data.")
        return 

    gs = gridspec.GridSpec(3, 1)

    # define graph
    fig = plt.figure(figsize=(20, 12))
    fig.suptitle("Data Analysis Result", fontsize=20)

    original = plt.subplot(gs[0, 0])
    original.set_title("Original Signal")
    original.set(xlabel="Samples", ylabel="mV")

    stft = plt.subplot(gs[1, 0])
    stft.set_title("Short-Term FFT")
    stft.set(xlabel="Time(sec)", ylabel="Frequency")

    fft_graph = plt.subplot(gs[2, 0])
    fft_graph.set_title("FFT")
    fft_graph.set(xlabel="Frequency", ylabel="Amplitude")

    plt.subplots_adjust(hspace=0.5)

여기까지가 데이터를 읽어오고, plot들을 설정하는 부분이다. 

3.1. FFT

    # plot fft
    Fs = 25600.0
    Ts = 1/Fs
    n = len(entire_data)

    # pytorch fft
    graph_data = torch.Tensor(entire_data)
    graph_data.to("cuda:0")
    fft_data = torch.fft.rfft(graph_data) / n

    frequency = torch.arange(0.0, Fs/2.0, Fs/n)
    fft_graph.clear()
    fft_graph.plot(frequency[1:len(fft_data)//2], np.abs(fft_data[1:len(fft_data)//2]))

fft는 위 실시간으로 처리하는 것과 코드가 같다. 자세한 설명은 생략...

3.2. STFT

    # pytorch stft
    graph_data = graph_data.unsqueeze(0)
    stft_ = torch.stft(graph_data, 25600)
    stft_data = torch.sqrt(stft_[:,:,:,0] ** 2 + stft_[:,:,:,1]**2)
    stft.pcolormesh(stft_data[0])

    # scipy stft
    f, t, Zxx = signal.stft(entire_data, fs=Fs, nperseg=len(entire_data)//100)
    stft.pcolormesh(t, f, np.abs(Zxx), shading="gouraud"

STFT도 PyTorch 공식문서를 보고 따라했다. PyTorch의 STFT는 데이터와 sampling rate를 요구한다. 결과를 plot할 수 있는 형태로 바꿔주고 나서 pcolormesh로 시각화하였다.

scipy의 STFT는 데이터, sampling rate, 그리고 window를 요구한다. window를 잘 정하는 방법에 대해서는 아직 잘...모르겠다. 딥러닝 학습할 때 hyperparameter 조정을 하는 기분이라고 해야하나...  window를 넓게 잡으면 주파수 영역에 대한 해상도가 올라가는 대신 시간영역에 대한 해상도가 떨어지고, 좁게 잡으면 시간 영역에 대한 해상도가 개선되는 대신 주파수 영역에 대한 해상도가 떨어진다.

 완성! 전체 코드는 아래 깃헙의 udp_server.py와 stft.py 코드에 자세히 나와있다!  

GitHub - melli0505/monitoring_system: FastAPI/PostgreSQL/MQTT 기반 센서 모니터링 시스템

여담으로 gpu를 사용해서 빠르게 연산하기 위해 torch나 tensorflow의 적용을 구현했지만, 현재로써는 데이터를 gpu 메모리로 복사하는 시간이 존재하기 때문에 체감되는 속도 차이는 별로 없었다. 데이터 양이 많으니 더 그런 듯. 그냥 numpy와 scipy를 사용해도 되겠다.