LabTalkによる信号処理

 

概要

Originは、信号処理、ノイズデータのスムージングからFFT、短時間 FFT（STFT）、コンボリューションと相関、FFTフィルタリング、ウェーブレット分析といったXファンクションを提供しています。このチュートリアルでは、データに対して1D FFTを実行するサンプルを紹介します。

学習する項目

このチュートリアルでは、LabTalk スクリプトを使用して次のことを行う方法を説明します。

  • 事前に定義したフィルタでファイルをインポート
  • 複数のファイルに対して FFT (高速フーリエ変換) を実行
  • すべてのファイルをループ
  • FFTピークの位置をサマリーシートに保存し、ファイル名からの変数の関数としてピーク位置をプロット
  • 線形フィットを実行し、傾きと切片の値を返す

 

ステップ

結果シートを作成

以下のスクリプトを実行してプロジェクトエクスプローラのディレクトリを変更し、pe_cdおよびnewbook Xファンクションを使用して新しいワークブックを開きます。次に、列の範囲を定義し、ロングネームと単位を割り当てます。

// プロジェクトエクスプローラフォルダをルートに設定
pe_cd /;
// FFT分析の結果を保持する新しいブックを作成
// 新しいブックのショートネームを文字列 '''ResultBook''' に保存
newbook result:=ResultBook$;

// 結果を後で保存するために、列1と2の範囲を定義
//複数の範囲変数をカンマで区切って同じ行に宣言できます。
range rx = 1, rp = 2;
// 列1 と 2のロングネームと単位を設定
rx[L]$ = "X";
rx[U]$ = "mm";
rp[L]$ = "Peak Frequency";
rp[U]$ = "kHz";

文字列でファイルを検索

<Origin EXE Folder>\Samples\Signal Processing\にあるファイルをワイルドカードとfindfiles Xファンクションを使って検索します。そしてファイル数を数えます。

// データファイルがある場所を指定
string path$ = system.path.program$ + "Samples\Signal Processing\";
// ''TR*.dat''というすべてのファイルを検索
findfiles ext:="TR*.dat";
// ファイル数取得。LF は改行のこと
int numFiles = fname.GetNumTokens(LF);

ループで複数ファイルのFFTを実行

1 つのファイルで FFT を実行
TR*.dat ファイルを見つけたら、それらをインポートしてFFTを実行し、FFTピークの位置をサマリーシートに保存します。このサンプルでは、1つのデータファイルに対してこれらの操作を実行する手順を詳しく説明します。

  1. データフォルダにあるインポートフィルタTR Data Files.oif を使用して、ファイル名とインポートデータファイルを取得します。フィルタは、 frequency列を作成して値を設定することで、データの後処理を行います。
    この場合、Xファンクションpe_mkdir および impFileが使用されます。 
    // ''ifile'' 整数変数をファイル インデックスとして宣言
// 最初のファイルでは、この変数を1に設定
int ifile = 1;

// 2 つの文字列変数を作成
string filepath$, file$;      
// ファイル名を取得
filepath$=fname.gettoken(ifile,LF)$;
// パスなし、拡張子なしのファイル名だけを解析
file$ = filepath.gettoken(filepath.getnumtokens(\),\)$;
file$ = file.gettoken(1,.)$;
        
// 新しいプロジェクトエクスプローラサブフォルダーを作成し、アクティブに設定
pe_mkdir file$ cd:=1;
        
// データファイルをインポートするための新しいブックを追加
newbook;
// データフォルダに存在する定義済みフィルタを含むファイルをインポート
impfile fname:=filepath$ filtername:="TR Data Files.oif" location:=0;
  2. FFTを実行します。 
    // インポートされたデータの列 2 と 3 で デフォルト設定でFFTを実行
fft1 (2,3);
  3. 最大ピーク位置を見つけ、関連する周波数値を出力し、ファイル名を解析して ResultBook に渡します。
    wks.ncolsを使用してワークシートの列数を取得し、limitコマンドでデータセットの制限値を見つけます。 
    // FFT結果シートをアクティブにする - シート2
page.active = 2;
// ワークシートに新しい列を追加
wks.addcol();
// 新しい列はこのシートの最後の列なので
// 列のインデックスはワークシートの列数と等しい
// ワークシートの列数で新しい整数変数 ''nmag'' を宣言
int nmag = wks.ncols;
// 列 には複素データが含まれ、それを使用してマグニチュードを計算し、ピーク位置を検索
// 列 2 の絶対値を最後の列に入力
wcol(nmag) = abs(col(2));

// limitコマンドを使用して、最後の列の最大ピーク位置を検索
limit wcol(nmag);
// 関連付けられた周波数値を列1からPEフォルダのResultBookシート1の列2に保存
// limit.imax最大Y値に対応するインデックス
rp[ifile] = col(1)[$(limit.imax)];
// 最終列を削除
delete wcol(nmag);
        
// ファイル名を解析し、文字列値 ''aa'' を取得
string aa$ = file.gettoken(2,'R')$;
aa$=aa.gettoken(1,'M')$;
// ''aa''値を ResultBook シート1列1に保存
rx[ifile] = %(aa$);

// PE フォルダを前のレベルに設定
pe_cd ..;


複数のファイルに対してFFTを実行
複数のファイルに対して操作を繰り返すためにforループコマンドを使用します。このセクションでは、ファイルをインポート後FFTを実行し、FFTピークの位置をサマリーシートに保存するための完全なループスクリプトを示します。
Note: スクリプトウィンドウでループスクリプトを実行する場合は、スクリプト全体を選択状態にしてから、Enterキーを押してください。

for(int ifile = 1; ifile <= numFiles; ifile++)
{
        string filepath$, file$;
        filepath$=fname.gettoken(ifile,LF)$;
        file$ = filepath.gettoken(filepath.getnumtokens(\),\)$;
        file$ = file.gettoken(1,.)$;
        pe_mkdir file$ cd:=1;
        newbook;
        impfile fname:=filepath$ filtername:="TR Data Files.oif" location:=0;
        
        fft1 (2,3);
        
        page.active = 2;
        wks.addcol();
        int nmag = wks.ncols;
        wcol(nmag) = abs(col(2));
        limit wcol(nmag);
        rp[ifile] = col(1)[$(limit.imax)];
        delete wcol(nmag);
        string aa$ = file.gettoken(2,'R')$;
        aa$=aa.gettoken(1,'M')$;
        rx[ifile] = %(aa$);
        
        pe_cd ..; 
}

結果シートで新しいグラフを作成

すべてのファイルに対してFFTを実行した後、plotxy Xファンクションを使用して、ファイル名から抽出されたX値とピーク周波数値を使用して散布図を作成します。

// ResultBook をアクティブにする
win -a %(ResultBook$);
// 列1と2でプロット
plotxy (1,2);

線形フィット

結果データに対して線形回帰を実行し、傾きと切片を出力します。この場合、 LRコマンドを使用します。

//ResultBookのシート1列2で線形フィットを実行
range bb = [ResultBook$]1!2;
Lr bb;
//''lr.b'' は傾き
type "Slope is $(lr.b)"; 
//''lr.a'' は切片
type "Intercept is $(lr.a)";
English | Deutsch | 日本語