So-net無料ブログ作成
  • ブログをはじめる
  • ログイン

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-fft.png
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アナライザとして一通りの機能をもったものになります。
ここまで来るとカーソルで値を読み取る機能も欲しいですが、
それはちょっとめんどくさそう・・・


nice!(0)  コメント(0)  トラックバック(0) 
共通テーマ:趣味・カルチャー

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

※ブログオーナーが承認したコメントのみ表示されます。

トラックバック 0

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。

×

この広告は180日以上新しい記事の更新がないブログに表示されております。