手書きモデルで固有値解析

以前、試した Google の機械学習サービス Google Cloud Vision API と iPhone の Pythonista3 を使って手書きモデルから固有値解析を実行するアプリを作成してみました。

動作の雰囲気は以下です。

Python ソースは以下のとおり。

# -*- coding: utf-8 -*-
import appex
import numpy as np
from matplotlib import pyplot as plt
from requests_toolbelt.multipart.encoder import MultipartEncoder

import base64
import json
from requests import Request, Session
from io import BytesIO
from PIL import Image

#PILで開いた画像をbase64形式に変換する
def pil_image_to_base64(pil_image):
    buffered = BytesIO()
    pil_image.save(buffered, format="PNG")
    str_encode_file = base64.b64encode(buffered.getvalue()).decode("utf-8")
    return str_encode_file


def pil_image_to_base642(fp):
    str_encode_file = MultipartEncoder(fields={'encoded_image': (fp.name, fp)})
    return str_encode_file


#PILで開いた画像をCloud Vision APIに投げる
def recognize_image(pil_image):
        str_encode_file = pil_image_to_base64(pil_image)
        str_url = "https://vision.googleapis.com/v1/images:annotate?key="
        str_api_key = "取得したAPIキー"
        str_headers = {'Content-Type': 'application/json'}
        str_json_data = {
            'requests': [
                {
                    'image': {
                        'content': str_encode_file
                    },
                    'features': [
                        {
                            #'type': "TEXT_DETECTION",
                            'type': "DOCUMENT_TEXT_DETECTION",
                            'maxResults': 10
                        }
                    ]
                }
            ]
        }

        obj_session = Session()
        obj_request = Request("POST",
                              str_url + str_api_key,
                              data=json.dumps(str_json_data),
                              headers=str_headers
                              )
        obj_prepped = obj_session.prepare_request(obj_request)
        obj_response = obj_session.send(obj_prepped,
                                        verify=True,
                                        timeout=60
                                        )

        if obj_response.status_code == 200:
            with open('data.json', 'w') as outfile:
                json.dump(obj_response.json(), outfile)
                text = get_fullTextAnnotation(obj_response.text)
            return text

        else:
            return "error"
#返ってきたjsonデータの"fullTextAnnotation"部分のテキストを抽出する
def get_fullTextAnnotation(json_data):
    text_dict = json.loads(json_data)
    try:
        text = text_dict["responses"][0]["fullTextAnnotation"]["text"]
        return text
    except:
        print(None)
        return None

# 手書き画像を共有Extensionで取得する
pil_image = appex.get_image()

# Google Cloud Vision API に画像を投げて取得したデータから重量とバネ剛性を取得する
res = recognize_image(pil_image)
sp = res.splitlines()

m = list()
k = list()

print('◼️Google Cloud Vision API の認識結果')
for v in sp:
    v = v.strip()
    v = v.replace(' ', '')
    v = v.replace(',', '.')
    d = v.split('=')
    print(v, d)
    if d[0] == 'm' or d[0] == 'M':
        m.append(float(d[1]))
    if d[0] == 'k' or d[0] == 'K':
        k.append(float(d[1]))

print("◼️バネマス情報")
for v in m:
    print('m =',v)
for v in k:
    print('k =',v)

M1 = len(m)
M2 = M1 + 1
M = np.zeros((M2, M2))
K = np.zeros((M2, M2))

# 質量マトリクスと剛性マトリクスを作成
for i in range(len(m)):
    M[i][i] = m[i] / 9.80665

for i in range(len(k)):
    i1 = i
    i2 = i + 1
    K[i1][i1] = K[i1][i1] + k[i]
    K[i2][i2] = K[i2][i2] + k[i]
    K[i1][i2] = K[i1][i2] - k[i]
    K[i2][i1] = K[i2][i1] - k[i]

# 最後のばねの端部が固定点として行と列を削除して拘束条件とする
M = np.delete(M, M1, 0)
M = np.delete(M, M1, 1)
K = np.delete(K, M1, 0)
K = np.delete(K, M1, 1)

# 質量マトリクスの逆行列を計算
M_inv = np.linalg.inv(M)

# 固有値と固有ベクトルを計算
omega, v = np.linalg.eig(np.dot(M_inv, K))

