信号処理

Signal-Processing

Originは、信号処理、ノイズデータのスムージングからFFT、逆FFT、ショートタイムFFT、コンボリューションと相関、FFtフィルタリング、ウェーブレット分析のXファンクションとLabTalk関数を提供しています。

これらのXファンクションは、信号処理カテゴリーで利用でき、次のコマンドを実行して、一覧表示できます。

lx cat:="signal processing*";

ショートタイムFFTやウェーブレットなどのいくつかの機能は、OriginProでのみ利用できます。

信号処理のLabTalk関数は、信号処理関数で使用でき、複数のデータセットの振幅、位相などのFFT結果を個別に計算し、結果を目的の列に配置するために使用されます。一方、これらの関数を使用して、異なるデータセット間でDCオフセット除去振幅を比較できます。

次のセクションでは、スクリプトから信号処理のXファンクションとLabTalk関数を呼び出すいくつかの簡単なサンプルを提供しています。

スムージング

ノイズデータのスムージング は、 smooth Xファンクションを使って実行することができます。

// 3次多項式でSavitzkyGolay法を使って
// ワークシートの列1,2のXYデータをスムージング
range r=(1,2); // ワークシートにXYがあり、アクティブである前提
smooth iy:=r meth:=sg poly:=3;

レイヤ内のすべてのプロットをスムージングするため、以下のプロットをループ

// レイヤ内のデータプロットの数を数え、結果を 
// 変数"count"に保存
layer -c;       
// このグラフページの名前を取得      
string gname$ = %H;   
// smoothという名前の新しいブックを作成 - 実際の名前はbkname$に保存されている
newbook na:=Smoothed; 
// 列なしで開始
wks.ncols=0;          
loop(ii,1,count) {
    // 入力範囲は'ii'番目のプロットを参照
    range riy = [gname$]!$(ii);            
    // 出力範囲は2つの新しい列を参照
   range roy = [bkname$]!($(ii*2-1),$(ii*2)); 
    // 3次多項式を使ったSavitsky-Golay スムージング
    smooth iy:=riy meth:=sg poly:=3 oy:=roy;   
}

FFT および IFFT

次のサンプルでは、信号処理関数にあるLabTalk関数によってFFT結果を個別に取得する方法を示しています。

newbook;
// 信号データをインポート
string fname$ = system.path.program$ + "Samples\Signal Processing\fftfilter1.DAT";
impASC fname:=fname$;

// 5列を追加して、異なる値を保存
worksheet -a 5;
// 区別するために列ラベルを設定
col(B)[L]$ = "Raw Signal";

// FFT分析からの複素数の結果を保存する前に、データ型をcomplexに設定
wks.col3.numerictype = 11;
col(C)[L]$ = "FFT Complex";

// FFT複素数結果を計算
col(C) = fftc(col(B)); //シフトなし

col(D)[L]$ = "Frequency";
wks.col4.type = 4; //X列に設定

// 列Aの時間データに基づいて周波数を計算
col(D) = fftfreq(col(A)[2]-col(A)[1], wks.col1.nrows);

col(E)[L]$ = "Magnitude";

// FFTマグニチュード結果
col(E) = fftmag(col(C)); //  fftmag(col(C), 2)で両側のシフト付きマグニチュードを取得

col(F)[L]$ = "Phase";
col(F) = fftphase(col(C)); // fftphase(col(C), 2, 1, 0)でラジアン単位で両側の巻きのない位相を取得

wks.col7.numerictype = 11;
col(G)[L]$ = "Shifted FFT Complex";
col(G) = fftshift(col(C));

// 計算結果のスパークラインを更新
sparklines sel:=0 c1:=4 c2:=6;
// コンテンツに合わせて列幅を自動調整
wautosize;

次の例は、シフト付きFFT複素数結果からIFFT結果を取得する方法を示しています。

//上記の例から続けて、G列のFFT結果をシフトしたと仮定
// さらに2列追加して、異なる値を保存
worksheet -a 2;
// 区別するために列にラベル付け
col(H)[L]$ = "Unshifted FFT Complex";
// データ型をcomplexに設定
wks.col8.numerictype = 11;
// IFFTを実行する前に、シフトしたFFT結果のシフトを解除
col(H) = ifftshift(col(G));

col(I)[L]$ = "IFFT result";
wks.col9.numerictype = 11;
// シフトされていないFFT結果から逆FFTを計算
col(I) = invfft( col(H) );

以下のサンプルでは、FFTの際にウィンドウを指定する方法を示します。

