TBS1052B + octave 4.0.1 ( + Vine linux 6.3 ) [Linux]
Vine linux でもOctave が3.6.4 から4.0.1 へアップグレードされましたね。
パッケージを作ってくださった方に感謝します。
FFT結果を判りやすく表示させるために findpeaks の関数を使いたかったのですが、
旧版ではサポートされておらず、あきらめいていました。
が、今回、新しいバージョンがリリースされ、findpeaksも使えるようになったので、
早速アップグレードした次第です。
FFTについては前の記事でTBS1052B用のoctaveスクリプトを作っていて、
これは新バージョンでもそのまま動作しましたが、
Octave を起動するためのbashスクリプトはそのままでは動きませんでした。
そういった訳で、ほんの少し引っかかった部分があったのでメモします。
3.6.4では octave の起動オプションに --eval とスクリプトファイル指定の併用が可能でしたが、
4.0 .1ではどうもだめっぽい。
bash のスクリプトで取得したデータのファイル名を octave に渡す必要があるため、
--eval オプションの方でデータのファイル名とoctave スクリプトの両方をを渡す形にしました。
そのための スクリプトは下記。スクリプト中の filename は TBS1052B で測定してUSB に保存した
波形データ(CSVファイル)で、octfft.m(後述の octave スクリプトファイル)の中で使用しています。
#!/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アナライザとして一通りの機能をもったものになります。
ここまで来るとカーソルで値を読み取る機能も欲しいですが、
それはちょっとめんどくさそう・・・
コメント 0