http://kojoryoshu83.edoblog.net/%E7%A7%91%E5%AD%A6%E6%8A%80%E8%A1%93/%E5%9C%B0%E9%9C%87%E6%B3%A2%E5%BD%A2%E8%A1%A8%E7%A4%BA%E3%83%97%E3%83%AD%E3%82%B0%E3%83%A9%E3%83%A0%EF%BC%92地震波形表示プログラム２
前回のW32BPFtoPSプログラムに続いて、今回はW32IIRtoPSプログラムの説明をする。
W32IIRtoPSプログラムはwin32データの中の選択されたチャネルにIIRフィルタを掛けてポストスクリプト・ファイル描画するプログラムである。フイルタの内容はテキストファイルで記述する(注3)。
機能：
・IIRフィルタは、チャンル毎に１～２次を最大16段まで可能。
・フィルタとwinチャンネルの指定は、テキストファイルで指定
・適当なダンピング定数のフィルタを段重ねすることにより、任意のフィルタを指定することが出来る。
書式：
W32IIRtoPS ch w32 ps IIR [t1]［tw] [s]
ch ：チャネルテーブル
w32 ：win32データファイル
ps ：作成されるポストスクリプト・ファイル
IIR ：IIRフィルタ設定ファイル、０の場合は設定ファイル無し
t1 ：記録開始時間［秒］
tw ：記録時間幅［秒］
s ：Sの場合作業ファイルを削除しない
"W32IIRtoPS" のみの場合は、上記help が表示される。
"W32IIRtoPS ？" とした場合は、各種フィルタの例が表示される。
"W32IIRtoPS ch ？" とした場合は、tmp.fltの名のバンドパスフイルタのテキストファイルが作成されるので、このテキストファイルに必要な事項を記入して使用することが出来る。
フィルタのテキストファイル書式：(チャンネル毎)
1行目：chID0,chID1,[g],[a],[t]
chID0 = チャンネルID（４バイトHEX)
chID0 = 変更するチャンネルID（４バイトHEX)。
このプログラムでは元のwinチャンネルIDのみ。
g = ゲインを意味する(波形をg倍する）。
小数点指定。
通常は1.0（×1倍)。
マイナスを記入した場合は、波形を反転する。
このプログラムでは１.0の方が良い。
a = オフセット値の取り方を指定。（A、T、S の何れかを指定）
A と指定した場合は全サンプルの平均値をとる（tは不要）。
T と指定した場合は、先頭からt カウントの平均値を取る。
s と指定した場合は、設定した値(tカウント値)を零レベルとする。
2行目以降：filt,n,fc,h (指定チャンネルに付、32段(行)記入可能)
filt = フィルタ選択（HPF、LPF、DIF、INT、BEF、IIRの内どれか)。
HPF は、ハイパスフィルタ。
LPF は、ローパスフィルタ。
DIF は、微分。
INT は、積分。
BEF は、バンドストップフィルタ。
IIR は、IIR定数を個別に指定(注１)。
n = 次数（2、1、0、-2の内どれか)。0はフィルタなし、
DIV、INTの次数は1次のみ可能。
-2 は、2次のHPFの逆フィルタを掛ける場合に使用(補足説明の2次式参照)。
fc = カットオブ周波数（Hz)、注１。
INI、DIFを指定した場合は、ダミーとして適当な値を記入。
h = ダンピング定数値(注２)。
INI、DIFを指定した場合は、ダミーとして適当な値を記入。
注１)、IIRを指定したときは、nより後の書式が異なる。
補足説明にあるABCDE順に値を記入。
注２)、n次のバタワース特性の定数は、Butterコマンドで得ることが出来る。
Butterコマンドについては一番下にCソースを掲載した。
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
使用例：
１，つくば広帯域速度計のフィルタ無しの波形図を作成
$ W32IIRtoPS TSKFch.tbl 01220108Bm12.w32 TSKFnon.ps 0 140 300
２，つくば短周期速度計のフィルタ無しの波形図を作成
$ W32IIRtoPS TSKHch.tbl 01220108VM12.w32 TSKHnon.ps 0 140 300
３，つくば広帯域速度計の1HzHPFの波形図を作成（ 短周期速度計化）
$ W32IIRtoPS TSKFch.tbl 01220108Bm12.w32 TSKFiHzHPF.ps TSKFiHzHPF.flt 140 300
４，つくば短周期速度計の1HzHPFの逆フィルタの波形図を作成（ 広帯域速度計化）
$ W32IIRtoPS TSKHch.tbl 01220108VM12.w32 TSKHwide.ps TSKHwide.filt 140 300
左がTSKFiHzHPF.flt ファイルの中身で、右がTSKHwide.filtファイルの中身
結果：
短周期速度計化した広帯域速度計（上)と短周期速度計生波形(下)
広帯域速度計生波形（上)と広帯域速度計化した短周期速度計波形(下)
短周期化した広帯域速度計(TSKF)と短周期速度計(TSKH)生波形、広帯域速度計生波形と広帯域化した短周期の波形は似ているようであるが多少違う。TSKH200mの地下埋設でTSKFは横孔と設置場所に違うがあるので当然とも言えるが、正否の判断はよく吟味する必要があるだろう。
でもこれはあくまで、W32IIRtoPSは、こようなことも出来るとのデモンステトレーションである。
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
補足説明：
W32IIRtoPSは、指定されたフィルタ定数を下記の式に当てはめて実行する。H(s)はラプラス変換による伝達関数式、H(z)はZ変換による伝達関数式である。つまり上のラプラス変換の式は下のZ変換の式と粗等しいものとしてZ変換されて実行される。従って、完全に等価でないことを注意する必要がある。
２次式：
C･s^2 + D･s + E
H(s) = --------------------------
s^2 + A･s + B
a =2/T , G = 1/(a^2 + A*a + B), b=C*a^2
(b + D･a + E) + 2･(E - b)･z^-1 + (b - D･a + E)･z^-2
H(z) =G ＊ -----------------------------------------------------------------------
1 + 2･G･(B - a^2)･z^-1 + G･(a^2 - A･a + B)･z^-2
a2･z^-2 + a1･z^-1 + a0
H(z) =G * -------------------------------------
b2･z^-2 + b1･z^-1 + b0
１次式：
Bs + C
H(s) = -----------------
s + A
1 (B*a + C)+(C - B*a)*z^-1
H(z) = ---------- ＊------------------------------------
A + a A - a
1 + -------------* z^-1
A + a5
a0 + a1*z^-1
H(z) = G * -----------------------
b0 + b1*z^-1
a = 2/T
ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー
Butter コマンドのソース：
/********************************************************************************/
// < Butter.c > バタワース・フィルタの次数に対するダンピング定数値を出力
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.14159265358979323846264338327950288419716939937510
int main(int argc, char *argv[])
{
int n0,i,j;
double a,k,n;
if(argc==1){
printf("Butter n\n");
printf(" n : 2.. oder\n");
printf("ver.1.2: 2016/05/19\n");
exit(0);
}
n0=atoi(argv[1]);
if(n<=0){ printf("bad oder %d\n",n0); exit(0);}
n=(double)n0;
if(n0-(n0/2*2)==0){
for(i=1; i<=n0/2; i++){
k=(double)i;
a=-1.0*cos((2.0*k+n-1.0)*PI/(2.0*n));
printf("%2d : 2: %40.38G\n",i,a);
}
}else{
printf("%2d :1: \n",1);
for(i=1; i<=n0/2; i++){
k=(double)i;
a=-1.0*cos((2.0*k+n-1.0)*PI/(2.0*n));
printf("%2d :2: %40.38G\n",i+1,a);
}
}
exit(0);
}
PR