newbook;
// 信号データをインポート
string fname$ = system.path.program$ + "Samples\Signal Processing\Chirp Signal.dat.";
impASC fname:=fname$;

// 入力信号の範囲変数を定義
range rSignal = col(B);
int wSize = rSignal.nrows;

col(C)[L]$ = "Apply Blackman window";
col(C) = col(B) * windata(6, wSize); // Blackmanウィンドウを適用

col(D)[L]$ = "Amplitude";
col(D) = fftamp( fftc(col(C)) ); //Originウィンドウ修正なし

次のサンプルは、ループを使用して複数の列でFFT分析を実行し、取得した振幅の列と位相の列を並べて配置する方法を示しています。ワークシートで複数の列を選択し、右クリックして複数列の値を設定...を選択することにより、ワークシートで複数の列のFFT結果を直接計算することもできます。

//FFT分析用に複数の列データを準備
newbook;
int nc = 5, ii;
wks.ncols = 4*nc+2;
//X列に時間データを入力
col(A) = data(0, 1-1e-3, 1e-3);
//5列をf1 Hz正弦波とf2 Hz正弦波の合計を入力
for( ii = 1; ii<=nc; ii++ )
{
double f1, f2;
f1 = 50 + 20*ii;
f2 = 100 + 20*ii;
wcol(ii+1)[L]$ = "Data$(ii)";
wcol(ii+1) = (2+ii)*sin(2*pi*f1*col(A))+(8-ii)*sin(2*pi*f2*col(A))+rnd();
}

//FFTの前に平均値を減算してDCオフセットを除去
for( ii = 1; ii<=nc; ii++ )
{
wcol(nc+ii+1)[C]$ = "Subtract mean of Data$(ii)";
wcol(nc+ii+1) = wcol(ii+1)-mean( wcol(ii+1) );
}

//FFT結果の周波数を計算
wks.col$(2*nc+2).type = 4; //X列に設定
wcol(2*nc+2)[L]$ = "Frequency";
wcol(2*nc+2) = fftfreq( col(A)[2]-col(A)[1], wks.col1.nrows );

//5列のFFTの振幅と位相を計算
for( ii = 1; ii<=nc; ii++ )
{
//FFTの振幅
wcol(2*nc+ii+2)[L]$ = "Amplitude";
wcol(2*nc+ii+2)[C]$ = "Data$(ii)";
wcol(2*nc+ii+2) = fftamp( fftc( wcol(nc+ii+1) ) );

//FFTの位相
wcol(3*nc+ii+2)[L]$ = "Phase";
wcol(3*nc+ii+2)[C]$ = "Data$(ii)";
wcol(3*nc+ii+2) = fftphase( fftc( wcol(nc+ii+1) ) );
}

//2つのレイヤに5列のFFT振幅と位相の結果をプロット
plotstack iy:=($(2*nc+2),$(2*nc+3):$(4*nc+2)) portrait:=0 order:=0 layer:=2 number:="5 5";

FFTとフィルタリング

次のサンプルは、 fft1 Xファンクションを使ったデータの1D FFTを実行する方法を示しています。

// サンプルファイルをインポート
newbook;
fname$ = system.path.program$ + "Samples\Signal Processing\fftfilter1.dat";
impasc;
// FFTを実行し、名前を付けたツリーに出力
Tree myfft;
fft1 ix:=2 rd:=myfft rt:=<none>; 
// コマンド list vtを使ってすべてのツリーを一覧表示

ツリーに結果を格納すると、その出力に対して次のように分析を行うことができます。

// 目的のツリーベクターノードをデータセットにコピー
// ピークを位置づけ、周波数成分を平均
dataset tmp_x=myfft.fft.freq;
dataset tmp_y=myfft.fft.amp;
// statsを実行し、結果を出力
percentile = {0:10:100};
diststats iy:=(tmp_x, tmp_y) percent:=percentile;
type "The mean frequency is $(diststats.mean)";

次のサンプルは、 fft_filters Xファンクションを使ったデータの信号 フィルタリング を実行する方法を示しています。

// ノイズ入りのデータをインポートする
newbook;
string filepath$ = "Samples\Signal Processing\";
string filename$ = "Signal with High Frequency Noise.dat";
fname$ = system.path.program$ + filepath$ + filename$;
impasc;
plotxy iy:=(1,2) plot:=line;
 
// ローパスフィルタを実行
fft_filters filter:=lowpass cutoff:=1.5;