振り子を作ってWebカメラを使って振り子の位置を検出して，グラフにしてみます．
まぁ今日のテストがイミフで，水曜日の第二外国語の勉強したくないだけですけど．
1.振り子
こんなの
文句はなしで行きましょう（笑)．
暇つぶしクオリティなので重りはキャップを2つくっつけただけです．
糸も細くて変形しにくい針金などではなく，裁縫道具に入ってたやつです．
2.検出プログラム
こんな感じ．output.csvに時間，x座標，y座標が保存されます．検出アルゴリズムとしては，HSV色空間にWebカメラの画像を変換し，緑色の成分のみをOpenCVのinRange関数を使って取り出しています．
なので以下のプログラムは「緑色の物体」の座標を認識しているので緑色のものがないところでやるようにしてください．
座標は緑色が一番たくさんあった行と列の要素番号をx,yの座標としてます．
まぁそのへんはプログラムを参照してください．
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # -*- coding: utf-8 -*- | |
| import numpy as np | |
| import cv2 | |
| import time | |
| cap = cv2.VideoCapture(0) | |
| f = open("output.csv","w") | |
| f.write("t,x,y\n") | |
| start = time.time() | |
| while(True): | |
| # フレームをキャプチャする | |
| ret, frame = cap.read() | |
| # hsv色空間に変換 | |
| hsv = cv2.cvtColor(frame,cv2.COLOR_BGR2HSV) | |
| # 0 <= h <= 179 (色相) OpenCVではmax=179なのでR:0(180),G:60,B:120となる | |
| # 0 <= s <= 255 (彩度) 黒や白の値が抽出されるときはこの閾値を大きくする | |
| # 0 <= v <= 255 (明度) これが大きいと明るく，小さいと暗い | |
| # ここでは緑色を抽出するので60±20を閾値とした | |
| low_color = np.array([40, 75, 75]) | |
| upper_color = np.array([80, 255, 255]) | |
| # 色を抽出する | |
| ex_img = cv2.inRange(hsv,low_color,upper_color) | |
| # 縦，横方向の和の中で最も大きいインデックスを抽出 | |
| x = np.argmax(np.sum(ex_img,axis=0)) | |
| y = np.argmax(np.sum(ex_img,axis=1)) | |
| # 時間と座標の保存 | |
| t = time.time() - start | |
| f.write(str(t)+","+str(x)+","+str(y)+"\n") | |
| print(str(t)+","+str(x)+","+str(y)+"\n") | |
| # 抽出した座標に丸を描く | |
| cv2.circle(frame,(x,y),10,(0,255,0)) | |
| # 画面に表示する | |
| cv2.imshow('frame',frame) | |
| if t > 120: | |
| break | |
| key = cv2.waitKey(1) & 0xFF | |
| # qが押された場合は終了する | |
| if key == ord('q'): | |
| break | |
| # 後始末 | |
| cap.release() | |
| cv2.destroyAllWindows() | |
| f.close() |
3.プロット
データがoutput.csvに入ったら次はmatplotlibで可視化してみます．可視化プログラムはこんな感じです．
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #coding:utf-8 | |
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| import pandas as pd | |
| # 最初のデータはWebカメラの調子が悪いため無視 | |
| data = pd.read_csv("output.csv") | |
| data = data[50:] | |
| # 計算しやすいように取り出す | |
| t = data['t'] | |
| x = data['x'] | |
| y = data['y'] | |
| # x,yは平均(振動の中心)を0に持ってこさせる | |
| x = x - np.average(x) | |
| y = y - np.average(y) | |
| # 時間はt[0]が0になるようにしとく | |
| # (pandasを使うとt[0]という表現ができないのでこうしました) | |
| t = t - np.min(t) | |
| # 凡例とともに表示 | |
| plt.plot(t,x,label="x") | |
| plt.plot(t,y,label="y") | |
| plt.legend() | |
| plt.show() |
今日は以上です．
振り子の長さは大体37cmだったので，T=2*pi*sqrt(l/g)で，周波数はおおよそ0.8[Hz]になっているはずで，本当はここでフーリエ変換で周波数を求めようと思ったのですが……これ一定間隔でサンプリングしてないのです．
numpyやscipyのフーリエ変換の関数はたぶんサンプリング間隔が一定という条件でやってる気がする……時間のデータを与えられないし．
複素数をPythonで扱ったりとかはできず，現状実装できる能力がないので近いうちにフーリエ変換を実装しなければ……．
