Signal Processing Toolbox | ![]() ![]() |
Yule-Walker AR法を使ったパワースペクトル密度の計算
表示
Pxx=
pyulear(x,p) [Pxx,w]=
pyulear(x,p,nfft) [Pxx,f]=
pyulear(x,p,nfft,fs) [Pxx,f]=
pyulear(x,p,nfft,fs,'range
') [Pxx,w]=
pyulear(x,p,nfft,'range
') pyulear(...)
詳細
Pxx
は、パラメトリックスペクトル推定法の一つのYule-Walkerアルゴリズムを使って、ベクトル =
pyulear(x,p)
x
のパワースペクトル密度(PSD)、Pxx
を出力します。x
の要素は、離散時間信号のサンプルを表しています。p
は、PSDを計算するときに使用する、信号の自己回帰(AR)予測モデルの次数を指定する整数です。
パワースペクトル密度は、パワー/ラジアン/サンプルの単位で計算します。入力が実数値の場合、(周波数軸で)片側PSD(デフォルト)を、複素数の場合両側PSDを出力します。
一般的に、FFTの長さNと入力x
の値(実数/複素数)により、Pxx
の長さと対応する正規化された周波数の範囲が決まります。この書式で、FFTの(デフォルト)の長さNは、256より長くなり、x
の長さを超える2のベキ乗数の内一番小さいものになります。つぎの表は、この書式用のPxx
の長さと対応する正規化された周波数の範囲を示しています。
実数/複素数入力データ |
Pxxの長さ |
対応する正規化された周波数の範囲 |
実数値 |
129 |
[0, ![]() |
複素数値 |
256 |
[0, 2![]() |
[Pxx,w]
は、PSDを計算した周波数をベクトル =
pyulear(x,p)
w
で出力します。Pxx
とw
は、同じ長さです。周波数の単位は、ラジアン/サンプルです。
[Pxx,w]
は、Yule-Walker法を使って、PSDを計算しますが、整数 =
pyulear(x,p,nfft)
nfft
を使って、FFTの長さを設定することができます。nfft
を空ベクトル[]
とすると、デフォルト値(256)が使われます。
Pxx
の長さとw
の周波数範囲は、nfft
と入力x
の値(実数/複素数/個数)に依存します。つぎの表は、この書式でのPxx
の長さとw
の周波数範囲を示しています。
実数/複素数入力データ |
nfft 偶数/奇数 |
Pxxの長さ |
w の範囲 |
実数値 |
偶数 |
(nfft /2 + 1) |
[0, ![]() |
実数値 |
奇数 |
(nfft + 1)/2 |
[0, ![]() |
複素数値 |
偶数または奇数 |
nfft |
[0, 2![]() |
[Pxx,f]
は、PSDベクトル( =
pyulear(x,p,nfft,fs)
Pxx
)と対応する周波数ベクトル(f
)を計算するためにヘルツ単位で設定したサンプリング周波数fs
を使います。計算されるスペクトル密度は、パワー/Hzの単位です。fs
を空ベクトル[]
とすると、サンプリング周波数は、デフォルトの1 Hzです。
f
に関する周波数範囲は、nfft
とfs
と入力x
の値(実数/複素数/個数)に依存します。Pxx
の長さは、上記の表の中のものと同じです。つぎの表は、この書式を使ったf
に対する周波数範囲を示しています。
実数/複素数入力データ |
nfft 偶数/奇数 |
fの範囲 |
実数値 |
偶数 |
[0 , fs/2 ] |
実数値 |
奇数 |
[0 , fs/2) |
複素数値 |
偶数または奇数 |
[0, fs ) |
[Pxx,f]
または、[Pxx,w] =
pyulear(x,p,nfft,fs,'range
')
=
pyulear(x,p,nfft,'range
') は、周波数の範囲をf
またはw
に出力します。この書式は、x
が実数の場合に有効です。'range
' には、つぎのいずれかを設定することができます。
'twosided'
:周波数範囲[0,fs
)で両側のPSDを計算します。これは、x
が複素数の場合の周波数範囲のデフォルトの設定です。fs
を空ベクトル[]
とする場合、周波数範囲は[0,1)
です。 fs
を設定していない場合、周波数範囲は[0, 2'onesided'
:実数x
に対して設定した周波数範囲上での片側PSDを計算します。これは、x
が実数の場合の周波数範囲のデフォルトの設定です。pyulear(...)
は、出力変数を設定しない場合、カレントフィギュアウインドウにPSDをプロット表示します。プロット上の周波数範囲は、与えられたパラメータに対して出力w
(または f
)の範囲と同じです。
例題
Yule-Walker法は、設定した次数のAR予測モデルを信号に近似することで、スペクトル密度を計算するために、まず、設定した次数もAR(全極)モデルから信号を生成しましょう。関数freqz
を使って、ユーザのARフィルタの周波数応答の大きさをチェックすることができます。関数pyulear
を使ってPSDを計算する場合、結果を推定することができます。
a =
[1 -2.2137 2.9403 -2.1697 0.9606]; % AR フィルタ係数
freqz(1,a) % AR フィルタの周波数応答
title('AR System Frequency Response')
さて、白色ノイズをARフィルタに適用して入力信号x
を作成します。4次のAR予測モデルをベースにx
のPSDを計算しましょう。ここで、オリジナルのARシステムモデルは、次数が4であることがわかっているとします。
randn('state',1);
x =
filter(1,a,randn(256,1)); % AR システム出力
pyulear(x,4) % 4次の推定
注意
パワースペクトル密度は、単位周波数に対するパワーの分布として計算できます。
このアルゴリズムは、ユーザの信号に対して、選択したモデル次数によって影響を受けます。
アルゴリズム
線形予測フィルタは、信号の2次の統計的な特徴をモデル化するのに使われます。予測フィルタ出力は、入力が白色ノイズの場合、信号をモデル化するのに使われます。
pyulear
は、Yule-Walker 法を用いて信号のPSDを推定します。この方法は、最小二乗の意味で、フォーワード予測誤差を最小化することにより、自己回帰(AR)線形予測フィルタモデルを信号に近似します。この式は、Yule-Walker方程式から導かれ、Levinson-Durbin再帰法を使って、解かれます。pyulear
で計算されるスペクトル推定は、このARモデルの周波数応答の大きさを二乗したものです。
参考
|
Yule-Walker法を使ったARモデルパラメータの推定 |
|
線形予測係数 |
|
Burg法を使って、パワースペクトル密度を計算 |
|
共分散法を使ったパワースペクトル密度の計算 |
|
固有ベクトル法を使った擬似スペクトルの推定 |
|
ピリオドグラムを使ったパワースペクトル密度の計算 |
|
修正共分散法を使ったパワースペクトル密度の計算 |
|
マルチテーパ法を使ったパワースペクトル密度の計算 |
|
MUSIC アルゴリズムを使った擬似スペクトルの推定 |
|
時間領域IIRフィルタ設計に対するProny法
|
|
パワースペクトル密度データのプロット |
|
Welch法を使ったパワースペクトル密度の計算 |
参考文献
[1] Marple, S.L., Digital Spectral Analysis, Prentice-Hall, 1987, Chapter 7.
[2] Stoica, P., and R.L. Moses, Introduction to Spectral Analysis, Prentice-Hall, 1997.
![]() | pwelch | rc2ac | ![]() |