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; }
次のサンプルでは、信号処理関数にある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";
次のサンプルは、 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;