# 固有値の順番を昇順にソート
freq = list()
omega_sort = np.sort(omega)
for o in omega_sort:
    freq.append(np.sqrt(o)/(2.0*np.pi))

# 固有値のソート時のインデックスを取得
# ⇒固有ベクトルと対応させるため
sort_index = np.argsort(omega)

# 固有値に対応する固有ベクトルをソート
v_sort = []
for i in range(len(sort_index)):
    v_sort.append(v.T[sort_index[i]])
v_sort = np.array(v_sort)

print('◼️固有値解析結果')
# 結果をコンソールに表示
for i in range(len(omega_sort)):
    s = '{0:d}次 {1:10.3f} Hz'.format(i+1,freq[i])
    print(s)
    for j in range(0,v_sort.shape[1]):
        s = '{:10.3f}'.format(v_sort[i][j])
        print(s)

# グラフ化のために自由度軸を作成
dof = np.linspace(len(sort_index), 0, len(sort_index)+1)

# ここからグラフ描画
# フォントの種類とサイズを設定する。
plt.rcParams['font.size'] = 9
plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。
fig = plt.figure(figsize=(3,6))
ax1 = fig.add_subplot(111)
ax1.yaxis.set_ticks_position('both')
ax1.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。
ax1.set_xlabel('Eigen Vector')
ax1.set_ylabel('Degree of freedom')

# データの範囲と刻み目盛を明示する。
ax1.set_xticks(np.arange(-2, 2, 0.5))
ax1.set_yticks(np.arange(0, dof[0]+1, 1))
ax1.set_xlim(-1, 1)
ax1.set_ylim(0, dof[0])

# データプロット 固有ベクトルの形を見る
for i in range(len(sort_index)):
    eigen_vector = np.concatenate([v_sort[i], [0]])
    s = 'Mode' + str(i+1) + ' {:10.3f}Hz'.format(freq[i])
    ax1.plot(eigen_vector, dof, lw=1, marker='o', label=s)

fig.tight_layout()
plt.legend(bbox_to_anchor=(0, 0.1), loc='lower left')

# グラフを表示する。
plt.show()
plt.close()

固有値解析のソースはWATLABさんのこちらを参考にさせていただきました。

iPhoneのExcelデータをPythonista3でフーリエ変換

Pythonista3はiOSで動作するPythonの開発環境です。
numpyが標準でインストールされていましたが、openpyxlがインストールされていたのは少し驚きました。

ということで、iPhoneのExcelデータを読み込んでfftスペクトルを計算するアプリを作成してみました。

動作の雰囲気はこんな感じです。

ソースファイルは以下です。

# -*- coding: utf-8 -*-
 import appex
 import numpy as np
 import matplotlib.pyplot as plt
 import openpyxl as px
 import ui
 

 def fft(sender):
     dt = float(sender.superview['textfield1'].text)
     print('時間刻み = ',dt)
     t = list()
     f = list()
     
     count = 0
     for i in range(len(d)):
         t.append(dt*float(count))
         count += 1
     f = np.linspace(0, 1.0/dt, count) 

     # 高速フーリエ変換
     F = np.fft.fft(d)/float(count)
 
     # 振幅スペクトルを計算
     Amp = np.abs(F)
     # 位相角を計算
     Phs = np.degrees(np.angle(F))
 
     fmax = 1.0/dt/2.0 #後半の共役複素数部は除外
 
     # グラフ表示
     plt.figure(figsize=(10,15))
     plt.subplot(311)
     plt.plot(t, d)
     plt.xlabel("time")
     plt.ylabel("signal")
 
     plt.subplot(312)
     plt.plot(f, Amp)
     plt.xlim(0, fmax)
     plt.xlabel("frequency")
     plt.ylabel("amplitude")
 
     plt.subplot(313)
     plt.scatter(f, Phs, s=5, c="blue", edgecolors="blue")
     plt.xlim(0, fmax)
     plt.ylim(-180, 180)
     plt.xlabel("frequency")
     plt.ylabel("phase")
     plt.show()
     v.close()
 
 
 d = list()
 myf = appex.get_file_path()
 wb = px.load_workbook(myf, data_only = True)
 ws = wb.worksheets[0]
 for cell in ws['A']:
     d.append(cell.value)
 
 v = ui.load_view()
 v.present('fullscreen')