オンライン周波数解析アプリを作ってみた


背景

今回、周波数解析アプリを作ってみた経緯として、現在私が所属している大学の研究室では、脈波や心電図について研究している班が複数あるのですが、それらを解析するためのソフト(WaveLab有料ソフト)が入っているPCが1台しかないという現状です。
この現状により一つの班がソフトを使用している時に別の班は解析できないので、自分で作ってみようと思いました。
また、オンラインで周波数解析ができるようになると、ネットに繋がってさえいれば、わざわざソフトをインストールせずにどこでも解析ができるようになり、とても便利になると考えています。

使用した言語について

このWebアプリケーションを作成するのにPythonを使用しました。理由として、Pythonで有名なライブラリのNumpyとScipyを使って作成しようと思ったからです。この2つのライブラリは、周波数解析をするための機能がたくさん入っているので、これらを使って作成していきます。
フレームワークでDjangoを使用した理由はPythonのFWと言えば、Djangoなのかなと思ったので使用しました。

周波数解析モジュールの作成

脈波や心電図を解析するためのモジュールを作成します。
モジュール名はanalysis.pyにしています。

analysis.py
import io
import pathlib
import numpy as np
from scipy import signal

import matplotlib.pyplot as plt

from .models import Pulse_Rate

def setPlt(pk):
    # 対象データを取得
    pulse_rate = Pulse_Rate.objects.get(pk=pk)
    path = pathlib.Path(pulse_rate.data.url)

    # データを読み込む
    f = open(path.resolve())
    lines = f.read().split()

    # -----変数の準備-----
    N = 4096
    lines_length = int(len(lines) / 2)
    if lines_length < N:
        N = 2048

    dt = float(lines[2]) - float(lines[0])
    pulse = []
    for num in range(N - 1):
        pulse.append(float(lines[num * 2 + 1]))
    # -------------------

    # -----サンプリング周波数計算-----
    t = np.arange(0, N * dt, dt)  # time
    freq = np.linspace(0, 1.0 / dt, N)  # frequency step
    # -----------------------------

    # -----波形の生成-----
    y = 0
    for pl in pulse:
        y += np.sin(2 * np.pi * pl * t)
    # -------------------

    # -----フーリエ変換-----
    yf = np.fft.fft(y)  # 高速フーリエ変換
    # --------------------

    # パワースペクトル算出
    yf_abs = np.abs(yf)

    # -----グラフの生成-----
    plt.figure()
    plt.subplot(211)
    plt.plot(t, y)
    plt.xlabel("time")
    plt.ylabel("amplitude")
    plt.grid()

    plt.subplot(212)
    plt.plot(freq, yf_abs)
    plt.xlim(0, 10)
    plt.xlabel("frequency")
    plt.ylabel("amplitude")
    plt.grid()
    plt.tight_layout()

こんな感じで脈波や心電図のグラフとパワースペクトルのグラフを生成します。
データはユーザーがアップロードしたテキストファイルから読み出します。

Webページにグラフを表示する

上のモジュールではplotでグラフを作成しているので、これをsvg形式でWebページに表示したいと思います。

analysis.py
def pltToSvg():
    buf = io.BytesIO()
    plt.savefig(buf, format='svg', bbox_inches='tight')
    s = buf.getvalue()
    buf.close()
    return s

これでplotをsvgに変換してWebページに表示することができます。

views.pyの設定

上記のモジュールを使用してviews.pyを書きます。

views.py
def get_svg(request, pk):
    setPlt(pk)       # create the plot
    svg = pltToSvg() # convert plot to SVG
    response = HttpResponse(svg, content_type='image/svg+xml')
    return response

これでHttpResponse(svg形式)でレスポンスを返します。

urls.pyの設定

URLの設定をします。

urls.py
urlpatterns = [
    ...
    path('pulse_rate-detail/<int:pk>/plot/', views.get_svg, name="plot"),
]

templatesの設定

今回グラフを表示するページは詳細画面にしたいと思います。

pulse_rate_detail.html
<img src="{% url 'pulse_rate:plot' object.pk %}" width=600, height=300>

これでグラフを表示することができます。

作成したグラフ

実際に作成できるグラフは以下のようになります。

フィルタ処理をしたグラフの生成

上の図はフィルタ処理をしていないデータをグラフ化しているので、このままでは何もわかりません。
そこでフィルタ処理の機能を追加してグラフを生成しました。(コードについては省きます)

フィルタをかけるとなんとなく脈波っぽい波形がでます。(この脈波はカメラで撮影したRGB値のG成分を利用して出しているので、完璧な波形にはなりません)
フィルタの下限周波数と上限周波数はWebページから自分で調整することができるようにしています。

まとめ

今回作成したアプリケーションは今後研究室で使用していく予定です。
また、私の所属している研究室では、心電図や脈波以外の研究をしている班もあるため、それらの班のためにも何か作成したいと考えています。
実際、私の研究テーマは人の鼻の温度からその人のストレス度を測定するというテーマをしております。

私は自分の所属している研究室を大学内で一番イケてる研究室にしたいと考えているので4年生になってもこのようなアプリケーションを作成したいです!