/*********************************************

Transcription

1 デイ ディスプ DADiSP WORKSHEET データ解析 / ディスプレイソフトウェアオプション Unit/C.R.F マニュアル 2010 年 8 月 Copyright (C) 2007 DSP Development Corporation Version 1.6

2 Unit / Custom Response Functions Version 1.6 CUSTOM RESPONSE FUNCTIONS ライセンス同意書... 1 保証および法的責任の放棄... 1 まえがき... 2 第 1 章 CUSTOM RESPONSE FUNCTIONS の基礎 CUSTOM RESPONSE FUNCTIONS とは インストール 使い方... 3 第 2 章 CUSTOM RESPONSE FUNCTIONS の概略 振幅応答 位相応答 コヒーレンス フィルター PSD の平均 サンプリング その他 第 3 章関数リファレンス 振幅応答 メニュー項目 : /* 振幅応答 ( 単純 ) */ メニュー項目 : /* 振幅応答 2 */ メニュー項目 : /* 振幅応答 3 */ メニュー項目 : /* 振幅応答 ( パルス用 ) */ メニュー項目 : /* 振幅応答 (om) */ メニュー項目 : /* 振幅応答 (im,om) */ メニュー項目 : /* 振幅応答 (im,op) */ メニュー項目 : /* 振幅応答 (ip,op) */ MULTI 出力振幅応答 メニュー項目 : /* Mlt 振幅応答 ( 単純 ) */ メニュー項目 : /* Mlt 振幅応答 2 */ メニュー項目 : /* Mlt 振幅応答 3 */ メニュー項目 : /* Mlt 振幅応答 ( パルス用 ) */ メニュー項目 : /* Mlt 振幅応答 (om) */ メニュー項目 : /* Mlt 振幅応答 (im,om) */ メニュー項目 : /* Mlt 振幅応答 (im,o P) */ メニュー項目 : /* Mlt 振幅応答 (ip,op) */ 位相応答 メニュー項目 : /* 位相応答 ( 単純 ) */ メニュー項目 : /* 位相応答 (om) */ メニュー項目 : /* 位相応答 (im,om) */ メニュー項目 : /* 位相応答 (im,op) */ メニュー項目 : /* 位相応答 (ip,op) */ メニュー項目 : /* 位相応答 ( 連続 ) */ MULTI 出力位相応答 メニュー項目 : /* Mlt 位相応答 ( 単純 ) */ メニュー項目 : /* Mlt 位相応答 (om) */ i

3 メニュー項目 : /* Mlt 位相応答 (im,om) */ メニュー項目 : /* Mlt 位相応答 (im,o P) */ メニュー項目 : /* Mlt 位相応答 (ip,op) */ メニュー項目 : /* Mlt 位相応答 ( 連続 ) */ コヒーレンス メニュー項目 : /* コヒーレンス ( 単純 ) */ メニュー項目 : /* コヒーレンス (om) */ メニュー項目 : /* コヒーレンス (im,om) */ メニュー項目 : /* コヒーレンス (im,op) */ メニュー項目 : /* コヒーレンス (ip,op) */ メニュー項目 : /* コヒーレンス (H, 単純 ) */ メニュー項目 : /* コヒーレンス (H,oM) */ メニュー項目 : /* コヒーレンス (H,iM,oM) */ メニュー項目 : /* コヒーレンス (H,iM,oP) */ メニュー項目 : /* コヒーレンス (H,iP,oP) */ MULTI 出力コヒーレンス メニュー項目 : /* Mltコヒーレンス ( 単純 ) */ メニュー項目 : /* Mltコヒーレンス (om) */ メニュー項目 : /* Mltコヒーレンス (im,om) */ メニュー項目 : /* Mltコヒーレンス (im,op) */ メニュー項目 : /* Mltコヒーレンス (ip,op) */ メニュー項目 : /* Mltコヒーレンス (H, 単純 ) */ メニュー項目 : /* Mltコヒーレンス (H,oM) */ メニュー項目 : /* Mltコヒーレンス (H,iM,oM) */ メニュー項目 : /* Mltコヒーレンス (H,iM,oP) */ メニュー項目 : /* Mltコヒーレンス (H,iP,oP) */ フィルター メニュー項目 : /* Filter(0*5) 作成 */ メニュー項目 : /* Filter(0*5) 適用 */ メニュー項目 : /* Filter(0*3) 作成 */ メニュー項目 : /* Filter(0*3) 適用 */ PSDの平均 (X(0)=0) メニュー項目 : /* PSDの平均 ( 窓なし ) */ メニュー項目 : /* PSDの平均 (hanning) */ メニュー項目 : /* PSDの平均 (hamming) */ メニュー項目 : /* PSDの平均 (kaiser) */ メニュー項目 : /* PSDの平均 ( 各列 )( 窓なし ) */ メニュー項目 : /* PSDの平均 ( 各列 )(hanning) */ メニュー項目 : /* PSDの平均 ( 各列 )(hamming) */ メニュー項目 : /* PSDの平均 ( 各列 )(kaiser) */ PSDの平均 (X(0)!=0) メニュー項目 : /* PSDの平均 (X0!=0) */ メニュー項目 : /* PSDの平均 ( 各列 ) (X0!=0) */ サンプリング メニュー項目 : /* サンプリング */ メニュー項目 : /* サンプリング ( 番号付 ) */ その他 メニュー項目 : /* 列データの拡大 シフト */ メニュー項目 : /* 前回平均値比較機能 */ メニュー項目 : /* 相互相関関数 ( 正規化 ) */ メニュー項目 : /* 移動積分 */ メニュー項目 : /* 振幅 /(2πf) */ メニュー項目 : /* 心拍数と血圧 */ ii

4 Custom Response Functions ライセンス同意書 ディスク入りの封筒を開ける前に重要な Custom Response Functions ライセンス同意書を読んで 内容を良く理解して下さい この封筒を開封することで以下の諸条件についてあなたの同意が得られたものとみなされます 本ライセンス受権者は提供するソフトウエアに対する非独占的使用権を与えられます 提供するプログラムを 本書の第 1 章に記載の インストール で述べる方法以外でコピーすることを禁止します 本ライセンス受権者が本プログラムまたはその文書のいかなる部分をもライセンスを保持しない第三者に公開したり 再販売することを禁止します 保証および法的責任の放棄本文書の内容には信頼性があり また正確であるものと確信しております しかし著者宇山靖政ならびにアストロデザイン株式会社ではいかなる誤り 抜け落ち また不正確さに対してその責任を負うことはありません 著者宇山靖政ならびにアストロデザイン株式会社は 本マニュアルの内容およびそこに記述するソフトウエアに関して特定の機能に対する市場性および適合性の保証を始めとして いかなる またあらゆる明示的または暗黙の保証を放棄します 本製品の品質および性能に関する危険 ( リスク ) は本製品の購入者またはユーザーに帰属することになります 著者宇山靖政およびアストロデザイン株式会社が事前にその可能性について通知を受けている場合においても かような情報およびソフトウエアの使用に起因するような特殊なまたは結果的な損害を始めとして いかなる損害に対しても法的責任を負うことはありません 1

5 まえがき これは DADiSP で提供される標準の関数以外のもので 過去に DADiSP の利用者からの要望に答えて作成した応答関数を中心に特別仕様の関数を集めたものです ( 特別仕様とはいえ あまりにも特殊なものは除きました ) 作成の方針は 必要なものだけを表示する です コンピュータでは計算誤差を避けることは不可能ですが 必要なものだけを表示することは可能です この方針で作成された関数ですので 他の関数とは違って特にシャープな表現力を持たせてあります 関数は DADiSP の SPL によって作成されていますので C 言語を理解できる方なら SPL File をご覧になればアルゴリズムなどはご理解いただけるものと思います とはいえ 自分で作るにはかなりの時間がかかります 特殊だけれど便利な関数を集めたものです 多くの方に利用していただきたいと思います また 新たな関数についてのご要望がありましたら作者又はアストロデザインまでお知らせ下さい この作品をよりよいものにしてゆきたいと思っておりますので 内容の不明な点 ご指摘などありましたら 著者の アドレス uyama@kc4.so-net.ne.jp へお知らせ お問い合わせをしていただければ幸いです 2001 年 11 月 宇山靖政 Version 1.6 について古い SPL ファイルが DADiSP2002 ではうまく動かないことが分かりました うまく動くように修正しました 2007 年 11 月 4 日 宇山靖政 2

6 第 1 章 Custom Response Functions の基礎 1.1 Custom Response Functions とは DADiSP には ユーザーによるカスタマイズを可能とする機能があります C 言語に似た SPL (Series Processing Language) を利用して必要な機能を DADiSP に付け加えることができます Custom Response Functions は SPL を利用して 応答関数を中心とする特別仕様の関数を提供するものです 関数の利用は メニュー方式で必要な機能を選択すれば データの解析ができるようになっています 提供されるファイルはテキストファイルですので もし必要ならば自分の手で改造することもできます この Custom Response Functions(Ver.1.6) を利用するには DADiSP2002 以上が必要です 1.2 インストール CD は CrfJ フォルダ 2.Unit CRFJ(DADiSP ラボブック ) 2.Dadisp.mac 3.Dadisp.spl 上のようになっています CrfJ フォルダを DADiSP ディレクトリ下にある SPL のフォルダにコピーをして下さい 階層的には C: program files\dsp2002 Spl CrfJ となります まだ Dadisp.mac と Dadisp.spl と Unit CRFJ を使用していない方は DADiSP のあるディレクトリ C: program files\dsp2002 にコピーをしてください Crf 用のツールボタンが DADiSP 起動時に自動的に作成されます また Unit CRFJ は ファイル 開く ラボブックで Unit CRFJ を選択してそれぞれのワークシートを選択し 使用して下さい * すでに Dadisp.mac Dadisp.spl をご使用の方はコピーしないで以下の指示に従って下さい Dadisp.spl は すでにウエーブレット 多変量解析を使っているものとして 次のようになっていますので 他の場合には splmain() ( toolbar(1,-1,4,red,"mva",'menufile("spl\mva\mvmenu.men")'," 多変量解析 "); toolbar(1,-1,4,green,"wlt",'menufile("spl\wltj\waveletj.men")',"wavelet"); toolbar(1,-1,4,green," 応答関数 ",'MENUFILE("spl\crfj\crfj.men")'," 応答関数 "); ) この最後の一行をご自分のファイルに追加して下さい Dadisp.mac は #include macros\matrix.mac となっているだけですので この行を書き加えて下さい 1.3 使い方 1. ファイル - 新規作成 - ラボブックを選択し フォルダ \dsp2002 の中に UnitCRFJ を作成します このとき 新規ワークシートが作成されます 2. 応答関数 - 関数読込みとして 解析に必要な関数を読み込みます 関数は DADiSP を起動してから一度だけ読み込めば十分です 3. ワークシートのウインドウ (W1 とします ) にデータを置きます 3

