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