https://jq1qnv.sakura.ne.jp/body/rotor_model/rotor_model.html
実験の様子の動画(YouTube: https://youtu.be/83qR7vTQLZk)
加速度センサを買ったのは2012年10月と、ずいぶん前で、その頃から構想はあったようなのですが、
しばらくほったらかしにしていて、2015年の春ごろから部品を少しずつそろえ始めてようやく完成しました。
この装置は実用的な何かではありませんが、
ここで取り上げているセンサー類と計測器、分析器を使って、
たとえば自転車等のホイールバランスに使えそうかなと思います。
#!/bin/sh
# TBS 1052B のデータを FFT して表示
OCTDIR="/home/hoge/octave"
csvfile=$(zenity --file-selection --title "TBS 1052B FFT Analysis" --text "send CSV file ")
octcmd="filename='$csvfile' ; cd $OCTDIR ; octfft"
/usr/bin/octave-cli --eval "$octcmd" --persist
Octave スクリプトのパスの設定をしてやるともう少しスマートになると思いますが、
動いているのでまぁいいかと妥協しています。
折角なので findpeaks を使ってピークの表示をしてみました。
表示例は下記。
octave スクリプト octfft.m は以下です。
# USB に保存した TBS1052B の csv データを読み込んで
# FFT 解析、表示させる octave スクリプト
TBS1052Bdata=dlmread( filename ,',');
Timeseri=[TBS1052Bdata(:,4) TBS1052Bdata(:,5)];
# 情報
infofile=fopen( filename, 'r')
TBS1052Binfo=fgetl(infofile);
for i= 1:20
infostr=fgetl(infofile);
TBS1052Binfo=[TBS1052Binfo; infostr];
endfor
fclose(infofile)
sourcech=strsplit(TBS1052Binfo( 7,:),",")
verut =strsplit(TBS1052Binfo( 8,:),",")
verscl =strsplit(TBS1052Binfo( 9,:),",")
veroff =strsplit(TBS1052Binfo(10,:),",")
horut =strsplit(TBS1052Binfo(11,:),",")
horscl =strsplit(TBS1052Binfo(12,:),",")
ptfmt =strsplit(TBS1052Binfo(13,:),",")
y0 =strsplit(TBS1052Binfo(14,:),",")
prbatt =strsplit(TBS1052Binfo(15,:),",")
mdlnum =strsplit(TBS1052Binfo(16,:),",")
selnum =strsplit(TBS1052Binfo(17,:),",")
firmver =strsplit(TBS1052Binfo(18,:),",")
source_ch =sourcech{[2]}
vertical_unit =verut{[2]}
vertical_scale =verscl{[2]}
vertical_offset =veroff{[2]}
horizontal_unit =horut{[2]}
horizontal_scale=horscl{[2]}
pt_fmt =ptfmt{[2]}
yzero =y0{[2]}
probe_atten =prbatt{[2]}
model_num =mdlnum{[2]}
serial_num =selnum{[2]}
firmware_ver =firmver{[2]}
# サンプリング数(Record Length) 2列の1行目
num_of_smpl=TBS1052Bdata(1,2)
# サンプリング周期(Sample Interval) 2列の2行目
smpl_period=TBS1052Bdata(2,2)
# サンプリング長
smpl_len=TBS1052Bdata(1,2)*TBS1052Bdata(2,2)
# FFT
FFTcomplex=fft( Timeseri(:,2))/num_of_smpl ;
# ゼロシフト(直流分を周波数範囲の中央にもってくる)
ZsFFTcomplex=fftshift( FFTcomplex ) ;
# 周波数軸設定
Ax_freq=linspace(-num_of_smpl/2/(smpl_period*num_of_smpl), (num_of_smpl/2-1)/(smpl_period*num_of_smpl), num_of_smpl);
# 周波数軸を含むFFT行列
AxFFTcomplex=[Ax_freq.' ZsFFTcomplex];
# 等価振幅と位相のFFT行列
AxFFTamph=[ AxFFTcomplex((num_of_smpl/2+1):num_of_smpl,1) abs( AxFFTcomplex((num_of_smpl/2+1):num_of_smpl,2))*2 angle( AxFFTcomplex((num_of_smpl/2+1):num_of_smpl,2))*360/(2*pi)];
AxFFTamph(1,:)=[ AxFFTcomplex((num_of_smpl/2+1),1) abs( AxFFTcomplex((num_of_smpl/2+1),2)) angle( AxFFTcomplex((num_of_smpl/2+1),2))*360/(2*pi)];
# 最大値
sorAxFFT=sort( AxFFTamph(:,2), "descend")
indxsf = find( AxFFTamph(:,2) >= sorAxFFT(5) )
max_list = AxFFTamph(indxsf,:)
# ピーク
[pksam idxpk] = findpeaks(AxFFTamph(:,2));
peak_list = [pksam idxpk]
[peakls idxpks] = sort( peak_list(:,1), "descend")
peaksize=size(idxpks)
for i=1:5
if ( i <= peaksize(1) )
peak_list5(i,:) = AxFFTamph(peak_list(idxpks(i),2),:)
else
peak_list5(i,:) = [0.0 0.0 0.0]
endif
endfor
# グラフ
figure('Position',[100,100,560,630])
subplot(4,1,1)
plot(Timeseri(:,1),Timeseri(:,2))
axis([min(Timeseri(:,1)), max(Timeseri(:,1)), -max(abs(Timeseri(:,2)))*1.1, max(abs(Timeseri(:,2)))*1.1])
title(filename)
ylabel(strcat( 'amp ','[',vertical_unit,']'))
xlabel(strcat( 'time ','[',horizontal_unit,']'))
subplot(4,1,2)
plot(AxFFTamph(:,1),AxFFTamph(:,2),AxFFTamph(idxpk,1),AxFFTamph(idxpk,2),'.m');
ylabel( strcat( 'amp ','[',vertical_unit,']'))
hx=xlabel( strcat( 'freq ','[1/',horizontal_unit,']'))
axis([ min(AxFFTamph(:,1)), max(AxFFTamph(:,1))])
ph=subplot(4,1,3)
plot(AxFFTamph(:,1),AxFFTamph(:,3));
set(ph,"ytick",[-180 -90 0 90 180])
ylabel('phase [deg]')
xlabel( strcat( 'freq ','[1/',horizontal_unit,']'))
axis([ min(AxFFTamph(:,1)), max(AxFFTamph(:,1)) ,-180,180])
#bottom_title(strcat( 'number of smple :',sprintf("%d",num_of_smpl)),'sampling period:',sprintf("%e ",smpl_period),"sec")
subplot(4,1,4)
plot(1)
axis([0, 1, 0, 1])
text( 0, 0.9, [model_num " " firmware_ver " S/N " serial_num ;"Record Length " sprintf("%d [samples] %1.3e [sec] ",num_of_smpl,smpl_len) ; "Sample Interval " sprintf("%1.3e [sec]",smpl_period)])
#text( 0, 0.5, [model_num " " firmware_ver ; "S/N " serial_num ;"Record Length " num_of_smpl " Sample Interval " smpl_period])
#xlabel( [model_num " " firmware_ver ; serial_num] )
#text( 0, 0.2, [ "freq Amp Pha"; sprintf("%1.4e %1.4e %1.4e", max_list(1,1), max_list(1,2), max_list(1,3)); sprintf("%1.4e %1.4e %1.4e", max_list(2,1), max_list(2,2), max_list(2,3)) ])
text( 0, 0.1, [ "peak list" ; "freq Amp Pha"; sprintf("%1.4e %1.4e %3.3f", peak_list5(1,1), peak_list5(1,2) , peak_list5(1,3)); \
sprintf("%1.4e %1.4e %3.3f", peak_list5(2,1), peak_list5(2,2) , peak_list5(2,3)); \
sprintf("%1.4e %1.4e %3.3f", peak_list5(3,1), peak_list5(3,2) , peak_list5(3,3)); \
sprintf("%1.4e %1.4e %3.3f", peak_list5(4,1), peak_list5(4,2) , peak_list5(4,3)); \
sprintf("%1.4e %1.4e %3.3f", peak_list5(5,1), peak_list5(5,2) , peak_list5(5,3))])
axis("off")
あとは窓関数の機能を盛り込めば 1ch の FFTアナライザとして一通りの機能をもったものになります。
ここまで来るとカーソルで値を読み取る機能も欲しいですが、
それはちょっとめんどくさそう・・・
表示例
zenity を使った octave 起動スクリプト
#!/bin/sh
# TBS 1052B のデータを FFT して表示
csvfile=$(zenity --file-selection --title "TBS 1052B FFT Analysis" --text "send CSV file ")
/usr/bin/octave --eval 'filename="'$csvfile'"' --persist /home/hoge/fuga/octfft.m
TBS1052B のデータをFFTするスクリプト octfft.m
# TBS1052B の CSV データを FFT してグラフ表示
TBS1052Bdata=dlmread( filename ,',');
Timeseri=[TBS1052Bdata(:,4) TBS1052Bdata(:,5)];
TBS1052Binfo=textread( filename, '%s')
# 情報
sourcech=strsplit(TBS1052Binfo{[22]},",")
verut =strsplit(TBS1052Binfo{[26]},",")
verscl =strsplit(TBS1052Binfo{[30]},",")
veroff =strsplit(TBS1052Binfo{[34]},",")
horut =strsplit(TBS1052Binfo{[38]},",")
horscl =strsplit(TBS1052Binfo{[42]},",")
ptfmt =strsplit(TBS1052Binfo{[46]},",")
y0 =strsplit(TBS1052Binfo{[49]},",")
prbatt =strsplit(TBS1052Binfo{[53]},",")
mdlnum =strsplit(TBS1052Binfo{[57]},",")
selnum =strsplit(TBS1052Binfo{[61]},",")
firmver =strsplit(TBS1052Binfo{[65]},",")
source_ch =sourcech{[2]}
vertical_unit =verut{[2]}
vertical_scale =verscl{[2]}
vertical_offset =veroff{[2]}
horizontal_unit =horut{[2]}
horizontal_scale=horscl{[2]}
pt_fmt =ptfmt{[2]}
yzero =y0{[2]}
probe_atten =prbatt{[2]}
model_num =mdlnum{[2]}
serial_num =selnum{[2]}
firmware_ver =firmver{[2]}
# サンプリング数(Record Length) 2列の1行目
num_of_smpl=TBS1052Bdata(1,2)
# サンプリング周期(Sample Interval) 2列の2行目
smpl_period=TBS1052Bdata(2,2)
# サンプリング長
smpl_len=TBS1052Bdata(1,2)*TBS1052Bdata(2,2)
# FFT
FFTcomplex=fft( Timeseri(:,2))/num_of_smpl ;
# ゼロシフト(直流分を周波数範囲の中央にもってくる)
ZsFFTcomplex=fftshift( FFTcomplex ) ;
# 周波数軸設定
Ax_freq=linspace(-num_of_smpl/2/(smpl_period*num_of_smpl), (num_of_smpl/2-1)/(smpl_period*num_of_smpl), num_of_smpl);
# 周波数軸を含むFFT行列
AxFFTcomplex=[Ax_freq.' ZsFFTcomplex];
# 等価振幅と位相のFFT行列
AxFFTamph=[ AxFFTcomplex((num_of_smpl/2+1):num_of_smpl,1) abs( AxFFTcomplex((num_of_smpl/2+1):num_of_smpl,2))*2 angle( AxFFTcomplex((num_of_smpl/2+1):num_of_smpl,2))*360/(2*pi)];
#
subplot(3,1,1)
plot(Timeseri(:,1),Timeseri(:,2))
axis([min(Timeseri(:,1)), max(Timeseri(:,1)), -max(abs(Timeseri(:,2)))*1.1, max(abs(Timeseri(:,2)))*1.1])
title(filename)
ylabel(strcat( 'amp ','[',vertical_unit,']'))
xlabel(strcat( 'time ','[',horizontal_unit,']'))
subplot(3,1,2)
plot(AxFFTamph(:,1),AxFFTamph(:,2));
ylabel( strcat( 'amp ','[',vertical_unit,']'))
xlabel( strcat( 'freq ','[1/',horizontal_unit,']'))
axis([ min(AxFFTamph(:,1)), max(AxFFTamph(:,1))])
subplot(3,1,3)
gtestset=plot(AxFFTamph(:,1),AxFFTamph(:,3));
ylabel('phase [deg]')
xlabel( strcat( 'freq ','[1/',horizontal_unit,']'))
axis([ min(AxFFTamph(:,1)), max(AxFFTamph(:,1)) ,-180,180])
#bottom_title(strcat( 'number of smple :',sprintf("%d",num_of_smpl)),'sampling period:',sprintf("%e ",smpl_period),"sec")
GUIというのであれば、アイコンから起動したいですよね。
gnome だとランチャーにスクリプトのパスを設定すればいいのですが、直接 octfft.sh を指定しても動いてくれません。
どうも一旦端末を起動して、そこから octfft.sh を起動する必要があるようです。
そのためのスクリプトが下記。このスクリプトをランチャーで起動すればOK。
#!/bin/shまぁ、この程度ならわざわざスクリプト書かなくても、ランチャーのコマンドに直接入れてもいいでしょうけどね。 ]]>
# TBS 1052B のデータを FFT して表示
/usr/bin/gnome-terminal -e "/home/hoge/fuga/octfftz.sh"
そんなQCADですが、使っている上で、意外と情報が見つからなかったものがあったので以下にメモります。
QCADで形状データから寸法を表示させる場合で、直径記号のΦを寸法値の前に入れたい場合は、
対象の寸法線を選択し、プロパティエディタの「Specific Properties」の「ラベル」に"%%c<>"と入れる。
(下図参照)
"%%c"は直径記号の「Φ」をつけるための呪文で"<>"は形状データからの読み取り値のようです。
ここらへんはAutoCADと同じらしいのですが、AutoCADは使ったことがないので・・・・
ちなみに角度の「°」は"%%d"のようです。
以下はQCADのスクリーンショット
げっ。これまでの作業が水の泡???。あきらめないとダメ?
と思ったけど、Fritzing ってユーザデータは XML なんですよね。
ユーザファイルの *.fzz というのはどうも zip 形式のファイルのようなので、
解凍して *.fz ファイル(XML形式のテキストファイル)を読みだして、
関連しそうな所(今回の私の事例では新しく作ったパーツの部分)を探して、
削除したら何とか立ち上がってくれた。助かった。^^;
削除した場所は関連するパーツの名称のある <instance> タグの
スコープ全部を削除しました。
悪さしてそうな<instance>から</instance>まで全部消す。
<instance moduleIdRef="hogehoge" ・・・・
・
・
</instance>
消した後、XMLファイル(*.fz)をzipファイルの*.fzzに戻してやって、
Fritzing で読み直したら悪くない部分がなんとか生き残っていました。^^;
今回使った Fritzing はバージョン0.8.7。
OS は Vine Linux 6.3 です。
]]>