7 4. 別のウインドウに移動しデータが入っているウインドウに関数を作用させます 関数の選択は ツールバー応答関数のクリック 関数の選択 関数のパネルでデータウインドウなどの引数の指定をします 5. 関数パネルの入力方法 上の図は心電図と血圧の例ですが 心電図データ : 心電図のデータが入っているウインドウを w1 の様にして指定します 基準値 : 実数値を 0.6 のように入力します Up-down : 整数値を 1 のように入力します 心拍数 : 心拍数が表示されるウインドウ番号の数値だけを 3 のように入力します 血圧データ : 血圧のデータが入っているウインドウを w2 のように指定します 最大値 : 血圧の最大値を表示するウインドウの番号の数値だけを 4 のように入力します 最大値 : 血圧の最小値を表示するウインドウの番号の数値だけを 5 のように入力します 最大値 : 血圧の平均値を表示するウインドウの番号の数値だけを 6 のように入力します 注意 : 計算対象のデータが入っているウインドウの指定には W3 のように W をつけ 4

8 計算結果を出力するウインドウの指定には 数値のみを入力します 第 2 章 Custom Response Functions の概略 2.1 振幅応答 振幅応答では 単純にスペクトラムの比を取るものから ピーク値を取る周波数成分についてだけ計算するものまで 8 種類を用意しました 特徴は次のようになります 1. 振幅応答 ( 単純 ) 各周波数について単純にスペクトラムの比をとったものです 2. 振幅応答 2 入力波形 出力波形のスペクトラムをとり それぞれの移動平均を取ってから 各周波数についての比をとったものです 3. 振幅応答 3 入力波形 出力波形のスペクトラムをとり 各周波数についての比をとってから 各周波数についての移動平均をとったものです 4. 振幅応答 ( ハ ルス用 ) 入力波形を完全なパルスに変換してから 入力波形 出力波形のスペクトラムを取り比をもとめたもの 必要ならその後移動平均を取ることも出来ます 5. 振幅応答 (om) 入力波形 出力波形のスペクトラムをとり 各周波数についての比をとったもの 必要ならその後移動平均を取ることも出来ます ただし 出力波形のスペクトラムを取ったときに その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱っています 6. 振幅応答 (im,om) 入力波形 出力波形のスペクトラムをとり 各周波数についての比をとったもの 必要ならその後移動平均を取ることも出来ます ただし 入力波形のスペクトラムを取ったときに ( 最大値 * 入力 0 化 ) 以上の強さを持ち しかも出力波形のスペクトラムを取ったときに ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱っています 7. 振幅応答 (im,op) 入力波形 出力波形のスペクトラムをとり 各周波数についての比をとったもの 必要ならその後移動平均を取ることも出来ます ただし 入力波形のスペクトラムを取ったときに ( 最大値 * 入力 0 化 ) 以上の強さを持ち さらに出力波形のスペクトラムを取ったときにピーク値となり ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱っています 8. 振幅応答 (ip,op) 入力波形 出力波形のスペクトラムをとり 各周波数についての比をとったもの 必要ならその後移動平均を取ることも出来ます ただし 入力波形のスペクトラムを取ったときにピーク値となりしかも ( 最大値 * 入力 0 化 ) 5

9 以上の強さを持ち さらに出力波形のスペクトラムを取ったときにピーク値となり ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱っています 入力データと出力データに関して 各周波数成分について それぞれの強さの比を求めるときに 扱う信号の性質によって 関数を使い分けると効果的な場合があります 最初に 入力がパルスに近い場合を考えてみます 上の図の w1 がパルスに近い性質を持つ入力 w2 が出力です W1: パルス状の入力 w2: 出力波形 w3: 振幅応答 ( 単純 ) w4: 振幅応答 2 w5: 振幅応答 3 w6: 振幅応答 ( パルス ) です 各振幅応答では 250Hz 付近のピークは正確ですが 高周波成分に関しては w3,w4,w5 では かなりの誤差が見られます 入力を純粋なパルスに変形してから計算した w6 では きれいなカーブが得られています パルスのような入力に対しては この 振幅応答 ( パルス ) で計算するのが良いでしょう 6

10 下の図のような場合は 入力 (w2) 出力 (w6) ともに 1Hz 3Hz 9Hz の周波数成分が主な成分です w1:gline(2^7,1/2^7,2*pi,0) w2:sin(w1)-sin(3*w1)+sin(9*w1) w3:gline(2^6,1,0) w4:gline(2^6,1/2^7,2*pi,0) w5:sin(9*w1) w6:sin(w1)-2*sin(3*w1)+concat(w3,w5) それに対する w7:magresp5(w2,w6,0.01,1) w8:magresp7(w2,w6,0.001,0.1,1) w9:magresp8(w2,w6,0.01,0.001,1) w10:spectrum(w2) w11:spectrum(w6) 振幅応答 (om) 振幅応答 (im,op) 振幅応答 (ip,op) の結果です w7 は不適当な値です この場合は 振幅応答を調べるときに 出力波形の specutrum を取ったときに現れる値のうちで小 7

11 さなものを計算誤差と見て完全な 0 に変えてから計算すべきでしょう このようにして計算したものが w8 w9 です 特定の周波数成分についての計算結果がきれいに表現されています このことから 入力データの性質 出力データの性質を考えて どの振幅応答関数を使うかを決定しなくてはならないことが分かります これは w10,w11 のようなスペクトラムの結果を見てから決定するべきであると言うことを示しています 入力が 1 チャンネルでそれに対する出力が多チャンネルの場合は次の関数をご利用下さい 9.Mlt 振幅応答 ( 単純 ) 10.Mlt 振幅応答 2 11.Mlt 振幅応答 3 12.Mlt 振幅応答 ( ハ ルス用 ) 13.Mlt 振幅応答 (om) 14.Mlt 振幅応答 (im,om) 15.Mlt 振幅応答 (im,op) 16.Mlt 振幅応答 (ip,op) これらは 1.~8. で示された関数と同じアルゴリズムを 入力データと出力の各列データに対して繰り返して適用したものです 8

12 2.2 位相応答 位相応答では 入力 出力に関してそれぞれの FFT() を取った後 各周波数成分ごとの商を求めその偏角を計算しているのですが 入出力の特徴に合わせていくつかの関数を用意しました 特徴は次のようになります 1. 位相応答 ( 単純 ) これは 入力 出力に関してそれぞれの FFT() を取った後 各周波数成分ごとの商を求めその偏角を計算しています 2. 位相応答 (om) これは 入力 出力に関してそれぞれの FFT() を取った後 各周波数成分ごとの商を求めその偏角を計算してから 出力のスペクトラムをとり その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います 3. 位相応答 (im,om) これは 入力 出力に関してそれぞれの FFT() を取った後 各周波数成分ごとの商を求めその偏角を計算してから 出力のスペクトラムをとり その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います さらに 入力波形のスペクトラムをとり その ( 最大値 * 入力 0 化 ) 以上の強さを持つ周波数成分以外も 0 として扱っています 4. 位相応答 (im,op) これは 入力 出力に関してそれぞれの FFT() を取った後 各周波数成分ごとの商を求めその偏角を計算してから 出力のスペクトラムをとり スペクトラムのピーク値となっていて しかもその ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います さらに 入力波形のスペクトラムをとり その ( 最大値 * 入力 0 化 ) 以上の強さを持つ周波数成分以外も 0 として扱っています 5. 位相応答 (ip,op) これは 入力 出力に関してそれぞれの FFT() を取った後 各周波数成分ごとの商を求めその偏角を計算してから 出力のスペクトラムをとり スペクトラムのピーク値となっていてしかも その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います さらに 入力波形のスペクトラムをとり スペクトラムのピーク値となっていて しかもその ( 最大値 * 入力 0 化 ) 以上の強さを持つ周波数成分以外も 0 として扱っています 6. 位相応答 ( 連続 ) これは 入力 出力に関してそれぞれの FFT() を取った後 各周波数成分ごとの商を求めその偏角を計算しています この偏角の変化が微小な周波数の変化に対して極端に大きく (2π に近い ) 場合はこれを連続化して扱っています 基本的な操作としては 入力波形と出力波形の位相差を FFT から計算するのですが 単純に計算す 9

13 ると計算誤差のためにあまりよい結果が得られません この場合は 出力の周波数成分のうちで 小さいもの ( 出力波の spectrum のうちで最大値に比べて 倍以下のもの ) を 完全な 0 に変えてから計算します これが 下の図の内容です w1:gline(2^10,1/2^10,2*pi,0) w2:cos(w1*4)-sin(35*w1)+sin(200*w1) w3:exp(-w1)*(2*sin(4*w1)-3*cos(35*w1)+4*cos(200*w1)) w4:phaseresp(w2,w3) 位相応答 ( 単純 ) w5:phaseresp2(w2,w3,0.1) 位相応答 (om) w6:phaseresp3(w2,w3,0.01,0.001) 位相応答 (im,om) w7:phaseresp4(w2,w3,0.0,0.001) 位相応答 (im,op) w8:phaseresp5(w2,w3,0.0,0.001) 位相応答 (ip,op) w9:cphaseresp(w2,w3,5.0) 位相応答 ( 連続 ) w10:spectrum(w3) w4 では w3 と w2 の位相差がはっきり現れています 入力出力ともに周波数成分を連続的に含んでいる場合は w5,w6 で使用した関数がよいと思います 扱う波形の周波数成分をそれぞれ調べてから どの関数を利用するかを決めなくてななりません 10

