SPECTRUM1D

spectrum1d − 1 つ[又は 2 つ]の時系列からパワースペクトル[及びクロススペ クトル]を計算する

spectrum1d [ x[y]file ] −Ssegment_size] [ −C[xycnpago] ] [ −Ddt ] [ −Nname_stem ] [ −V ] [ −W ] [ −b[i|o][s|S|d|D][ncol] ] [ −f[i|o]colinfo ]

spectrum1d は X [及び Y ]の値を標準入力[又は x[y]file ]の 1 列目[及び 2 列 目] か ら 読 み込む。これらの値は等しい時間間隔 dt を持つ時系列 X(t) [Y(t)] として扱われる。入力の行数は何行でも構わない 。 spectrum1d は Bendat and Piersol の標準誤差の推定法を用いて複数の重なり合う窓のアンサ ンブル平均をとる Welch の方法によるパワースペクトル密度[及びクロスス ペ クトル密度]の推定値を含むファイルを作る。

出 力 ファイルには 3 列のデータ、 f 又は w, p, 及び e が入る。 f 又は w は周波数又は波長、 p は推定スペクトル密度、 e は 1 標準偏差でエラーバー の サ イズである。これらのファイルは name_stem に基づいて名付けられる。 −C オプションが使われたときは、最大 8 つのファイルが作られる。そうで な いときは 1 つ( xpower )だけが書き出される。ファイル( −bo が設定されない 限りアスキー)は以下のようになる。

name_stem.xpower

X(t) のパワースペクトル密度。 X * X * dt の単位。

name_stem.ypower

Y(t) のパワースペクトル密度。 Y * Y * dt の単位。

name_stem.cpower

コヒーレント出力のパワースペクトル密度。 ypower と同じ単位。

name_stem.npower

ノイズ出力のパワースペクトル密度。 ypower と同じ単位。

name_stem.gain

ゲイン出力、すなわち伝達関数の絶対値。 (Y / X) の単位。

name_stem.phase

位相スペクトル、すなわち伝達関数の位相。単位はラジアン。

name_stem.admit

アドミッタンススペクトル、すなわち伝達関数の実部。 (Y / X) の単 位。

name_stem.coh

(2乗)コヒーレンシースペクトル、すなわち周波数の関数としての線形 相関係数。 [0, 1] の範囲の無次元数。信号雑音比 (SNR) は coh / (1 - coh) である。 coh = 0.5 のとき SNR = 1 となる。

x[y]file

最 初の 1 [又は 2 ]列に X(t) [Y(t)] のサンプルを持つアスキー(又 はバイナリ、 −bi 参照)ファイル。ファイルを指定しないと き は 、 spectrum1d は標準入力から読み込む。

−S

segment_size はアンサンブル平均の窓に含まれるサンプルの数で 2 の冪乗の数である。推定される最小の周波数は 1.0/(segment_size * dt) で、最大は 1.0/(2 * dt) である。パワースペクトル密度の標 準偏差はおよそ 1.0 / sqrt(n_data / segment_size) になるので、も し segment_size = 256 なら、 10% の標準エラーバーを得るためには 25,600 のデータが必要になる。クロススペクトルのエラーバーはより 大きくより複雑で、コヒーレンシーの関数にもなっている。

−C

入力の最初の 2 列を 2 つの時系列 X(t) と Y(t) のサンプルと して読み込む。ノイズのある線形システムにおいて Y(t) が 出 力 で X(t) が入力であると考える。ノイズ出力が最小になりコヒーレント出 力とノイズ出力が相関を持たないような、最適な周波数応答関数を 最 小自乗法で推定する。随意で任意の順序で文字をセット{ x y c n p a g o }から最大 8 つまで指定し、デフォルト[全て]の代わりにこれ ら の 出力ファイルだけを作ることができる。 x = xpower, y = ypower, c = cpower, n = npower, p = phase, a = admit, g = gain, o = coh.

−D

dt 時系列のサンプルの間隔を設定する[デフォルト = 1]。

−N

name_stem 出力ファイルに使われる名前の主部を与える[デフォル ト = "spectrum" ]。

−V

冗長モードを選択する。標準エラー出力に経過報告を送る[デフォ ルトでは「黙って」実行する]。

−W

出力ファイルの 1 列目に周波数でなく波長を書き出す [デフォル ト = 周波数 (周期 / dt )]。

−bi

バイナリ入力を選択する。単精度に対しては s を付けること [ デ フォルトは d (倍精度)]。大文字の S (又は D) にするとバイトス ワップを行う。随意で、バイナリファイルの列数がこのプログラム で 必 要な列数を越えるときは、バイナリファイルの列数 ncol を付ける こと [デフォルトは 2 入力列]。

−bo

バイナリ出力を選択する。単精度に対しては s を付けること [ デ フォルトは d (倍精度)]。大文字の S (又は D) にするとバイトス ワップを行う。随意で、バイナリ出力ファイルの必要な列数 ncol を 付けること [デフォルトは 2 出力列]。

−f

入 出 力列の特殊フォーマットを使う(時間又は地理データ)。 i(nput) 又は o(utput) を指定すること [デフォルトでは入力と出 力 の 両方]。コンマ区切りで 1 列以上(又は列の範囲)を与えること。各 列又は列の範囲の項目に対して T (絶対暦時刻)、 t (選 択 さ れ た TIME_EPOCH に対する相対時間)、 x (経度)、 y (緯度)、 f (小数) のいずれかを付けること。 −f[i|o]g は −f[i|o]0x,1y (地理座標) の 省略形である。

数値データのアスキー出力形式は .gmtdefaults4 ファイルのパラメータによっ て決まる。緯度経度は OUTPUT_DEGREE_FORMAT に従ってフォーマットされ、 一 方 、他の値は D_FORMAT に従ってフォーマットされる。フォーマットの結果と して出力の精度が落ちて、その先の処理に影響する可能性があることに注意 す ること。十分な精度で出力されていないと思ったら、バイナリ出力に変更する( もし使えるなら −bo とする)か、 D_FORMAT の設定を用いてより多くの桁数 を 指定すること。

data.g は mGal 単位の重力データで、 1.5 km 間隔にサンプルされているもの とする。このパワースペクトルを単位 mGal**2-km でファイル data.xpower に 書き出す。

spectrum1d data.g −S256 −D1.5 −Ndata

data.g に 加えて data.g と同じ位置でサンプルされた m 単位の地形データ data.t があるとする。 data.t を入力、 data.g を出力と考えて、その伝達関 数の様々な特徴を推定する。

paste data.t data.g | spectrum1d −S256 −D1.5 −Ndata −C

GMT(l), grdfft(l)

Bendat, J. S., and A. G. Piersol, 1986, Random Data, 2nd revised ed., John Wiley & Sons.
Welch, P. D., 1967, The use of Fast Fourier Transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms, IEEE Transactions on Audio and Electroacoustics, Vol AU-15, No 2.