pspect
Welchの平均ペリオドグラム法による 2つの離散時間信号の間の両面相互スペクトル推定.
呼び出し手順
sm = pspect(sec_step, sec_leng, wtype, x) sm = pspect(sec_step, sec_leng, wtype, x, y) sm = pspect(sec_step, sec_leng, wtype, nx) sm = pspect(sec_step, sec_leng, wtype, nx, ny) [sm, cwp] = pspect(.., wpar)
引数
- x
ベクトル, 最初の信号の時間領域標本.
- y
ベクトル, 2番目の信号の時間領域標本.
y
が 省略された場合,x
に等しいと仮定されます(自己相関). 指定された場合,x
と同じ数の要素を有している必要があります.- nx
スカラー :
x
信号の標本数. この場合,x
信号のセグメントはgetx
という名前の ユーザ定義の関数によりロードされます (下記参照).- ny
スカラー :
y
信号の標本数. この場合,x
信号のセグメントはgety
という名前の ユーザ定義の関数によりロードされます (下記参照). 指定された場合,ny
はnx
と等しい必要があります.- sec_step
各データウインドウのオフセット. オーバーラップ
D
はisec_leng -sec_step
により指定されます.sec_step==sec_leng/2
の場合, 50% オーバーラップとなります.- sec_leng
ウインドウの点の数.
- wtype
ウインドウの種類
're'
: 矩形'tr'
: 三角形'hm'
: ハミング'hn'
: ハニング'kr'
: カイザー,この場合,wpar引数を指定する必要があります'ch'
: チェビシェフ, この場合,wpar引数を指定する必要があります
- wpar
カイザーおよびチェビシェフウインドウ
に関する オプションのパラメータ:'kr':
wpar
は厳密に正の数とする必要があります'ch':
wpar
は要素数2のベクトル[main_lobe_width,side_lobe_height]
とする必要があります. ただし,0<main_lobe_width<.5
, およびside_lobe_height>0
です.
- sm
正規化された周波数の範囲
[0,1]
における両面スペクトル推定. これは,sec_len
個の要素を有する行配列です. この配列は自己相関の場合には実数,相互相関の場合には複素数となります.関係する正規化された周波数配列は
linspace(0,1,sec_len)
です.- cwp
チェビシェフウインドの場合は チェビシェフウインドウの未指定のパラメータ, または空の行列.
説明
2つの信号x
および y
が両方共指定
された場合,相互スペクトル推定値,
そうでない場合に自己相関スペクトルを計算します.
修正ペリオドグラム法により得られたスペクトル推定値.
2つの信号x
および y
の
相互スペクトルは以下のように定義されます
スペクトル推定の修正ペリオドグラム法は,
x
および y
に含まれる
データのウインドウ処理されたサブセクションのペリオドグラムを
繰り返し計算します.
これらのペリオドグラムは,最終的なスペクトル推定値を得るために平均化され,
適当な定数により正規化されます.
これは推定値の分散を減らすための平均化処理です.
バッチ処理の場合, x
および
y
データはユーザ定義関数
getx
および gety
により
セグメント毎に読まれます.
これらの関数の呼び出し手順は以下のようになります:
xk = getx(ns,offset)
および
yk = gety(ns,offset)
. ただし,ns
は
セグメントの大きさ, offset
は
完全な信号のセグメントの最初の要素のインデックスです.
参考文献
Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing, Upper Saddle River, NJ: Prentice-Hall, 1999
例
rand('normal');rand('seed',0); x = rand(1:1024-33+1); // make low-pass filter with eqfir nf = 33; bedge = [0 .1;.125 .5]; des = [1 0]; wate = [1 1]; h =eqfir(nf, bedge, des, wate); // filter white data to obtain colored data h1 = [h 0*ones(1:max(size(x))-1)]; x1 = [x 0*ones(1:max(size(h))-1)]; hf = fft(h1,-1); xf = fft(x1,-1); y = real(fft(hf.*xf,1)); // plot magnitude of filter h2 = [h 0*ones(1:968)]; hf2 = fft(h2,-1); hf2 = real(hf2.*conj(hf2)); hsize = max(size(hf2)); fr = (1:hsize) / hsize; plot(fr, log(hf2)); // pspect example sm = pspect(100,200,'tr',y); smsize = max(size(sm)); fr = (1:smsize) / smsize; plot(fr, log(sm)); rand('unif');
参照
Report an issue | ||
<< mese | Spectral estimation | Transforms >> |