14 1 チャンネルの入力データに対して多チャンネルの出力データがある場合は次の関数をご利用下さい 7.Mlt 位相応答 ( 単純 ) 8.Mlt 位相応答 (om) 9.Mlt 位相応答 (im,om) 10. Mlt 位相応答 (im,op) 11. Mlt 位相応答 (ip,op) 12.Mlt 位相応答 ( 連続 ) これらは 1.~6. で示された関数と同じアルゴリズムを 入力データと出力の各列データに対して繰り返して適用したものです 11

15 2.3 コヒーレンス コヒーレンスの計算でも 計算誤差が問題になります そこで 次のような関数を用意いたしました 特徴は次のようになります 1. コヒーレンス ( 単純 ) 入力信号 出力信号の長さは同じにして下さい 指定された長さで 入出力信号を切り出します このとき 指定された長さの半分だけ重なるようにして切り出します 切り出された信号から コヒーレンスを計算します 2. コヒーレンス (om) 入力信号 出力信号の長さは同じにして下さい 指定された長さで 入出力信号を切り出します このとき 指定された長さの半分だけ重なるようにして切り出します 切り出された信号から コヒーレンスを計算します ただし 出力信号のスペクトラムを取って その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います 3. コヒーレンス (im,om) 入力信号 出力信号の長さは同じにして下さい 指定された長さで 入出力信号を切り出します このとき 指定された長さの半分だけ重なるようにして切り出します 切り出された信号から コヒーレンスを計算します ただし 出力信号のスペクトラムを取って その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います さらに 入力信号のスペクトラムをとって その ( 最大値 * 入力 0 化 ) 以上の強さを持つ周波数成分だけを計算の対象にします それ以外は 0 として扱います 4. コヒーレンス (im,op) 入力信号 出力信号の長さは同じにして下さい 指定された長さで 入出力信号を切り出します このとき 指定された長さの半分だけ重なるようにして切り出します 切り出された信号から コヒーレンスを計算します ただし 出力信号のスペクトラムを取り さらにピーク値となっている周波数成分のうちで その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います さらに 入力信号のスペクトラムをとって その ( 最大値 * 入力 0 化 ) 以上の強さを持つ周波数成分だけを計算の対象にします それ以外は 0 として扱います 5. コヒーレンス (ip,op) 入力信号 出力信号の長さは同じにして下さい 指定された長さで 入出力信号を切り出します このとき 指定された長さの半分だけ重なるようにして切り出します 切り出された信号から コヒーレンスを計算します ただし 出力信号のスペクトラムを取り さらにピーク値となっている周波数成分のうちで その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います さらに 入力信号のスペクトラムをとり さらに ピーク値となっている周波数成分で その ( 最大値 * 入力 0 化 ) 以上の強さを持つ周波数成分だけを計算の対象にします それ以外は 0 として扱います 6. コヒーレンス (H, 単純 ) 入力信号 出力信号の長さの半分の長さの信号を切り出します このとき 長さの半分だけ重なるようにして切り出します 切り出された信号から コヒーレンスを計算します 12

16 7. コヒーレンス (H,oM) 入力信号 出力信号の長さの半分の長さの信号を切り出します このとき 長さの半分だけ重なるようにして切り出します 切り出された信号から コヒーレンスを計算します ただし 出力信号のスペクトラムを取って その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います 8. コヒーレンス (H,iM,oM) 入力信号 出力信号の長さの半分の長さの信号を切り出します このとき 長さの半分だけ重なるようにして切り出します 切り出された信号から コヒーレンスを計算します ただし 出力信号のスペクトラムを取って その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います さらに 入力信号のスペクトラムをとって その ( 最大値 * 入力 0 化 ) 以上の強さを持つ周波数成分だけを計算の対象にします それ以外は 0 として扱います 9. コヒーレンス (H,iM,oP) 入力信号 出力信号の長さの半分の長さの信号を切り出します このとき 長さの半分だけ重なるようにして切り出します 切り出された信号から コヒーレンスを計算します ただし 出力信号のスペクトラムを取り さらにピーク値となっている周波数成分のうちで その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います さらに 入力信号のスペクトラムをとって その ( 最大値 * 入力 0 化 ) 以上の強さを持つ周波数成分だけを計算の対象にします それ以外は 0 として扱います 10. コヒーレンス (H,iP,oP) 入力信号 出力信号の長さの半分の長さの信号を切り出します このとき 長さの半分だけ重なるようにして切り出します 切り出された信号から コヒーレンスを計算します ただし 出力信号のスペクトラムを取り さらにピーク値となっている周波数成分のうちで その ( 最大値 * 出力 0 化 ) 以上の強さを持つ周波数成分以外は 0 として扱います さらに 入力信号のスペクトラムをとり さらに ピーク値となっている周波数成分で その ( 最大値 * 入力 0 化 ) 以上の強さを持つ周波数成分だけを計算の対象にします それ以外は 0 として扱います 13

17 コヒーレンスの計算ですが ここでも計算誤差が問題になります w1:gline(2^10,1/2^10,2*pi,0) w2:sin(w1*4)-sin(35*w1)+sin(200*w1) w3:2*cos(w1*4)-3*cos(35*w1)+4*cos(200*w1) w2 と w3 の波に対してコヒーレンスを計算しています w4:ucohh(w2,w3) コヒーレンス (H, 単純 ) w5:ucoh3h(w2,w3,0.0,0.001) コヒーレンス (H,iM,oM) w6:ucoh4h(w2,w3,0.0,0.001) コヒーレンス (H,iM,oP) w7:ucoh5h(w2,w3,0.0,0.001) コヒーレンス (H,iP,oP) w8:spectrum(w3,256) となっていますが 特定の周波数成分しか持っていない場合は w6,w7 の方が w4 w5 よりも正確なものとなっています 原因は 切り出す長さが短く w8 にあるようにスペクトラムの計算結果がシャープではないからです 注意 : 周波数軸上の目盛りは もとの信号に対して単純に FFT や spectrum を取ったときよりも大きな幅になります これは信号を短い長さに切り出しているからです 14

18 1 チャンネルの入力に対して 多チャンネルの出力データが得られる場合は次の関数をご利用下さい 11.Mlt コヒーレンス ( 単純 ) 12.Mlt コヒーレンス (om) 13.Mlt コヒーレンス (im,om) 14.Mlt コヒーレンス (im,op) 15.Mlt コヒーレンス (ip,op) 16.Mlt コヒーレンス (H, 単純 ) 17.Mlt コヒーレンス (H,oM) 18.Mlt コヒーレンス (H,iM,oM) 19.Mlt コヒーレンス (H,iM,oP) 20.Mlt コヒーレンス (H,iP,oP) これらは 1.~10. で示された関数と同じアルゴリズムを 入力データと出力の各列データに対して繰り返して適用したものです 15

19 2.4 フィルター 電源からの雑音は 50Hz(60Hz) およびその高調波からなっています この雑音を取り除くためのフィルタとして作成したのですが 少し拡張して任意の周波数およびその高調波を取り除くようにしてあります 特徴はつぎのようになります 1. Filter(0*5) 作成入力信号の FFT() をとります これに対してフィルタを作成します 0 が 5 つずつ並んだ構造になっています 2. Filter(0*5) 適用入力信号から特定の周波数およびその高調波を取り除きます このときのフィルタは 0 が 5 つ並んでいます 3. Filter(0*3) 作成入力信号の FFT() をとります これに対してフィルタを作成します 0 が 3 つずつ並んだ構造になっています これによって 他の周波数成分に関する影響が少なくなっています 4.Filter(0*3) 適用入力信号から特定の周波数およびその高調波を取り除きます このときのフィルタは 0 が 3 つ並んでいます これによって 他の周波数成分に関する影響が少なくなっています 交流電源は 50Hz または 60Hz です この電源による雑音を除去するためのフィルターです 16

20 与えられた周波数とその自然数倍の周波数を持つ波を除去します 下の図では 1Hz 50Hz 60Hz の波を加えたものが w2 です この波から 50Hz 成分を取り除いたものが W5 W6 です w1:gline(1024,1/1024,2*pi,0) w2:sin(w1)-sin(50*w1)+sin(60*w1) w3:fft(w2) w4:filter(0*5) 作成 w5:filter(0*5) 適用 w6:ifft(w3*w4) w7:sin(w1)+sin(60*w1) w8:w5-w7 w8 から 50Hz の成分だけがきれいに除去されていることが分かります 注意 : 周波数の指定を変えれば 他の周波数の波を取り除くことができます ただし その自然数倍の周波数を持つ波も除去しますので注意してください 17

21 2.5 PSD の平均 PSD の平均 (X(0)=0) これは PSD の計算に際して 平均値を引いてから計算を始めます 従って 直流成分が取り除かれることになります 1 チャンネルデータに対しての関数です PSD の平均 ( 窓なし ) PSD の平均 (hanning) PSD の平均 (hamming) PSD の平均 (kaiser) 多チャンネルデータに対して 各列ごと n それぞれの平均値を引いてから計算を始めます PSD の平均 ( 各列 )( 窓なし ) PSD の平均 ( 各列 )(hanning) PSD の平均 ( 各列 )(hamming) PSD の平均 ( 各列 )(kaiser) PSD の平均 (X(0)!=0) これは 平均値を引かないで計算をします 従って直流成分が残ります 1チャンネルデータに対しての関数です PSD の平均 (X(0)!=0) 多チャンネルデータに対しての関数です PSD の平均 ( 各列,X(0)!=0) これらの計算では 入力データを幾つかに分割して それぞれの PSD を取ってから平均をとっています 2.6 サンプリング サンプリングサンフ リンク ( 番号付 ) これは 与えられたデータからランダムにサンプリングするものです その結果は最初に並んでいた順番を保ちながら返ってきます サンプリングする時に利用する乱数は種数によって決まる疑似乱数です 乱数の発生は最初のデータ数と同じ数だけですので サンプリング数が最初のデータ数の半分以上になると失敗する確率が高まります 失敗したときは 0 の列が返ります 2.7 その他 列データの ( 拡大 シフト ) 前回平均値比較機能相互相関 ( 正規化 ) 移動積分 これについては それぞれの関数の説明をご覧下さい 18

