第4章 ディジタル信号処理パッケージ

1. ISPP の概要

ISPP (Interactive Signal Processing Package) は,SATELLITE の中核をなすモジュールで,ウィンドウ等の前処理,FFT,線形予測によるスペクトル推定,フィルタリング,ケプストラム解析等,ディジタル信号処理の方法論を網羅する処理機能が関数として実装されています.こうした様々な関数は,全て複数のデータ格納領域 (バッファ) またはデータファイルを対象として処理を実行することができ,関数を組み合わせて処理を記述することによって,計測されたデータの信号処理,統計処理をはじめ多角的なデータ解析を行うことができます.

2. ISPP のコマンド体系

ISPP の関数は 表 4.1. 「ISPP 関数一覧」 に示すように,データの生成,操作,補間,算術演算,解析に分類されます.こうした基本的な機能を有した関数を組み合わせることによって,複雑なデータ解析も行うことができます.

表 4.1. ISPP 関数一覧

データ生成
argenAR モデルによるデータを生成する
gauss22 次元ガウス関数を生成する
arand任意の確率分布に従う乱数データを生成する
mnrand多次元正規乱数データを生成する
nrand正規乱数データを生成する
urand一様乱数データを生成する
データ操作
dccutバッファの直流成分を除去する
norm指定したバッファの内容を正規化する
データ補間
interp22 次元データの補間を行う
算術演算
det2 次元データを行列とみなし,行列式を求める
eigen2 次元データを行列とみなし,固有値・固有ベクトルを求める
sumバッファの内容を積分する
inv2 次元データを行列とみなし,逆行列を求める
mul2 つの 2 次元データを行列とみなし,積を求める
trans2 次元データを行列とみなし,転置行列を求める
nmeq正規方程式を解く
データ解析
burgBurg 法よりパワースペクトルと AIC を求める
cepデータの複素ケプストラムを求める
fftc複素数高速 (逆) フーリエ変換を求める
fftn多次元データに対しての複素数高速 (逆) フーリエ変換を求める
firデータに対して FIR 型フィルタ処理を行う
firmakeFIR 型フィルタを設計する
hilヒルベルト変換を行う
butwmakeバタワース特性の IIR 型フィルタを設計する
icep逆ケプストラム解析を行う
iirデータに対して IIR 型フィルタ処理を行う
iirmake零点,極から IIR 型フィルタの係数を求める
levinLevinson-Durbin 法によりパワースペクトルと AIC を求める
phase2 つのバッファの内容により位相を求める
poleAR 係数から極の位置を求める
power2 つのバッファの内容を 2 乗して加算する
rankデータの度数分布,ガウス分布を求める
spcfデータのパワースペクトル,位相を求める
windowデータに対してウィンドウ処理を行う

3. 使用例

ISPP の基本的な使用例として,主にフーリエ変換,フィルタ処理,行列演算について,以下にその例を示します.

3.1. フーリエ変換

ここでは,2 つの正弦波を合成した信号に雑音を重畳させた信号をフーリエ変換する一連の手順を,実際の計算結果とともに示します.

データの生成

正弦波の生成

[]SATELLITE[]~/home/demo:[1]% sampling(4000);
[]SATELLITE[]~/home/demo:[2]% t = (0~1999)/2000;

サンプリング周波数,時刻データを設定し,series 変数 (t) に格納します.

[]SATELLITE[]~/home/demo:[3]% a = 5*sin(2*PI*20*t) + 3;
[]SATELLITE[]~/home/demo:[4]% b = 3*sin(2*PI*50*t) + 3;

周波数,振幅が異なり,直流成分を持つ正弦波を series 変数 (a, b) に格納します.

データ切り出し生成された正弦波の 0 点目から 1023 点目までの,計 1024 点のデータを切り出します.

acut = 
bcut = 

切り出したデータを格納する series 変数を決めます.