22 第 3 章 関数リファレンス 3.1 振幅応答 メニュー項目 : /* 振幅応答 ( 単純 ) */ 機能 : 入力と出力の周波数成分ごとの強さの比を求めます 書式 :magresp(indat, outdat) indat: 入力データ outdat: 出力データ ( データサイズは 2^k) ( データサイズは 2^k) 戻り値 : 周波数成分ごとの強さの比 例 : w1:gline(2^10,1/1024,2*pi,0) w2:10*sin(w1*5)-5*sin(30*w1)+25*sin(95*w1) w3:w2*exp(-w1) w4:concat(gser(0,0,0,0,0,0,0,0,50,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0),gline(1000,1/1024,0,0));setdeltax(1/1024) w5:spectrum(w3) w6: 振幅応答 ( 単純 ) 説明 : 単純に spectrum の比を取っているので 計算誤差などの問題を含んでいます これは 入力が完全なパルスなのできれいに表示されますが w4 の信号が変化するときれいに表示されなくなります メニュー項目 : /* 振幅応答 2 */ 機能 : 入力と出力の周波数成分ごとの強さの比を求めます 書式 :magresp2(indat, outdat, wid) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) wid: 移動平均の幅 戻り値 : 周波数成分ごとの波の強さの比 説明 : 入力 出力の spectrum をとり それぞれに対して移動平均を取ってから比を求めています 本来 入力 出力ともに周波数成分の強さが連続的に変化すると思える場合に使用して下さい 19

23 メニュー項目 : /* 振幅応答 3 */ 機能 : 入力と出力の周波数成分ごとの強さの比を求めます 書式 :magresp3(indat, outdat, wid) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) wid: 移動平均の幅 戻り値 : 周波数成分ごとの波の強さの比 説明 : 入力 出力の spectrum をとり その比を取ってから移動平均を取っています 本来 振幅応答が連続的であると思える場合に使用して下さい メニュー項目 : /* 振幅応答 ( パルス用 ) */ 機能 : 入力 ( パルス ) と出力の周波数成分ごとの強さの比を求めます 書式 :magresp4(indat, outdat, wid) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) wid: 移動平均の幅 戻り値 : 周波数成分ごとの波の強さの比 説明 : 最初に 入力波形を理想的なパルスに変えます そして 理想化された入力と出力の spectrum をとり その比を取ってから移動平均を取っています 入力が パルス状態の信号である場合に使用して下さい メニュー項目 : /* 振幅応答 (om) */ 機能 : 入力と出力の周波数成分ごとの強さの比を求めます 書式 :magresp5(indat, outdat, outc, wid) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) outc: 出力のスペクトラムで最大値 *outc 以下の大きさの周波数成分は完全に0とする wid: 最後にとる移動平均の幅 戻り値 :0 にならなかった周波数成分についての 周波数ごとの波の強さの比 20

24 メニュー項目 : /* 振幅応答 (im,om) */ 機能 : 入力と出力の周波数成分ごとの強さの比を求めます 書式 :magresp6(indat, outdat, inc, outc, wid) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) inc: 入力のスペクトラムで最大値 *inc 以下の大きさの周波数成分は完全に0とする outc: 出力のスペクトラムで最大値 *outc 以下の大きさの周波数成分は完全に0とする wid: 最後にとる移動平均の幅 戻り値 :0 にならなかった周波数成分についての 周波数ごとの波の強さの比 メニュー項目 : /* 振幅応答 (im,op) */ 機能 : 入力と出力の周波数成分ごとの強さの比を求めます 書式 :magresp7(indat, outdat, inc, outc, wid) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) inc: 入力のスペクトラムで ( 最大値 *inc) 未満の大きさの周波数成分は完全に0とする outc: 出力のスペクトラムについて ピーク値となっている周波数成分のうちで ( 最大値 *outc) 未満の大きさの周波数成分は完全に0とする wid: 最後にとる移動平均の幅 戻り値 :0 にならなかった周波数成分についての 周波数ごとの波の強さの比 メニュー項目 : /* 振幅応答 (ip,op) */ 機能 : 入力と出力の周波数成分ごとの強さの比を求めます 書式 :magresp8(indat, outdat, inc, outc, wid) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) inc: 入力のスペクトラムについて ピーク値となっている周波数成分のうちで ( 最大値 *inc) 以上の大きさの周波数成分以外は完全に0とする outc: 出力のスペクトラムについて ピーク値となっている周波数成分のうちで ( 最大値 *outc) 以上の大きさの周波数成分以外は完全に0とする wid: 最後にとる移動平均の幅 戻り値 :0 にならなかった周波数成分についての 周波数ごとの波の強さの比 21

25 例 : w1:gline(2^10,1,5*2*pi/2^10,0) w2:sin(w1)-sin(3*w1)+sin(9*w1) w3:2*cos(w1)-3*cos(3*w1)+4*cos(9*w1) w4:magresp(w2,w3) 振幅応答 ( 単純 ) w5:spectrum(w3) w6:magresp5(w2,w3,0.0,0.001,1) 振幅応答 (im,om) 説明 :w6 では 入力 出力の spectrum をとります その後で その比を取ってから移動平均を取っています 入力に関しては スペクトラムの ( 最大値 *inc) よりも小さな値は完全に 0 に置き換え 弱い周波数成分は 0 として扱っています 出力に関しては スペクトラムの ( 最大値 *outc) よりも小さな値は完全に 0 に置き換え 弱い周波数成分は 0 として扱っています 周波数成分が少なく 本来 0 になるべき値が計算誤差で 0 にならず そのような誤差と誤差の割り算の結果として 振幅応答の計算に問題が生じる場合に使用して下さい これに対して w4 では計算誤差の影響がかなり強く出ています w5 の spectrum() に現れた周波数成分だけについて表示すればよいのです 22

26 3.2 Multi 出力振幅応答 メニュー項目 : /* Mlt 振幅応答 ( 単純 ) */ 機能 :1 チャンネル入力と多チャンネル出力の周波数成分ごとの強さの比を求めます 書式 :mmagresp(indat, outdat) indat: 入力データ (1-ch) outdat: 出力データ (k-ch) ( データサイズは 2^k) ( 各チャンネルのデータサイズは 2^k) 戻り値 : 周波数成分ごとの強さの比 例 : w1:gline(2^10,1/1024,2*pi,0) w2:10*sin(w1*5)-5*sin(30*w1)+25*sin(95*w1) w3:w2*exp(-w1) w4:concat(gser(0,0,0,0,0,0,0,0,50,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0),gline(1000,1/1024,0,0));setdeltax(1/1024) w5:spectrum(w3) w6:mlt 振幅応答 ( 単純 ) 説明 : 単純に spectrum の比を取っているので 計算誤差などの問題を含んでいます これは 入力が完全なパルスなのできれいに表示されますが w4 の信号が変化するときれいに表示されなくなります メニュー項目 : /* Mlt 振幅応答 2 */ 機能 :1 チャンネル入力と多チャンネル出力の周波数成分ごとの強さの比を求めます 書式 :mmagresp2(indat, outdat, wid) indat: 入力データ (1-ch) outdat: 出力データ (k-ch) wid: 移動平均の幅 ( データサイズは 2^k) ( 各チャンネルのデータサイズは 2^k) 戻り値 : 周波数成分ごとの波の強さの比 説明 : 入力 出力の spectrum をとり それぞれに対して移動平均を取ってから比を求めています 本来 入力 出力ともに周波数成分の強さが連続的に変化すると思える場合に使用して下さい 23

27 メニュー項目 : /* Mlt 振幅応答 3 */ 機能 :1 チャンネル入力と多チャンネル出力の周波数成分ごとの強さの比を求めます 書式 :mmagresp3(indat, outdat, wid) indat: 入力データ (1-ch) outdat: 出力データ (k-ch) wid: 移動平均の幅 ( データサイズは 2^k) ( 各チャンネルのデータサイズは 2^k) 戻り値 : 周波数成分ごとの波の強さの比 説明 : 入力 出力の spectrum をとり その比を取ってから移動平均を取っています 本来 振幅応答が連続的であると思える場合に使用して下さい メニュー項目 : /* Mlt 振幅応答 ( パルス用 ) */ 機能 : パルス状の 1 チャンネル入力と 多チャンネル出力の周波数成分ごとの強さの比を求めます 書式 :mmagresp4(indat, outdat, wid) indat: 入力データ (1-ch) outdat: 出力データ (k-ch) wid: 移動平均の幅 ( データサイズは 2^k) ( 各チャンネルのデータサイズは 2^k) 戻り値 : 周波数成分ごとの波の強さの比 説明 : 最初に 入力波形を理想的なパルスに変えます そして 理想化された入力と出力の spectrum をとり その比を取ってから移動平均を取っています 入力が パルス状態の信号である場合に使用して下さい メニュー項目 : /* Mlt 振幅応答 (om) */ 機能 :1 チャンネル入力と多チャンネル出力の周波数成分ごとの強さの比を求めます 書式 :mmagresp5(indat, outdat, outc, wid) indat: 入力データ (1-ch) ( データサイズは 2^k) outdat: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) outc: 出力のスペクトラムで最大値 *outc 以下の大きさの周波数成分は完全に0とする wid: 最後にとる移動平均の幅 戻り値 :0 にならなかった周波数成分についての 周波数ごとの波の強さの比 24

28 メニュー項目 : /* Mlt 振幅応答 (im,om) */ 機能 :1 チャンネル入力と多チャンネル出力の周波数成分ごとの強さの比を求めます 書式 :mmagresp6(indat, outdat, inc, outc, wid) indat: 入力データ (1-ch) ( データサイズは 2^k) outdat: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) inc: 入力のスペクトラムで最大値 *inc 以下の大きさの周波数成分は完全に0とする outc: 出力のスペクトラムで最大値 *outc 以下の大きさの周波数成分は完全に0とする wid: 最後にとる移動平均の幅 戻り値 :0 にならなかった周波数成分についての 周波数ごとの波の強さの比 メニュー項目 : /* Mlt 振幅応答 (im,op) */ 機能 :1 チャンネル入力と多チャンネル出力の周波数成分ごとの強さの比を求めます 書式 :mmagresp7(indat, outdat, inc, outc, wid) indat: 入力データ (1-ch) ( データサイズは 2^k) outdat: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) inc: 入力のスペクトラムで ( 最大値 *inc) 未満の大きさの周波数成分は完全に0とする outc: 出力のスペクトラムについて ピーク値となっている周波数成分のうちで ( 最大値 *outc) 未満の大きさの周波数成分は完全に0とする wid: 最後にとる移動平均の幅 戻り値 :0 にならなかった周波数成分についての 周波数ごとの波の強さの比 メニュー項目 : /* Mlt 振幅応答 (ip,op) */ 機能 :1 チャンネル入力と多チャンネル出力の周波数成分ごとの強さの比を求めます 書式 :mmagresp8(indat, outdat, inc, outc, wid) indat: 入力データ (1-ch) ( データサイズは 2^k) outdat: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) inc: 入力のスペクトラムについて ピーク値となっている周波数成分のうちで ( 最大値 *inc) 以上の大きさの周波数成分以外は完全に0とする outc: 出力のスペクトラムについて ピーク値となっている周波数成分のうちで ( 最大値 *outc) 以上の大きさの周波数成分以外は完全に0とする wid: 最後にとる移動平均の幅 戻り値 :0 にならなかった周波数成分についての 周波数ごとの波の強さの比 25

29 3.3 位相応答 メニュー項目 : /* 位相応答 ( 単純 ) */ 機能 : 入力と出力の周波数成分ごとの位相差を求めます 書式 :phaseresp(indat,outdat) indat: 入力データ outdat: 出力データ ( データサイズは 2^k) ( データサイズは 2^k) 戻り値 : 周波数成分ごとの位相の差からなる数列 説明 : 位相差は -π から +π 間での間で与えられます メニュー項目 : /* 位相応答 (om) */ 機能 : 入力と出力の周波数成分ごとの位相差を求めます 書式 :phaseresp2(indat,outdat,outc) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) outc:0にする出力成分の 最大値に対する比率 戻り値 : 周波数成分ごとの位相の差からなる数列ただし 出力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *outc) 以上の大きさを持つ周波数成分以外は 0 としています メニュー項目 : /* 位相応答 (im,om) */ 機能 : 入力と出力の周波数成分ごとの位相差を求めます 書式 :phaseresp3(indat,outdat,inc,outc) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) inc:0にする入力成分の 最大値に対する比率 outc:0にする出力成分の 最大値に対する比率 戻り値 : 周波数成分ごとの位相の差からなる数列 ただし 出力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *outc) 以上の大きさを持つ周波数成分以外は 0 としています さらに入力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *inc) 以上の大きさを持つ周波数成分以外は 0 としています 26

30 メニュー項目 : /* 位相応答 (im,op) */ 機能 : 入力と出力の周波数成分ごとの位相差を求めます 書式 :phaseresp4(indat,outdat,inc,outc) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) inc:0にする入力成分の 最大値に対する比率 outc:0にする出力成分の 最大値に対する比率 戻り値 : 周波数成分ごとの位相の差からなる数列 ただし 出力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *outc) 以上の大きさを持つ周波数成分で ピーク値となっているもの以外は 0 としています さらに入力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *inc) 以上の大きさを持つ周波数成分以外は 0 としています メニュー項目 : /* 位相応答 (ip,op) */ 機能 : 入力と出力の周波数成分ごとの位相差を求めます 書式 :phaseresp5(indat,outdat,inc,outc) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) inc:0にする入力成分の 最大値に対する比率 outc:0にする出力成分の 最大値に対する比率 戻り値 : 周波数成分ごとの位相の差からなる数列 ただし 出力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *outc) 以上の大きさを持つ周波数成分で ピーク値となっているもの以外は 0 としています さらに入力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *inc) 以上の大きさを持つ周波数成分でピーク値となっているもの以外は 0 としています 例 : w1:gline(2^10,1/2^10,2*pi,0) w2:cos(w1*4)-sin(35*w1)+sin(200*w1) w3:2*sin(4*w1)-3*cos(35*w1)+4*cos(200*w1) w4:phaseresp2(w2,w3,0.0,0.001) 位相応答 (om) w5:phaseresp(w2,w3) 位相応答 w6:cphaseresp(w2,w3,5.0) 位相応答 ( 連続 ) 27

31 説明 : 位相の計算は -π から +π 間での間で与えられます w4 では 出力の各周波数成分のうち spectrum の最大値に比べて弱いものを 0 に置き換えて計算します w5 は単純な位相応答 ここでは計算誤差の影響が強くでています w6 は w5 の結果を連続化したもの この例の場合は不適当ですが 周波数に関して連続的な現象を扱うときには役立ちます メニュー項目 : /* 位相応答 ( 連続 ) */ /* 位相差を連続化する */ 機能 : 入力と出力の周波数成分ごとの位相差連続化します 書式 :chphaseresp(indat,outdat,dlt) indat: 入力データ ( データサイズは 2^k) outdat: 出力データ ( データサイズは 2^k) dlt: dlt より大きな位相差の変化は連続化します 戻り値 : 周波数成分ごとの位相の差を連続化したもの 説明 : 位相差は π から少しでも増加すると 急に減少し -π に変化します そのときに この変化量は 2π に近い値となります このような場合は -π にするのではなく π+α となるようにしています その結果 急激な変化を連続化してとらえています 位相の変化が dlt 以上の量で急激に変化するときは その変化量を修正します たとえば から に変化したときは *pi = に変えるのです これによって 3.00 から に連続的に変化するようになります 28

32 3.4 Multi 出力位相応答 メニュー項目 : /* Mlt 位相応答 ( 単純 ) */ 機能 : 入力と出力の周波数成分ごとの位相差を求めます 書式 :mphaseresp(indat,outdat) indat: 入力データ (1-ch) outdat: 出力データ (k-ch) ( データサイズは 2^k) ( 各チャンネルのデータサイズは 2^k) 戻り値 : 周波数成分ごとの位相の差からなる数列 説明 : 位相差は -π から +π 間での間で与えられます メニュー項目 : /* Mlt 位相応答 (om) */ 機能 : 入力と出力の周波数成分ごとの位相差を求めます 書式 :mphaseresp2(indat,outdat,outc) indat: 入力データ (1-ch) ( データサイズは 2^k) outdat: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) outc:0にする出力成分の 最大値に対する比率 戻り値 : 周波数成分ごとの位相の差からなる数列ただし 出力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *outc) 以上の大きさを持つ周波数成分以外は 0 としています メニュー項目 : /* Mlt 位相応答 (im,om) */ 機能 : 入力と出力の周波数成分ごとの位相差を求めます 書式 :mphaseresp3(indat,outdat,inc,outc) indat: 入力データ (1-ch) ( データサイズは 2^k) outdat: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) inc:0にする入力成分の 最大値に対する比率 outc:0にする出力成分の 最大値に対する比率 戻り値 : 周波数成分ごとの位相の差からなる数列 ただし 出力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *outc) 以上の大きさを持つ周波数成分以外は 0 としています さらに入力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *inc) 以上の大きさを持つ周波数成分以外は 0 としています 29

33 メニュー項目 : /* Mlt 位相応答 (im,op) */ 機能 : 入力と出力の周波数成分ごとの位相差を求めます 書式 :mphaseresp4(indat,outdat,inc,outc) indat: 入力データ (1-ch) ( データサイズは 2^k) outdat: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) inc:0にする入力成分の 最大値に対する比率 outc:0にする出力成分の 最大値に対する比率 戻り値 : 周波数成分ごとの位相の差からなる数列 ただし 出力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *outc) 以上の大きさを持つ周波数成分で ピーク値となっているもの以外は 0 としています さらに入力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *inc) 以上の大きさを持つ周波数成分以外は 0 としています メニュー項目 : /* Mlt 位相応答 (ip,op) */ 機能 : 入力と出力の周波数成分ごとの位相差を求めます 書式 :mphaseresp5(indat,outdat,inc,outc) indat: 入力データ (1-ch) ( データサイズは 2^k) outdat: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) inc:0にする入力成分の 最大値に対する比率 outc:0にする出力成分の 最大値に対する比率 戻り値 : 周波数成分ごとの位相の差からなる数列 ただし 出力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *outc) 以上の大きさを持つ周波数成分で ピーク値となっているもの以外は 0 としています さらに入力波形に含まれる周波数成分のうち ( スペクトラムの最大値 *inc) 以上の大きさを持つ周波数成分でピーク値となっているもの以外は 0 としています メニュー項目 : /* Mlt 位相応答 ( 連続 ) */ /* 位相差を連続化する */ 機能 : 入力と出力の周波数成分ごとの位相差連続化します 書式 :mchphaseresp(indat,outdat,dlt) indat: 入力データ (1-ch) ( データサイズは 2^k) outdat: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) dlt: dlt より大きな位相差の変化は連続化します 戻り値 : 周波数成分ごとの位相の差を連続化したもの 説明 : 位相差は π から少しでも増加すると 急に減少し -π に変化します そのときに この変化量は 2π に近い値となります このような場合は -π にするのではなく π+α となるようにしています 30

34 その結果 急激な変化を連続化してとらえています 位相の変化が dlt 以上の量で急激に変化するときは その変化量を修正します たとえば から に変化したときは *pi = に変えるのです これによって 3.00 から に連続的に変化するようになります 逆に から に急激に変化したときは *pi = に変えるのです これによって から に連続的に変化するようになります 31

35 3.5 コヒーレンス メニュー項目 : /* コヒーレンス ( 単純 ) */ 機能 : 入力と出力のコヒーレンスを求めます 書式 :ucoh(sx,sy,len) sx: 入力データ sy: 出力データ len: 切り出す長さ ( データサイズは 2^k) ( データサイズは 2^k) 戻り値 : コヒーレンス 説明 : データの長さは全て 2 のべき乗にして下さい 条件 : /* (sx,sy は同じ長さ len は切出す長さ ) */ /* 条件 : length(sx) >= 2*len, len = 2^k */ に合うようにして下さい 周波数成分が連続的に変化するような場合にこれを使用すると良いと思います メニュー項目 : /* コヒーレンス (om) */ 機能 : 入力と出力のコヒーレンスを求めます 書式 :ucoh2(sx,sy,len,yz) sx: 入力データ ( データサイズは 2^k) sy: 出力データ ( データサイズは 2^k) len: 切り出す長さ yz:0にする出力成分の最大値に対する比率 戻り値 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 説明 : データの長さは全て 2 のべき乗にして下さい 条件 : /* (sx,sy は同じ長さ len は切出す長さ ) */ /* 条件 : length(sx) >= 2*len, len = 2^k */ に合うようにして下さい 32