acut = cut(a,
bcut = cut(b,

切り出すデータを設定します.

acut = cut(a,0
bcut = cut(b,0

切り出すデータの始点を設定します.

[]SATELLITE[]~/home/demo:[5]% acut = cut(a,0,1023);
[]SATELLITE[]~/home/demo:[6]% bcut = cut(b,0,1023);

切り出すデータの終点を設定します.

以上で切り出したデータが series 変数 (acut, bcut) に格納されます.

乱数の生成

正規乱数を発生させる場合には,nrand 関数を用います.

nois = 

発生させた乱数を格納する series 変数を決めます.

nois = nrand(1024,

発生させるデータ点数を設定します.

nois = nrand(1024,1

乱数の初期値 (奇数) を設定します.

nois = nrand(1024,1,0

乱数の平均値を設定します.

[]SATELLITE[]~/home/demo:[7]% nois = nrand(1024,1,0,1);

乱数の分散を設定します.

以上でデータ点数 1024 点,平均 0,分散 1 の正規乱数データが series 変数 (nois) に格納されます.なお,一様乱数を発生させる場合には,urand 関数を用います.

信号の合成

以上の作業により得られた信号を合成します.

data = 

合成信号を格納する series 変数を決めます.

[]SATELLITE[]~/home/demo:[8]% data = acut + bcut + nois;

各信号を合成します.

以上で 2 つの正弦波データと正規乱数データを合成したデータが series 変数 (data) に格納されます.

[]SATELLITE[]~/home/demo:[9]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[10]% x = 0~1023;
[]SATELLITE[]~/home/demo:[11]% size(100,50);
[]SATELLITE[]~/home/demo:[12]% scale("N","F","N","F",0,1000,-20,20);
[]SATELLITE[]~/home/demo:[13]% graph(data,x,0,0,0,0,0);
[]SATELLITE[]~/home/demo:[14]% title(1,"Time[sec]","Amplitude");
[]SATELLITE[]~/home/demo:[15]% axis(1,1,"XY","XY",5,0,3,200,20,1);
[]SATELLITE[]~/home/demo:[16]% frame();

合成した信号の波形を表示します.これらの関数はグラフィック・システム (GPM) で詳しく説明します.

合成した信号の波形を 図 4.1. 「2 つの正弦波および雑音を足し合わせた信号の波形」 に示します.

2 つの正弦波および雑音を足し合わせた信号の波形

図 4.1. 2 つの正弦波および雑音を足し合わせた信号の波形

データの前処理

ある信号に FFT を施す前に行うべき前処理として,以下に直流分除去とウィンドウ処理の方法を示します.

直流成分除去

データの直流成分を除去する場合には,dccut 関数を用います.

data1 =

直流成分を除去したデータを格納する series 変数を決めます.

[]SATELLITE[]~/home/demo:[17]% data1 = dccut(data);

直流成分除去対象データを設定します.

以上で直流成分が除去されたデータが series 変数 (data1) に格納されます.

直流成分を除去した信号の表示

直流成分を除去した信号波形を表示します.

[]SATELLITE[]~/home/demo:[18]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[19]% size(100,50);
[]SATELLITE[]~/home/demo:[20]% scale("N","F","N","F",0,1000,-20,20);
[]SATELLITE[]~/home/demo:[21]% graph(data1,x,0,0,0,0,0);
[]SATELLITE[]~/home/demo:[22]% title(1,"Time[sec]","Amplitude");
[]SATELLITE[]~/home/demo:[23]% axis(1,1,"XY","XY",5,0,3,200,20,1);
[]SATELLITE[]~/home/demo:[24]% frame();

直流成分を除去した信号を表示します.直流成分を除去した信号波形を 図 4.2. 「直流成分除去後の波形」 に示します.

直流成分除去後の波形

図 4.2. 直流成分除去後の波形

ウィンドウ処理

データにウィンドウ処理を施す場合には,window 関数を用います.

data2 =

ウィンドウ処理後のデータを格納する series 変数を決めます.

data2 = window(data1,

ウィンドウ処理対象データを設定します.

data2 = window(data1,1,

ウィンドウの種類 (1 : ハミング窓,2 : ハニング窓,3 : ブラックマン窓,4 : トライアングル窓) を設定します.

[]SATELLITE[]~/home/demo:[25]% data2 = window(data1,1,0);

データの補正 (補正を行うと,ウィンドウ処理前後の積分値が等しくなります) をしない場合は0,する場合は 1 と設定します.

以上でウィンドウ処理されたデータが series 変数 (data2) に格納されます.

ウィンドウ処理した信号の表示

ウィンドウ処理をした信号波形を表示します.

[]SATELLITE[]~/home/demo:[26]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[27]% size(100,50);
[]SATELLITE[]~/home/demo:[28]% scale("N","F","N","F",0,1000,-20,20);
[]SATELLITE[]~/home/demo:[29]% graph(data2,x,0,0,0,0,0);
[]SATELLITE[]~/home/demo:[30]% title(1,"Time[sec]","Amplitude");
[]SATELLITE[]~/home/demo:[31]% axis(1,1,"XY","XY",5,0,3,200,20,1);
[]SATELLITE[]~/home/demo:[32]% frame();

ウィンドウ処理後の信号波形を 図 4.3. 「ウィンドウ処理後の波形」 に示します.

ウィンドウ処理後の波形

図 4.3. ウィンドウ処理後の波形

フーリエ変換

fftc 関数により,複素数データのフーリエ変換および,逆フーリエ変換を行うことができます.なお,フーリエ変換には FFT 法を用いるため,変換対象データの点数は,2 のべき乗である必要があります.

[]SATELLITE[]~/home/demo:[33]% series Rout,Iout;

出力時系列の実部と虚部を格納する series 変数を宣言します.

[]SATELLITE[]~/home/demo:[34]% rei = (0~1023)*0;

入力信号に虚部データが存在しない場合 0 を入力するため,あらかじめ series 変数 (rei) に 1024 点の 0 を格納します.

fftc("P",

計算方式フラグ (P : フーリエ変換,I : 逆フーリエ変換) を設定します.

fftc("P",data2,

変換対象データの実部を入力します.

fftc("P",data2,rei,

変換対象データの虚部を入力します.

fftc("P",data2,rei,Rout,

フーリエ変換後の実部データを格納する series 変数を設定します.

[]SATELLITE[]~/home/demo:[35]% fftc("P",data2,rei,Rout,Iout);

フーリエ変換後の虚部データを格納する series 変数を設定します.

以上でフーリエ変換されたデータが series 変数 (Rout, Iout) に格納されます.

[]SATELLITE[]~/home/demo:[36]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[37]% size(100,50);
[]SATELLITE[]~/home/demo:[38]% scale("N","F","N","F",0,100,-1000,1000);
[]SATELLITE[]~/home/demo:[39]% graph(Rout,"F",0,0,0,0,0);
[]SATELLITE[]~/home/demo:[40]% ltype(5,1);
[]SATELLITE[]~/home/demo:[41]% graph(Rout,"F",0,0,0,0,0);
[]SATELLITE[]~/home/demo:[42]% title(1,"Frequency[sec]","Amplitude");
[]SATELLITE[]~/home/demo:[43]% axis(1,1,"XY","XY",5,0,3,20,1000,1);
[]SATELLITE[]~/home/demo:[44]% frame();

フーリエ変換後の信号を表示します.

フーリエ変換後の信号波形を 図 4.4. 「FFT 処理後の信号波形 (実部 : 実線,虚部 : 破線)」 に示します.

FFT 処理後の信号波形 (実部 : 実線,虚部 : 破線)

図 4.4. FFT 処理後の信号波形 (実部 : 実線,虚部 : 破線)

パワースペクトル

2 つの変数の内容を 2 乗して加算する場合には,power 関数を用います.この関数により,fftc 関数により得られたフーリエ変換後の実部データおよび虚部データから,入力時系列のパワースペクトルを求めることができます.また,入力時系列から直接パワースペクトルを求めることのできる,spcf 関数もあります.spcf 関数については,後に説明します.

pw =

求めたパワースペクトルを格納する series 変数を決めます.

pw = power(Rout,

フーリエ変換後の実部データを入力します.

[]SATELLITE[]~/home/demo:[45]% pw = power(Rout,Iout);

フーリエ変換後の虚部データを入力します.

以上で,パワースペクトルデータが series 変数 (pw) に格納されます.

[]SATELLITE[]~/home/demo:[46]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[47]% size(100,50);
[]SATELLITE[]~/home/demo:[48]% scale("N","F","N","F",0,100,-200000,2000000);
[]SATELLITE[]~/home/demo:[49]% ltype(1,1);
[]SATELLITE[]~/home/demo:[50]% graph(pw,"F",0,0,0,0,0);
[]SATELLITE[]~/home/demo:[51]% title(1,"Frequency[sec]","Power");
[]SATELLITE[]~/home/demo:[52]% axis(1,1,"XY","XY",5,0,3,20,1000000,0);
[]SATELLITE[]~/home/demo:[53]% frame();

求めたパワースペクトルを表示します.

求めたパワースペクトルを 図 4.5. 「パワースペクトル」 に示します.

パワースペクトル

図 4.5. パワースペクトル

位相特性

2 つの変数の内容により位相を求める場合には,phase 関数を用います.この関数により,fftc 関数により得られたフーリエ変換後の実部データおよび虚部データから,その位相特性を求めることができます.

phs = 

求めた位相を格納する series 変数を決めます.

phs = phase(Rout,

フーリエ変換後の実部データを入力します.

phs = phase(Rout,Iout

フーリエ変換後の虚部データを入力します.

phs = phase(Rout,Iout,"D",

出力データ形式 ("D" : degree,"O" : radian) を設定します.

[]SATELLITE[]~/home/demo:[54]% phs = phase(Rout,Iout,"D","U");

位相戻しをする場合は "U",しない場合は "O" と設定します.

以上で位相データが series 変数 (phs) に格納されます.

入力系列からパワースペクトルと位相を 1つの関数で求める方法

spcf 関数により,ある入力時系列から,そのパワースペクトルおよび位相を,同時に求めることができます.以下にその具体的な手順を示します.

[]SATELLITE[]~/home/demo:[55]% series pw,phs;

パワースペクトルおよび位相を格納する series 変数を宣言します.

spcf(data2,

対象時系列を設定します.

spcf(data2,pw,

求めたパワースペクトルを格納する series 変数を設定します.

[]SATELLITE[]~/home/demo:[56]% spcf(data2,pw,phs);

求めた位相を格納する series 変数を設定します.

以上でパワースペクトル,位相データが series 変数 (pw,phs) に格納されます.

3.2. フィルタ処理

MAフィルタ

移動平均法を用いたフィルタにより,信号を平滑化する手順を以下に示します.なおフィルタの次数は 5 次とします.

[]SATELLITE[]~/home/demo:[1]% coef = (1/5,1/5,1/5,1/5,1/5);

フィルタ係数を series 変数に格納します.

output =

平滑化された信号データを格納する変数を決めます.

output = fir(coef,

fir 関数により,時系列信号を平滑化します.まず,フィルタの係数を設定します.

[]SATELLITE[]~/home/demo:[2]% output = fir(coef, input);

平滑化の対象となる時系列信号 (input) を設定します.

以上で平滑化されたデータが series 変数 (output) に格納されます.なお,fir 関数は,図 4.6. 「fir 関数のフィルタ処理概念図」 ( は,フィルタ係数) の様な構成で処理を行います.したがって,そのフィルタ係数となるパラメータ系列の長さは,奇数値 (以下,2n+1 とする) でなければなりません.また,処理対象となる時系列データの最初から n 点までのデータについては,適切な処理を行うことができないため,n+1 点目のデータを採用します.時系列データの最後から n+1 点までのデータについても同様で,最後から n+1 点目のデータを採用します.

fir 関数のフィルタ処理概念図

図 4.6. fir 関数のフィルタ処理概念図

FIRフィルタ

ローパス,ハイパス,バンドパスの FIR ファイルを,窓関数法により設計する場合には,まず firmake 関数によりフィルタ係数を求め,fir 関数によりフィルタリングを行います.以下に,ローパスフィルタの場合について,その具体的な手順を示します.

[]SATELLITE[]~/home/demo:[3]% sampling(1024);

サンプリング周波数を設定します.

coef =

求めたフィルタ係数を格納する series 変数を入力します.

coef = firmake(1,

設計するフィルタ特性をローパスフィルタに設定します.(1 : LowPass,2 : HighPass,3 : BandPass)

coef = firmake(1,11

フィルタの次数を設定します.ただし,実際のフィルタリングには fir 関数を用いるため,次数は奇数値でなければなりません.

coef = firmake(1,11,100

遮断周波数を設定します.

[]SATELLITE[]~/home/demo:[4]% coef = firmake(1,11,100,3);

使用する窓関数の種類を設定します.(0 : 矩形窓,1 : ハニング窓,2 : ハミング窓,3 : ブラックマン窓,4 : カイザー窓)

以上で,遮断周波数 100Hz の 11 次 FIR 型ローパスフィルタのフィルタ係数が series 変数 (coef) に格納されます.なお,ハイパスの場合には,同様に,

coef = firmake(2,11,400,3);

と設定できます.この例では,遮断周波数 400Hz の 11 次 FIR 型ハイパスフィルタが設計されます.また,バンドパスの場合には,遮断周波数を 2 つ設定する必要があるので,series で遮断周波数を設定します.例えば,遮断周波数が 100Hz,400Hz で,フィルタ次数が 11 次の場合には,

coef = firmake(3,11,(100,400),3);

とします.

次に,coef を用いて fir 関数によりフィルタリングします.

[]SATELLITE[]~/home/demo:[5]% output = fir(coef,input);

IIRフィルタ

バタワース特性を持った,ローパス,ハイパス,バンドパスに IIR フィルタを設計する場合には,まず,butwmake 関数により伝達関数の零点,極,利得を求めます.次に,iirmake 関数により零点,極をフィルタ係数に変換し,fir,iir 関数によりフィルタリングを行います.以下に,ローパスフィルタの場合について,その具体的な手順を示します.

[]SATELLITE[]~/home/demo:[6]% sampling(1024);

サンプリング周波数を設定します.

[]SATELLITE[]~/home/demo:[7]% series zr,zi,pr,pi;

伝達関数の零点,極の座標を格納する series 変数を宣言します.

gain =

設計したフィルタの利得を格納する scalar 変数を入力します.

gain = butwmake(1,100,

ローパスフィルタの指定と遮断周波数を設定します.

gain = butwmake(1,100,13

フィルタの次数を設定します.ただし,FIR フィルタの場合と同様に,奇数値を設定しなければなりません.また,最大次数は 101 次です.

gain = butwmake(1,100,13,zr,zi,

伝達関数の零点の座標 (実部と虚部) を設定します.

[]SATELLITE[]~/home/demo:[8]% gain = butwmake(1,100,13,zr,zi,pr,pi);

伝達関数の極の座標 (実部と虚部) を設定します.

以上でバタワース特性を持った遮断周波数 100Hz の 13 次 IIR 型ローパスフィルタにおいて,その伝達関数の零点 (zr,zi),極 (pr,pi),利得 (gain) が series 変数に格納されます.なお,ハイパスの場合には,同様に,

gain = butwmake(2,400,13,zr,zi,pr,pi);

設計されます.また,バンドパスの場合は遮断周波数を 2 つ設定する必要があるので,例えば遮断周波数が 100Hz,400Hz で,フィルタ次数が 13 次の場合には,

gain = butwmake(3,(100,400),13,zr,zi,pr,pi);

とします.

次に,zr,zi,pr,pi から,iirmake 関数によりフィルタ係数を算出します.

[]SATELLITE[]~/home/demo:[9]% series a,b;

伝達関数の分母の係数,及び分子の係数を格納する series 変数を宣言します.

iirmake(zr,zi,pr,pi,

伝達関数の零点,極の座標を設定します.

iirmake(zr,zi,pr,pi,a,

伝達関数の分母の係数を設定します.

[]SATELLITE[]~/home/demo:[10]% iirmake(zr,zi,pr,pi,a,b);

伝達関数の分子の係数を設定します.

以上でフィルタ係数が fir,iir 関数にそのまま使用できるフォーマットで出力されます.後は,実際にフィルタリングを行います.

[]SATELLITE[]~/home/demo:[11]% temp = fir(b,input);

伝達関数の分子の計算を行います.

[]SATELLITE[]~/home/demo:[12]% output = iir(a,temp)*gain;

伝達関数の分母の計算を行い,利得を乗ずると出力が求まります.

以上で,フィルタリングされた結果が,series 変数 (output) に格納されます.

iir 関数のフィルタ処理概念図

図 4.7. iir 関数のフィルタ処理概念図

実際の計算例

実際に,バタワース特性を持った遮断周波数 100Hz の 13 次 IIR 型ローパスフィルタを設計し,インパルス応答を求めることでフィルタ特性を求めてみます.以下に,その手順を示します.

[]SATELLITE[]~/home/demo:[13]% sampling(512);

サンプリング周波数を設定します.

[]SATELLITE[]~/home/demo:[14]% series zr,zi,pr,pi;
[]SATELLITE[]~/home/demo:[15]% series a,b;
[]SATELLITE[]~/home/demo:[16]% series u,v;

butwmake 関数,iirmake 関数,spcf 関数で必要な変数をそれぞれ宣言します.

[]SATELLITE[]~/home/demo:[17]% delay = 50;
[]SATELLITE[]~/home/demo:[18]% datp = 511;
[]SATELLITE[]~/home/demo:[19]% impulse = (1,(1~datp)*0);
[]SATELLITE[]~/home/demo:[20]% d_impulse = ((0~delay)*0,impulse);

インパルス信号を生成します.

[]SATELLITE[]~/home/demo:[21]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[22]% x = 0~205;
[]SATELLITE[]~/home/demo:[23]% size(100,50);
[]SATELLITE[]~/home/demo:[24]% scale("N","F","N","A",0,205);
[]SATELLITE[]~/home/demo:[25]% graph(impulse,x,0,0,0,0,0);
[]SATELLITE[]~/home/demo:[26]% title(1,"Data Point","Amplitude");
[]SATELLITE[]~/home/demo:[27]% axis(1,1,"XY","XY",5,0,1,50,0.5,0);
[]SATELLITE[]~/home/demo:[28]% origin(120,20);
[]SATELLITE[]~/home/demo:[29]% graph(d_impulse,x,0,0,0,0,0);
[]SATELLITE[]~/home/demo:[30]% title(1,"Data Point","Amplitude");
[]SATELLITE[]~/home/demo:[31]% axis(1,1,"XY","XY",5,0,1,50,0.5,0);

生成したインパルス信号を出力します.

ここで,impulse (図 4.8. 「impulse 信号」) の前に delay 個の 0 を詰め,d_impulse (図 4.9. 「d_impulse 信号」)を生成した理由は,fir 関数の仕様上,フィルタの次数以内の出力が決定できないためです.つまり,インパルス応答を求めるためには,d_impulse をフィルタリング (図 4.10. 「firtemp 信号」) した後に,delay 分だけシフト (図 4.11. 「firinp 信号」) することにより求めます.

また,図 4.6. 「fir 関数のフィルタ処理概念図」 に示した様に,fir 関数は現在の出力を求めるのに未来の入力を用いています.これでは,因果性を持たないので,iirmake 関数は伝達関数の分子の係数に,(フィルタの次数-1) 個の 0 を詰めた形でその係数を出力します.このため,フィルタリングの対象になるのは,((2×フィルタの次数)-1) 個目より後のデータとなります.この場合 2×13-1=25 ですから,25 よりも大きな値 (50) を delay として定義し,前述の方法でフィルタリングすることにより,その正確な出力を求めることができます.

impulse 信号

図 4.8. impulse 信号

d_impulse 信号

図 4.9. d_impulse 信号

firtemp 信号

図 4.10. firtemp 信号

firinp 信号

図 4.11. firinp 信号

[]SATELLITE[]~/home/demo:[32]% gain = butwmake(1,100,13,zr,zi,pr,pi);
[]SATELLITE[]~/home/demo:[33]% iirmake(zr,zi,pr,pi,a,b);

フィルタの伝達関数を求め,フィルタを設計します.

[]SATELLITE[]~/home/demo:[34]% firtemp = fir(b,d_impulse);
[]SATELLITE[]~/home/demo:[35]% firinp = cut(firtemp,delay+1,datp+delay+1);

伝達関数の分子の部分の計算を行い,impulse の前に詰めた 0 を取り除くことにより,delay 分だけシフトします.

[]SATELLITE[]~/home/demo:[36]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[37]% x = 0~205;
[]SATELLITE[]~/home/demo:[38]% size(100,50);
[]SATELLITE[]~/home/demo:[39]% scale("N","F","N","A",0,205);
[]SATELLITE[]~/home/demo:[40]% graph(firtemp,x,0,0,0,0,0);
[]SATELLITE[]~/home/demo:[41]% title(1,"Data Point","Amplitude");
[]SATELLITE[]~/home/demo:[42]% axis(1,1,"XY","XY",5,0,1,50,900,0);
[]SATELLITE[]~/home/demo:[43]% origin(120,20);
[]SATELLITE[]~/home/demo:[44]% graph(firinp,x,0,0,0,0,0);
[]SATELLITE[]~/home/demo:[45]% title(1,"Data Point","Amplitude");
[]SATELLITE[]~/home/demo:[46]% axis(1,1,"XY","XY",5,0,1,50,900,0);

iir 関数によるインパルス応答を表示します.

[]SATELLITE[]~/home/demo:[47]% output = iir(a,firinp)*gain;

伝達関数の分母の部分の計算を行い,利得を乗じ,設計したフィルタのインパルス応答を求めます.

[]SATELLITE[]~/home/demo:[48]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[49]% x = 0~160;
[]SATELLITE[]~/home/demo:[50]% size(100,50);
[]SATELLITE[]~/home/demo:[51]% scale("N","F","N","A",0,160);
[]SATELLITE[]~/home/demo:[52]% graph(output,x,0,0,0,0,0);
[]SATELLITE[]~/home/demo:[53]% title(1,"Data Point","Amplitude");
[]SATELLITE[]~/home/demo:[54]% axis(1,1,"XY","XY",5,0,1,20,0.3,0);

フィルタのインパルス応答を表示します (図 4.12. 「設計したフィルタのインパルス応答」)

設計したフィルタのインパルス応答

図 4.12. 設計したフィルタのインパルス応答

[]SATELLITE[]~/home/demo:[55]% spcf(output,u,v);

output をフーリエ変換し,パワースペクトル (u) を求めることにより,設計したフィルタの振幅特性を求めます.

[]SATELLITE[]~/home/demo:[56]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[57]% size(100,50);
[]SATELLITE[]~/home/demo:[58]% scale("N","A","N","A");
[]SATELLITE[]~/home/demo:[59]% graph(u,"F",0,0,0,0,0);
[]SATELLITE[]~/home/demo:[60]% title(1,"Frequency[Hz]","Power");
[]SATELLITE[]~/home/demo:[61]% axis(1,1,"XY","XY",5,0,1,50,0.000004,0);

設計したフィルタの振幅特性を表示します (図 4.13. 「設計したフィルタの振幅特性」)

設計したフィルタの振幅特性

図 4.13. 設計したフィルタの振幅特性

3.3. 2 次元 series の扱い方

ISPP には,1 次元データの他に 2 次元データを扱うことができる,fir,fftn といった関数があります.これらの関数での 2 次元データの扱い方について説明します.

2次元FIRフィルタ

[]SATELLITE[]~/home/demo:[1]% x = (1~441)*0;
[]SATELLITE[]~/home/demo:[2]% x:[220] = 1;
[]SATELLITE[]~/home/demo:[3]% input = reform(x,(21,21));

まず,2 次元データを作成します.

[]SATELLITE[]~/home/demo:[4]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[5]% size(80,80);
[]SATELLITE[]~/home/demo:[6]% scale("N","A","N","A");
[]SATELLITE[]~/home/demo:[7]% gsolm(input,0.3,0.4,0,0,0,7,0,1,1,"X",1,0,5);

作成したデータを表示します (図 4.14. 「作成した 2 次元データ」).

作成した 2 次元データ

図 4.14. 作成した 2 次元データ

[]SATELLITE[]~/home/demo:[8]% coeftemp = (1~25)*0+1;
[]SATELLITE[]~/home/demo:[9]% coef = reform(coeftemp,(5,5));

次に,フィルタ係数を決定します.ここでは簡単な例として,全ての係数が 1 の場合について考えます.また,係数データは奇数×奇数である必要があります.

output = fir(coef,

fir 関数により,時系列信号を平滑化します.まず,フィルタ係数を設定します.

[]SATELLITE[]~/home/demo:[10]% output = fir(coef,input);

平滑化の対象となる時系列信号 (input) を設定します.

[]SATELLITE[]~/home/demo:[11]% output2 = fir(coef,input,0.0);

このとき,3 つ目の引数として scalar 型の値を与えた場合には,入力データからはみ出した部分はその値で補間し,2 つ目しか引数を与えなかった場合には,入力全体の平均値で補間します.

[]SATELLITE[]~/home/demo:[12]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[13]% size(80,80);
[]SATELLITE[]~/home/demo:[14]% scale("N","A","N","A");
[]SATELLITE[]~/home/demo:[15]% gsolm(output2,0.3,0.4,0,0,0,7,0,1,1,"X",1,0,5);

求めたフィルタの応答を表示します (図 4.15. 「2 次元 FIR フィルタの応答」)

2 次元 FIR フィルタの応答

図 4.15. 2 次元 FIR フィルタの応答

2次元フーリエ変換

[]SATELLITE[]~/home/demo:[16]% series u[32],v[32];

出力時系列の実部と虚部を格納する series 変数を宣言します.

[]SATELLITE[]~/home/demo:[17]% datp = 15;
[]SATELLITE[]~/home/demo:[18]% data = ((1~datp)*0,1,1,(1~datp)*0);
[]SATELLITE[]~/home/demo:[19]% x = mul(trans(data),data);

入力時系列の実部の 2 次元のデータを作成します.

[]SATELLITE[]~/home/demo:[20]% y = x*0;

入力信号に虚部データが存在しない場合 0 を入力するため,あらかじめ series 変数に 1024×1024 点の 0 を格納します.

[]SATELLITE[]~/home/demo:[21]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[22]% size(80,80);
[]SATELLITE[]~/home/demo:[23]% scale("N","A","N","A");
[]SATELLITE[]~/home/demo:[24]% gsolm(x,0.3,0.4,0,0,0,7,0,1,1,"X",1,0,5);

入力時系列の実部のデータを表示します (図 4.16. 「2 次元 FFT への入力波形 (実部)」)

2 次元 FFT への入力波形 (実部)

図 4.16. 2 次元 FFT への入力波形 (実部)

fftn("P",

計算方式フラグ ("P" : フーリエ変換,"I" : 逆フーリエ変換) を設定します.

fftn("P",x,

変換対象データの実部を入力します.

fftn("P",x,y,

変換対象データの虚部を入力します.

fftn("P",x,y,u,

フーリエ変換後の実部データを格納する series 変数を設定します.

[]SATELLITE[]~/home/demo:[25]% fftn("P",x,y,u,v);

フーリエ変換後の虚部データを格納する series 変数を設定します.power 関数により,fftn 関数により得られたフーリエ変換後の実部データおよび虚部データから,入力時系列のパワースペクトルを求めます.

pw2 =

求めたパワースペクトルを格納する series 変数を決めます.

pw2 = power(u,

フーリエ変換後の実部データを入力します.

[]SATELLITE[]~/home/demo:[26]% pw2 = power(u,v);

フーリエ変換後の虚部データを入力します.

[]SATELLITE[]~/home/demo:[27]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[28]% size(80,80);
[]SATELLITE[]~/home/demo:[29]% scale("N","A","N","A");
[]SATELLITE[]~/home/demo:[30]% gsolm(pw2,0.3,0.4,0,0,0,7,0,1,1,"X",1,0,5);

求めたパワースペクトルを表示します (図 4.17. 「パワースペクトル」)

パワースペクトル

図 4.17. パワースペクトル

3.4. 行列演算

ISPP の特徴の一つとして,行列演算が挙げられます.n 次元 1 次方程式の解を,正規方程式を導くことで求めてみます.

一般に連立方程式は,行列の積により表すことができ,n 次元 1 次方程式の場合は,次式のように表すことができます.

式 (4.1) は,以下のような行列演算を行うことにより,解くことができます

正規方程式とは,式 (4.2) の両辺に左側から をかけた,次式で示される式でした.

ここで,行列 は正則行列と仮定すると,

となり,式 (4.3) を解くことで,式 (4.1) の解を求めることができます.では,行列 に,次の様な数値を代入して,実際に ISPP により式 (4.3) の解を求める手順を以下に示します.

ISPP には,行列演算用の関数として,積を求める mul 関数,逆行列を求める inv 関数,転置行列を求める trans 関数などがあり,これらの関数を組み合わせることによって行列演算を行います.

[]SATELLITE[]~/home/demo:[1]% tempx = (1,1,1,1,2,2,2,1,3,3,3,1,2,3,2);
[]SATELLITE[]~/home/demo:[2]% x = reform(tempx,(5,3));

まず行列 を生成します.

[]SATELLITE[]~/home/demo:[3]% tempy = (10,20,25,27,21);
[]SATELLITE[]~/home/demo:[4]% y = reform(tempy,(5,1));

同様にして,行列 を生成します.

[]SATELLITE[]~/home/demo:[5]% xt = trans(x);

の転置行列を求め,演算結果を xt に格納します.

[]SATELLITE[]~/home/demo:[6]% xtx = mul(xt,x);

の積を求め,演算結果を xtx に格納します.

[]SATELLITE[]~/home/demo:[7]% ixtx = inv(xtx);

の逆行列 を求め,演算結果を ixtx に格納します.

[]SATELLITE[]~/home/demo:[8]% ixtxxt = mul(ixtx, xt);

の積を求め,演算結果を ixtxxt に格納します.

[]SATELLITE[]~/home/demo:[9]% theta = mul(ixtxxt,y);

の積を求めれば,解が求まります.

なお,次のように記述すれば,行列 と行列 から直接に を求めることができます.

theta = mul(mul(inv(mul(trans(x),x)),trans(x)),y);

以上の演算の結果,

と求まります.

3.5. データの統計処理

ここで,正規乱数を発生する ISPP の nrand 関数で生成されたデータの検定を,rank 関数を用いて行います.

まず,データを用意します(ここでは正規乱数).

[]SATELLITE[]~/home/demo:[1]% x = nrand(4096,1,0,1.0);

これより,平均 0 分散 1.0 の正規分布に従う乱数が 4096 点生成されます.次に,rank 関数の結果を出力する変数を宣言します.

[]SATELLITE[]~/home/demo:[2]% series u, v, w;

これより,series 型の変数 u,v,w が定義されます.ここで,rank 関数を用います.

rank(x,

入力するデータ系列を関数に与えます.

rank(x,u,

Y 軸の出力時系列 : 度数を格納する変数を設定します.

rank(x,u,v,

X 軸の出力時系列 : 小分割されたデータの範囲を格納する変数を設定します.

rank(x,u,v,w,

Y 軸の出力時系列 : フィットした正規分布の密度関数値を格納する変数を設定します.

rank(x,u,v,w,min(x),

階級の最小値を設定します.

rank(x,u,v,w,min(x),max(x),

階級の最大値を設定します.

[]SATELLITE[]~/home/demo:[3]% rank(x,u,v.w,min(x),max(x),100);

階級の個数を設定します.

この関数を実行すると,u,v,w に計算結果が格納され

   ** MEAN VALUE    =    -0.018528
      VARIANCE      =     0.99103
      STANDARD D    =     0.9955
      3-ORDER M     =     0.0070914
      SKEWNESS      =     0.0071879
      4-ORDER M     =     3.0125
      KURTOSIS      =     0.067269
  -0.01852805

という出力が返ってきます.上から,データの平均値,データの分散,データの標準偏差,3 次のモーメント,skewness (歪度),4 次のモーメント,kurtosis (尖度) を表しています.skewness は分布の左右対称性の目安としてよく使われます.これは,零に近いほど平均を基準に分布の左右が対象であることを示しています.また,kurtosis は確率分布がどの程度尖ったものとなるかを示す尺度で,データが正規分布に従う場合には 3 に近くなります.rank 関数の仕様では,主に正規性の検定に用いることを想定しているため,定義式から 3 を減じた値を kurtosis として表示していますので注意して下さい.

次に,計算した結果をグラフに表示する例を示します.

[]SATELLITE[]~/home/demo:[4]% wopen(1,"A4",0,0);
[]SATELLITE[]~/home/demo:[5]% size(80,80);
[]SATELLITE[]~/home/demo:[6]% origin(30,30);
[]SATELLITE[]~/home/demo:[7]% title(1,"Data Point","histogram");
[]SATELLITE[]~/home/demo:[8]% scale("N","F","N","F",min(v),max(v),min(u),max(u));
[]SATELLITE[]~/home/demo:[9]% graph(u,v,0,0,0,1,3.5);
[]SATELLITE[]~/home/demo:[10]% graph(w,v,0,0,0,1,3.5);
[]SATELLITE[]~/home/demo:[11]% axis(1,1,"XY","XY",5,0,0,0,0,0);

これらの関数を実行すると,図 4.18. 「正規乱数のヒストグラムとその正規分布」 のように正規分布およびヒストグラムが表示されます.

正規乱数のヒストグラムとその正規分布

図 4.18. 正規乱数のヒストグラムとその正規分布

Last updated: 2005/03/31