36 メニュー項目 : /* コヒーレンス (im,om) */ 機能 : 入力と出力のコヒーレンスを求めます 書式 :ucoh3(sx,sy,len,xz,yz) sx: 入力データ ( データサイズは 2^k) sy: 出力データ ( データサイズは 2^k) len: 切り出す長さ xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 説明 : データの長さは全て 2 のべき乗にして下さい 条件 : /* (sx,sy は同じ長さ len は切出す長さ ) */ /* 条件 : length(sx) >= 2*len, len = 2^k */ に合うようにして下さい メニュー項目 : /* コヒーレンス (im,op) */ 機能 : 入力と出力のコヒーレンスを求めます 書式 :ucoh4(sx,sy,len,xz,yz) sx: 入力データ ( データサイズは 2^k) sy: 出力データ ( データサイズは 2^k) len: 切り出す長さ xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 説明 : データの長さは全て 2 のべき乗にして下さい 条件 : /* (sx,sy は同じ長さ len は切出す長さ ) */ /* 条件 : length(sx) >= 2*len, len = 2^k */ に合うようにして下さい 33

37 メニュー項目 : /* コヒーレンス (ip,op) */ 機能 : 入力と出力のコヒーレンスを求めます 書式 :ucoh5(sx,sy,len,xz,yz) sx: 入力データ ( データサイズは 2^k) sy: 出力データ ( データサイズは 2^k) len: 切り出す長さ xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 説明 : データの長さは全て 2 のべき乗にして下さい 条件 : /* (sx,sy は同じ長さ len は切出す長さ ) */ /* 条件 : length(sx) >= 2*len, len = 2^k */ に合うようにして下さい メニュー項目 : /* コヒーレンス (H, 単純 ) */ 機能 : 入力と出力のコヒーレンスを求めます 書式 :ucohh(sx,sy) sx: 入力データ sy: 出力データ ( データサイズは 2^k) ( データサイズは 2^k) 戻り値 : コヒーレンス 説明 : 長さは全て 2 のべき乗にして下さい 切り出す長さは 入力信号の長さの半分です 条件 : sx,sy は同じ長さで 2 のべき乗 切出す長さはその半分です 周波数成分が連続的に変化するような場合にこれを使用すると良いと思います 34

38 メニュー項目 : /* コヒーレンス (H,oM) */ 機能 : 入力と出力のコヒーレンスを求めます 書式 :ucoh2h(sx,sy,yz) sx: 入力データ ( データサイズは 2^k) sy: 出力データ ( データサイズは 2^k) yz:0にする出力成分の最大値に対する比率 戻り値 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 切り出す長さは 入力信号の長さの半分です 説明 : データの長さは全て 2 のべき乗にして下さい 条件 : sx,sy は同じ長さで 2 のべき乗 切出す長さはその半分になります メニュー項目 : /* コヒーレンス (H,iM,oM) */ 機能 : 入力と出力のコヒーレンスを求めます 書式 :ucoh3h(sx,sy,xz,yz) sx: 入力データ ( データサイズは 2^k) sy: 出力データ ( データサイズは 2^k) xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 切り出す長さは 入力信号の長さの半分です 説明 : データの長さは全て 2 のべき乗にして下さい 条件 : sx,sy は同じ長さで 2 のべき乗 切出す長さはその半分になります メニュー項目 : /* コヒーレンス (H,iM,oP) */ 機能 : 入力と出力のコヒーレンスを求めます 書式 :ucoh4h(sx,sy,xz,yz) sx: 入力データ ( データサイズは 2^k) sy: 出力データ ( データサイズは 2^k) xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 35

39 切り出す長さは 入力信号の長さの半分です 説明 : データの長さは全て 2 のべき乗にして下さい 条件 : sx,sy は同じ長さで 2 のべき乗 切出す長さはその半分になります メニュー項目 : /* コヒーレンス (H,iP,oP) */ 機能 : 入力と出力のコヒーレンスを求めます 書式 :ucoh5h(sx,sy,xz,yz) sx: 入力データ ( データサイズは 2^k) sy: 出力データ ( データサイズは 2^k) xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 切り出す長さは 入力信号の長さの半分です 説明 : データの長さは全て 2 のべき乗にして下さい 条件 : sx,sy は同じ長さで 2 のべき乗 切出す長さはその半分になります 例 : w1:gline(2^10,1/2^10,2*pi,0) w2:sin(w1*4)-sin(35*w1)+sin(200*w1) w3:2*cos(w1*4)-3*cos(35*w1)+4*cos(200*w1) w4:ucoh(w2,w3,256) コヒーレンス ( 単純 ) w5:ucoh2(w2,w3,256,0.0,0.001) コヒーレンス (im,om) w5:magresp5(w2,w3,0.0,0.001,1) 振幅応答 (im,om) 説明 :w4 では 計算誤差のために いろいろな周波数に対して値が計算され 不正確なグラフになっています w5 では 本来存在しない周波数成分については 0 として扱っていますので 4Hz 35Hz 200Hz の部分だけが表示されています 切り出した入力信号の第 1 列のスペクトラムをとり その最大値の xz 倍よりも小さい周波数 36

40 成分は 0 に変えます 切り出した出力信号の第 1 列のスペクトラムをとり その最大値の yz 倍よりも小さい周波数成分は 0 に変えます 条件 : 信号の長さは全て 2 のべき乗 sx,sy は同じ長さにしてください len は切出す長さ 37

41 3.6 Multi 出力コヒーレンス メニュー項目 : /* Mlt コヒーレンス ( 単純 ) */ 機能 :1 チャンネル入力と多チャンネル出力の間のコヒーレンスを求めます 書式 :mucoh(sx,sy,len) sx: 入力データ (1-ch) sy: 出力データ (k-ch) len: 切り出す長さ ( データサイズは 2^k) ( 各チャンネルのデータサイズは 2^k) 戻り値 : 多チャンネルのコヒーレンス 説明 : データの長さは全て 2 のべき乗にして下さい 条件 : /* (sx,sy は同じ長さ len は切出す長さ ) */ /* 条件 : length(sx) >= 2*len, len = 2^k */ に合うようにして下さい 周波数成分が連続的に変化するような場合にこれを使用すると良いと思います メニュー項目 : /* Mlt コヒーレンス (om) */ 機能 :1 チャンネル入力と多チャンネル出力のコヒーレンスを求めます 書式 :mucoh2(sx,sy,len,yz) sx: 入力データ (1-ch) ( データサイズは 2^k) sy: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) len: 切り出す長さ yz:0にする出力成分の最大値に対する比率 戻り値 : 多チャンネルのコヒーレンス 説明 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています データの長さは全て 2 のべき乗にして下さい 条件 : /* (sx,sy は同じ長さ len は切出す長さ ) */ /* 条件 : length(sx) >= 2*len, len = 2^k */ に合うようにして下さい 38

42 メニュー項目 : /* Mlt コヒーレンス (im,om) */ 機能 :1 チャンネル入力と多チャンネル出力の間のコヒーレンスを求めます 書式 :mucoh3(sx,sy,len,xz,yz) sx: 入力データ (1-ch) ( データサイズは 2^k) sy: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) len: 切り出す長さ xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 多チャンネルのコヒーレンス 説明 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています データの長さは全て 2 のべき乗にして下さい 条件 : /* (sx,sy は同じ長さ len は切出す長さ ) */ /* 条件 : length(sx) >= 2*len, len = 2^k */ に合うようにして下さい メニュー項目 : /* Mlt コヒーレンス (im,op) */ 機能 :1 チャンネル入力と多チャンネル出力の間のコヒーレンスを求めます 書式 :mucoh4(sx,sy,len,xz,yz) sx: 入力データ (1-ch) ( データサイズは 2^k) sy: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) len: 切り出す長さ xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 多チャンネルのコヒーレンス 説明 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています データの長さは全て 2 のべき乗にして下さい 条件 : /* (sx,sy は同じ長さ len は切出す長さ ) */ /* 条件 : length(sx) >= 2*len, len = 2^k */ に合うようにして下さい 39

43 メニュー項目 : /* Mlt コヒーレンス (ip,op) */ 機能 :1 チャンネル入力と多チャンネル出力の間のコヒーレンスを求めます 書式 :mucoh5(sx,sy,len,xz,yz) sx: 入力データ (1-ch) ( データサイズは 2^k) sy: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) len: 切り出す長さ xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 多チャンネルのコヒーレンス 説明 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています データの長さは全て 2 のべき乗にして下さい 条件 : /* (sx,sy は同じ長さ len は切出す長さ ) */ /* 条件 : length(sx) >= 2*len, len = 2^k */ に合うようにして下さい メニュー項目 : /* Mlt コヒーレンス (H, 単純 ) */ 機能 :1 チャンネル入力と多チャンネル出力の間のコヒーレンスを求めます 書式 :mucohh(sx,sy) sx: 入力データ (1-ch) sy: 出力データ (k-ch) ( データサイズは 2^k) ( 各チャンネルのデータサイズは 2^k) 戻り値 : 多チャンネルのコヒーレンス 説明 : 長さは全て 2 のべき乗にして下さい 切り出す長さは 入力信号の長さの半分です 条件 : sx,sy は同じ長さで 2 のべき乗 切出す長さはその半分です 周波数成分が連続的に変化するような場合にこれを使用すると良いと思います 40

44 メニュー項目 : /* Mlt コヒーレンス (H,oM) */ 機能 :1 チャンネル入力と多チャンネル出力の間のコヒーレンスを求めます 書式 :mucoh2h(sx,sy,yz) sx: 入力データ (1-ch) ( データサイズは 2^k) sy: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) yz:0にする出力成分の最大値に対する比率 戻り値 : 多チャンネルのコヒーレンス 説明 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 切り出す長さは 入力信号の長さの半分です データの長さは全て 2 のべき乗にして下さい 条件 : sx,sy は同じ長さで 2 のべき乗 切出す長さはその半分になります メニュー項目 : /* Mlt コヒーレンス (H,iM,oM) */ 機能 :1 チャンネル入力と多チャンネル出力の間のコヒーレンスを求めます 書式 :mucoh3h(sx,sy,xz,yz) sx: 入力データ (1-ch) ( データサイズは 2^k) sy: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 多チャンネルのコヒーレンス 説明 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 切り出す長さは 入力信号の長さの半分です データの長さは全て 2 のべき乗にして下さい 条件 : sx,sy は同じ長さで 2 のべき乗 切出す長さはその半分になります メニュー項目 : /* Mlt コヒーレンス (H,iM,oP) */ 機能 :1 チャンネル入力と多チャンネル出力の間のコヒーレンスを求めます 書式 :mucoh4h(sx,sy,xz,yz) sx: 入力データ (1-ch) ( データサイズは 2^k) sy: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 多チャンネルのコヒーレンス 41

45 説明 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 切り出す長さは 入力信号の長さの半分です データの長さは全て 2 のべき乗にして下さい 条件 : sx,sy は同じ長さで 2 のべき乗 切出す長さはその半分になります メニュー項目 : /* Mlt コヒーレンス (H,iP,oP) */ 機能 :1 チャンネル入力と多チャンネル出力の間のコヒーレンスを求めます 書式 :mucoh5h(sx,sy,xz,yz) sx: 入力データ (1-ch) ( データサイズは 2^k) sy: 出力データ (k-ch) ( 各チャンネルのデータサイズは 2^k) xz:0にする入力成分の最大値に対する比率 yz:0にする出力成分の最大値に対する比率 戻り値 : 多チャンネルのコヒーレンス 説明 : 入力 出力ともに特別の周波数成分しか含まない場合のコヒーレンスの計算に適しています 切り出す長さは 入力信号の長さの半分です データの長さは全て 2 のべき乗にして下さい 条件 : sx,sy は同じ長さで 2 のべき乗 切出す長さはその半分になります 42

46 3.7 フィルター メニュー項目 : /* Filter(0*5) 作成 */ 機能 : 特定の周波数の波と その高調波を取り除くためのフィルターを作成します 書式 :mkfilt05(sx,ff) sx:fft データ ff: 除去する周波数 ( データサイズは 2^k) 戻り値 : ある周波数成分を除くためのフィルター 例 : w1:gline(1024,1/1024,2*pi,0) w2:sin(w1)-sin(50*w1)+sin(60*w1) w3:fft(w2) w4:mkfilt05(w3,50.0) (Filetr(0*5) 作成 ) w5:ifft(w3*w4) 説明 : この w4 が 50Hz 用のフィルタです これを FFT の結果にかけてから ifft() をとります この結果 50Hz の成分が取り除かれた波が現れます電源 (50Hz 60Hz) からの雑音を取り除くためのフィルタです 0 が 5 個並ぶので 削除する周波数成分がやや多くなります なお 周波数は 自由に設定できます メニュー項目 : /* Filter(0*5) 適用 */ 機能 : 特定の周波数の波とその高調波を取り除きます 書式 :adfilt05(sx,ff) sx: 入力データ ff: 除去する周波数 ( データサイズは 2^k) 戻り値 : 特定の周波数成分を除いたもの 例 : w1:gline(1024,1/1024,2*pi,0) w2:sin(w1)-sin(50*w1)+sin(60*w1) w6:adfilt05(w2,50.0) (Filetr(0*5) 適用 ) 説明 : この w6 が 50Hz の成分を除去した結果です 電源 (50Hz 60Hz) からの雑音を取り除くのに役立ちます 0 が 5 個並ぶので 削除する周波数成分がやや多くなります なお 周波数は 自由に設定できます 43

47 メニュー項目 : /* Filter(0*3) 作成 */ 機能 : 特定の周波数の波とその高調波を取り除くためのフィルタをつくる 書式 :mkfilt03(sx,ff) sx:fft データ ff: 取り除く周波数 ( データサイズは 2^k) 戻り値 : 特別の周波数成分を除くためのフィルター 例 : w1:gline(1024,1/1024,2*pi,0) w2:sin(w1)-sin(50*w1)+sin(60*w1) w3:fft(w2) w4:mkfilt03(w3,60.0) (Filetr(0*3) 作成 ) w5:ifft(w3*w4) 説明 : 特定の周波数成分とその高調波を取り除くためのフィルタです 電源 (50Hz 60Hz) からの雑音を取り除くのに適します 0 が 3 個並ぶだけですから 削除する周波数成分はより少なくなっています FFT データに対して使用します この結果できたフィルタを FFT データに掛けてから IFFT() をとれば もとの波から特定の周波数成分および その高調波を取り除いたものが得られます なお 周波数は自由に設定できます メニュー項目 : /* Filter(0*3) 適用 */ 機能 : 入力データから 特定の周波数の波とその高調波を取り除きます 書式 :adfilt03(sx,ff) sx: 入力データ ff: 取り除く周波数 ( データサイズは 2^k) 戻り値 : 特定の周波数成分と その高調波を除いたもの 交流 (50Hz 60Hz) 電源からの雑音を取り除くのに役立ちます 0 が 3 個並ぶだけですから 削除する周波数成分はより少なくなっています なお 周波数は自由に設定できます 44

48 例 : w1:gline(1024*8,1/(1024*8),2*pi,0) w2:-sin(22*w1)+sin(25*w1) w3:adfilt05(w2,22.0) w4:adfilt03(w2,24.0) w5:adfilt05(w2,24.0) 説明 :w2 は 22Hz と 25Hz の波からできています w3 はそれから 22Hz の成分を除いたものです w4 は 60Hz フィルタを使って 24Hz の成分を取り除いたものです w5 は 50Hz フィルタを使って 24hz の成分を取り除いたものです FFT の成分の 5 つずつを 0 にしている 50Hz フィルターでは 22Hz25Hz ともに除去されてしまいますが 3 個ずつ 0 にしている 60Hz フィルタでは 22Hz25Hz の両方とも残っています 残念ながら 23Hz 成分だけを除こうととすると 22Hz の成分も取り除かれてしまいす 45

49 3.8 PSD の平均 (X(0)=0) メニュー項目 : /* PSD の平均 ( 窓なし ) */ /* (sx,sy は同じ長さ len は切出す長さ ) */ /* length(sx) >= 2*len, len = 2^k */ 機能 : 入力データを長さ len で切り出し len/2 の重なりを持たせて何度か切り出します それぞれの PSD を取ってから平均しています 書式 :meanpsd(sx,len) sx: 入力データ len: 切り出す長さ ( データサイズは 2^k) 戻り値 :PSD の平均 説明 : 窓関数はかけていません 平均値は最初に引いています メニュー項目 : /* PSD の平均 (hanning) */ /* (sx,sy は同じ長さ len は切出す長さ, hanning 窓付き */ /* 条件 : length(sx) >= 2*len, len = 2^k */ 機能 : 入力データを長さ len で切り出し len/2 の重なりを持たせて何度か切り出します それぞれの PSD を取ってから平均しています 書式 :meanpsd2(sx,len) sx: 入力データ len: 切り出す長さ ( データサイズは 2^k) 戻り値 :PSD の平均 説明 : 窓関数 hanning 窓です 平均値は最初に引いています メニュー項目 : /* PSD の平均 (hamming) */ /* (sx,sy は同じ長さ len は切出す長さ ) hamming 窓付き */ /* 条件 : length(sx) >= 2*len, len = 2^k */ 機能 : 入力データを長さ len で切り出し len/2 の重なりを持たせて何度か切り出します それぞれの PSD を取ってから平均しています 書式 :meanpsd3(sx,len) sx: 入力データ len: 切り出す長さ ( データサイズは 2^k) 戻り値 :PSD の平均 説明 : 窓関数 hamming 窓です 平均値は最初に引いています 46

50 メニュー項目 : /* PSD の平均 (kaiser) */ /* (sx,sy は同じ長さ len は切出す長さ ) kaiser 窓付き */ /* 条件 : length(sx) >= 2*len, len = 2^k */ 機能 : 入力データを長さ len で切り出し len/2 の重なりを持たせて何度か切り出します それぞれの PSD を取ってから平均しています 書式 :meanpsd4(sx,len) sx: 入力データ len: 切り出す長さ ( データサイズは 2^k) 戻り値 : 窓関数 kaiser 窓をつけて PSD の平均を計算したもの 平均値は最初に引いています 例 : w1:gline(2^7,1,2*pi/2^7,0) w2:sin(w1)-sin(3*w1)+sin(9*w1) w3:meanpsd(w2,256) PSD の平均 w4:meanpsd2(w2,256) PSD の平均 (hanning) w5:meanpsd3(w2,256) PSD の平均 (hamming) w6:meanpsd4(w2,256) PSD の平均 (kaiser) w7:psd(w2) 普通の PSD 説明 : 窓関数の影響で値が小さくなっています メニュー項目 : /* PSD の平均 ( 各列 )( 窓なし ) */ /* (len は切出す長さ len = 2^k */ 47

51 機能 : 多チャンネルデータの各列に対して 長さ len で切り出し olp の重なりを持たせてさらに何度か切り出します 切り出したそれぞれのデータの PSD を取ってから平均しています 書式 :colpsd(md,len,olp) md: 入力データ ( 多チャンネルデータ ) len: 切り出す長さ olp: 重なる長さ 戻り値 : 各列ごとの PSD の平均 からなる多チャンネルデータ 説明 : 窓関数はかけていません 平均値は最初に引いています メニュー項目 : /* PSD の平均 ( 各列 )(hanning) */ /* (len は切出す長さ, hanning 窓付き */ /* 条件 : len = 2^k */ 機能 : 多チャンネルデータの各列に対して 長さ len で切り出し olp の重なりを持たせてさらに何度か切り出します 切り出したそれぞれのデータの PSD を取ってから平均しています 書式 :colpsd2(md,len,olp) md: 入力データ ( 多チャンネルデータ ) len: 切り出す長さ olp: 重なる長さ 戻り値 : 各列ごとの PSD の平均 からなる多チャンネルデータ 説明 : 窓関数 hanning 窓です 平均値は最初に引いています メニュー項目 : /* PSD の平均 ( 各列 )(hamming) */ /* (len は切出す長さ ) hamming 窓付き */ /* 条件 : len = 2^k */ 機能 : 多チャンネルデータの各列に対して 長さ len で切り出し olp の重なりを持たせてさらに何度か切り出します 切り出したそれぞれのデータの PSD を取ってから平均しています 書式 :colpsd3(md,len,olp) md: 入力データ ( 多チャンネルデータ ) len: 切り出す長さ olp: 重なる長さ 戻り値 : 各列ごとの PSD の平均 からなる多チャンネルデータ 説明 : 窓関数 hamming 窓です 平均値は最初に引いていますメニュー項目 : /* PSD の平均 ( 各列 )(kaiser) */ /* (len は切出す長さ 条件 : len = 2^k *) kaiser 窓付き */ 48

52 機能 : 多チャンネルデータの各列に対して 長さ len で切り出し olp の重なりを持たせてさらに何度か切り出します 切り出したそれぞれのデータの PSD を取ってから平均しています 書式 :colpsd4(md,len,olp) md: 入力データ ( 多チャンネルデータ ) len: 切り出す長さ olp: 重なる長さ 戻り値 : 各列ごとの PSD の平均 からなる多チャンネルデータ 説明 : 窓関数 kaiser 窓です 平均値は最初に引いています 3.9 PSD の平均 (X(0)!=0) メニュー項目 : /* PSD の平均 (X0!=0) */ /* (len は切出す長さ ) length(sx) >= 2*len, len = 2^k */ 機能 : 入力データを長さ len で切り出し len/2 の重なりを持たせて何度か切り出します それぞれの PSD を取ってから平均しています 書式 :meanpsdx0(sx,len) sx: 入力データ len: 切り出す長さ ( データサイズは 2^k) 戻り値 :PSD の平均 説明 : 窓関数はかけていません 切り出した波形の直流成分はそのまま残してあります 従って直流成分の強さが大きな場合は波の強さが見ずらくなります メニュー項目 : /* PSD の平均 ( 各列 ) (X0!=0) */ /* (len は切出す長さ len = 2^k */ 機能 : 多チャンネルデータの各列に対して 長さ len で切り出し olp の重なりを持たせてさらに何度か切り出します 切り出したそれぞれのデータの PSD を取ってから平均しています 書式 :colpsdx0(md,len,olp) md: 入力データ ( 多チャンネルデータ ) len: 切り出す長さ ( この長さは 2 のべき乗にしてください ) olp: 重なる長さ 戻り値 : 各列ごとの PSD の平均 からなる多チャンネルデータ 説明 : 窓関数はかけていません 切り出した波形の直流成分はそのまま残してあります 従って直流成分の強さが大きな場合は波の強さが見ずらくなります 49

53 3.10 サンプリング メニュー項目 : /* サンプリング */ 機能 : 入力データから pt 個のサンプルをランダムに取り出します 書式 :usampling(s0,pt,seed) s0: 与えられたデータ pt: サンプリング数 seed: 乱数発生のための種数 ( 自然数 ) 戻り値 : 最初のデータから pt 個だけ抽出し最初に並んでいた順に並べてかえします 説明 : 必要な個数 (pt) だけサンプリングして抜き出す のですが サンプリングする数が全体数に比べて多いと失敗します 失敗したときは 0 の列が返ります 目安は 全体の 1/2 くらいまでですが これが限界というわけではありません 乱数発生のための種数は自然数を適当に入れて下さい 乱数は疑似乱数ですので 種数によって乱数の列は決定されてしまいます メニュー項目 : /* サンプリング ( 番号付 ) */ 機能 : 入力データから pt 個のサンプルをランダムに取り出します 書式 :usampling2(s0,pt,seed) s0: 与えられたデータ pt: サンプリング数 seed: 乱数発生のための種数 ( 自然数 ) 戻り値 : 最初のデータから pt 個だけ抽出し最初に並んでいた順に並べてかえします 左側には 最初に並んでいたときの番号がつきます 説明 : 必要な個数 (pt) だけサンプリングして抜き出す のですが サンプリングする数が全体数に比べて多いと失敗します 失敗したときは 0 の列が返ります 目安は 全体の 1/2 くらいまでですが これが限界というわけではありません 乱数発生のための種数は自然数を適当に入れて下さい 乱数は疑似乱数ですので 種数によって乱数の列は決定されてしまいます サンプリングと同じですが 最初に並んでいたときの番号が付いています 50

54 3.11 その他 メニュー項目 : /* 列データの拡大 シフト */ 機能 : 複数列データのある列の値だけを拡大 シフト 書式 :shiext(m1,nb,ex,sh){ m1: 行列形テ ータ nb: 列の番号 (1,2,3,..) ex: 拡大率 sh: シフト量 戻り値 :nb で指定された列だけが ex 倍され さらに sh を加えられた行列 2 列以上の列を持つ行列形のテ ータは一つの窓で複数の曲線として表示されます そのうちの特定の曲線 ( これはテーブル形の表示での列に対応します ) について グラフの y 軸方向に関する拡大とシフトを行います 例 : w1:gline(2^10,1,2*pi/2^10,0) w2:sin(w1)-sin(3*w1)+sin(9*w1) w3:2*cos(w1)-3*cos(3*w1)+4*cos(9*w1) w4:ravel(w2,w3,-2*w2) w5:shiext(w4,2,2.0,10) 列データの拡大 シフト 説明 :w5 では 第 2 列のデータが 2 倍に拡大され さらに 10 だけ上にシフトされています 51

55 メニュー項目 : /* 前回平均値比較機能 */ 機能 : 前回までの測定の平均値と今回の測定値との差を比較判定し 急激な変動分のみをとらえます 書式 :pmc(sx) sx: 入力データ 戻り値 : 1 列目は 番号 2 列目は w1 の値 ( 入力値 D1,D2,...) 3 列目は N1,N2,... 4 列目は 最初は必ず 0 2 つ目が D2-N1,... となっている行列 説明 : 前回までの測定の平均値と今回の測定値との差を比較判定し 急激な変動分のみをとらえます 計算式は以下のようになっています 測定回数 入力値 次回入力時の比較値 表示値 1 D1 N1=D1 0 2 D2 N2=0.5*(N1+D2) D2-N1 3 D3 N3=0.5*(N3+D3) D3-N2 N Dn Nn=0.5*(Nn-1+Dn) Dn-Nn-1 w1 に データ D1,D2,D3.. があるとします w2 で pmh(w1) とすれば W2 の各列は 1 列目は 番号 ( 測定回数 ) 2 列目は 入力値 D1,D2,...) 3 列目は N1,N2,... 4 列目は 最初は必ず 0 2 つ目が D2-N1,... になります この 4 列目が 目的の数値です 52

56 メニュー項目 : /* 相互相関関数 ( 正規化 ) */ 機能 : データサイズが同じ場合 二つのデータを正規化して相互相関を求めます 書式 :nzxcorr(sx,sy) sx: 第 1 のデータ sy: 第 2 のデータ 戻り値 : 同じデータについて計算すると最大値が 1 でグラフは紡錘形になります 例 : w1:gline(2^10,1/2^10,2*pi,0) w2:sin(w1)-sin(3*w1)+sin(9*w1) w3:2*sin(w1)-3*sin(3*w1)+4*sin(9*w1) w4:crosscor(w2,w3) 相互相関関数 w5:mvinteg(w3,155) 移動積分 w6:nzxcorr(w2,w2) 相互相関関数 ( 正規化 ) 説明 : 最初にデータを正規化 ( 平均を引いて 標準偏差で割る ) しています 計算の時には 両端に 0 の列を加えた形で計算していますので 結果のグラフは紡錘形になります メニュー項目 : /* 移動積分 */ 機能 : 与えられた幅で移動しながら積分して行きます 書式 :mvinteg(s1,pt) s1: 入力データ pt: 積分する幅 戻り値 : 与えられた幅で積分をして その積分結果からなる数列を返します 説明 : 移動平均の積分版です 53

57 メニュー項目 : /* 振幅 /(2πf) */ 機能 : 各周波数成分ごとの振幅を ( 周波数 *2π) で割ったものを求めます 書式 :spde2pif(s1) s1: 入力データ 戻り値 : 各周波数成分ごとの振幅を ( 周波数 *2π) で割ったものからなる数列を返します 説明 : 上のグラフは w2: -sin(22*w1)+32.6*sin(250*w1)+45*cos(2500*w1) w3: spectrum(w2) w4: spde2pif(w2) ( 振幅 /(2πf)) であり w2のスペクトラムでは 各周波数成分の振幅が とそのまま現れています w4では 1/(2*π* 周波数 ) の値が並びますので 1/(2*π*22)= /(2*π*250)= /(2*π*2500)= となっています なお計算は 平均を引いてからスペクトラムをとり その結果を変形しています 54

58 メニュー項目 : /* 心拍数と血圧 */ 機能 : 心電図データと血圧のデータに対して 心拍数と血圧の最大値 最小値 平均値を同時に求めます 書式 : 心拍数 : heart1(s1, s4, s5) 最大値 : blood2(s2) 最小値 : blood3(s2) 平均値 : blood4(s2) 戻り値 : 連続した 4 つのウインドウのうち最初にものから順に心拍数 血圧の最大値 最小値 平均値が得られます 説明 : 関数の引数などを指定するパネルの項目は次のような意味を持ちます 心電図データ : ここは w1 の形で心電図データが入っているウインドウを指定します 基準値 : 心電図のピークをとらえるための基準値です この値を超えたものはピーク値として認識されます up-down(0,1,2): 心電図のグラフが Y= 基準値のグラフと交差する回数を求めます 0 のときは < グラフ値 = 基準値 > となる回数 1 のときは グラフが増加しながら < グラフ値 = 基準値 > となる回数 2 のときは グラフが減少しながら < グラフ値 = 基準値 > となる回数 心拍数 W: ここは 3 の形で心拍数を出力するウインド番号の数値のみを記入 血圧データ : ここは w2 の形で血圧データが入っているウインドウを指定します 最大値 W: ここは 4 の形で最大値を出力するウインドウ番号の数値のみを記入 最小値 W: ここは 5 の形で最小値を出力するウインドウ番号の数値のみを記入 平均値 W: ここは 6 の形で平均値を出力するウインドウ番号の数値のみを記入 結果は下図のようになります 55

59 東京都千代田区飯田橋 TEL FAX fluid.co.jp