Signal Processing Toolbox    
pmtm

マルチテーパ法(MTM)を使ったパワースペクトル推定

表示

詳細

pmtm は、参考文献[1]のマルチテーパ法(MTM)を使って、時系列xのパワースペクトル密度(PSD)を計算します。この方法は、PSDの計算に修正したピリオドグラムの線形結合または非線形結合のいずれかを使います。これらのピリオドグラムは、離散扁長回転楕円体列(関数dpssを参照)から設定される直交テーパ(周波数領域の中のウインドウ)列を使って計算されます。

[Pxx,w] = pmtm(x,nw) は、マルチテーパ推定法の中のデータテーパとして使われる 2nw-1個の離散扁長回転楕円体列を使って、入力信号xに対するPSD Pxxを計算します。nwは、離散扁長回転楕円体列に対する時間- 帯域幅の積です。nwを空ベクトル[]に設定すると、デフォルトの4を使います。他の取り得る値は、2, 5/2, 3, 7/2のいずれかです。関数pmtm は、PSDを計算する周波数点のベクトルwを出力します。Pxxwは、同じ長さです。周波数の単位は、ラジアン/サンプルです。

パワースペクトル密度は、パワー/ラジアン/サンプルの単位で計算されます。xが実数値の場合、片側PSDが計算され、複素数の場合、両側PSDが計算されます。

一般的に、FFTの長さNと入力xの値(実数/複素数)により、Pxxの長さと対応する正規化された周波数の範囲が決まります。この書式で、FFTの(デフォルト)の長さNは、256より長くなり、xの長さを超える2のベキ乗数の内一番小さいものになります。つぎの表は、この書式用のPxxの長さと対応する正規化された周波数の範囲を示しています。

表 7-26:  長さNのFFTのPSD ベクトル特性
実数/複素数入力データ
Pxxの長さ
対応する正規化された周波数の範囲
実数値
(N/2) +1
[0, ]
複素数値
N
[0, 2)

[Pxx,w] = pmtm(x,nw,nfft) は、FFTの長さに整数nfftを使って、PSDをマルチテーパ法を使って、計算します。nfftに空ベクトル[]を設定する場合、前述の書式で記述されたNに対するデフォルト値を使います。

Pxxの長さとwに関する周波数範囲は、nfftと入力xの値(実数/複素数/個数)に依存します。つぎの表は、この書式に対するPxxの長さとwに関する周波数の範囲を示しています。

表 7-27:  PSD と周波数ベクトル特性
実数/複素数入力データ
nfft 偶数/奇数
Pxxの長さ
wの範囲
実数値
偶数
(nfft/2 + 1)
[0, ]
実数値
奇数
(nfft + 1)/2
[0, )
複素数値
偶数または奇数
nfft
[0, 2)

[Pxx,f] = pmtm(x,nw,nfft,fs) は、PSDベクトル(Pxx)と対応する周波数ベクトル(f)を計算するためにヘルツ単位で設定したサンプリング周波数fsを使います。計算されるスペクトル密度は、パワー/Hzの単位です。fsを空ベクトル[]とすると、サンプリング周波数は、デフォルトの1 Hzです。

fに関する周波数範囲は、nfftfsと入力xの値(実数/複素数/個数)に依存します。Pxxの長さは、上記の表の中のものと同じです。つぎの表は、この書式を使ったfに対する周波数範囲を示しています。

表 7-28:  fsを指定して、
PSD と周波数ベクトルの特性
実数/複素数入力データ
nfft 偶数/奇数
fの範囲
実数値
偶数
[0fs/2]
実数値
奇数
[0fs/2)
複素数値
偶数または奇数
[0, fs)

[Pxx,Pxxc,f] = pmtm(x,nw,nfft,fs) は、Pxxの95%の信頼区間 Pxxcを出力します。信頼区間は、カイ二乗アプローチを使って計算されます。Pxxcは、Pxxと同じ行数をもつ、2列からなる行列です。Pxxc(:,1)は、信頼区間の下限、Pxxc(:,2)は上限です。

[Pxx,Pxxc,f] = pmtm(x,nw,nfft,fs,p) は、Pxxp*100%の信頼区間Pxxcを出力します。ここで、pは0と1の間のスカラです。pを設定しない場合、または、空ベクトル[]の場合、デフォルトの95%が使われます。

[Pxx,Pxxc,f] = pmtm(x,e,v,nfft,fs,p) は、行列eの列に含まれるデータテーパとベクトルvの集中度から、PSDの計算値Pxxと信頼区間Pxxc、周波数ベクトルfを出力します。vの長さは、eの中の列数と同じです。これらの引数を与えて、関数dpssの出力からこれらのデータを得ることもできます。

[Pxx,Pxxc,f] = pmtm(x,dpss_params,nfft,fs,p) は、関数dpssへの入力引数を含むセル配列dpss_paramsを使って、データテーパを計算することができます。たとえば、pmtm(x,{3.5,'trace'},512,1000)は、nfft = 512, fs = 1000 を使って、nw = 3.5に対するSlepian列を計算し、この計算に使った関数dpssの方法を表示します。他のオプションについては、関数dpssを参照してください。

[...] = pmtm(...,'method') は、個々のスペクトル推定を組み合わせるときに使用するアルゴリズムを設定します。文字列'method'は、つぎのいずれかを使用することができます。

[...] = pmtm(...,'range') は、周波数の範囲をfまたはwに出力します。この書式は、xが実数の場合に有効です。'range'には、つぎのいずれかを設定することができます。

出力引数をもたない pmtm(...) は、カレントフィギュアウインドウにPSDと信頼区間をプロットします。fsを設定していない場合、95%の信頼区間がプロットされます。fsを設定すると、pの値に依存した信頼区間がプロットされます。

例題

この例題は、白色ノイズを含んだ正弦波を解析します。

参考
dpss
離散扁長回転楕円体列の計算
pburg
Burg法を使ったパワースペクトル密度を計算
pcov
共分散法を使ったパワースペクトル密度の計算
peig
固有ベクトル法を使った擬似スペクトルの推定
periodogram
ピリオドグラムを使ったパワースペクトル密度の計算
pmcov
修正共分散法を使ったパワースペクトル密度の計算
pmusic
MUSIC アルゴリズムを使った擬似スペクトルの推定
psdplot
パワースペクトル密度データのプロット
pwelch
Welch法を使ったパワースペクトル密度の計算
pyulear
Yule-Walker AR法を使ったパワースペクトル密度の計算

参考文献

[1] Percival, D.B., and A.T. Walden, Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques, Cambridge University Press, 1993.

[2] Thomson, D.J., "Spectrum estimation and harmonic analysis," Proceedings of the IEEE, Vol. 70 (1982), pp. 1055-1096.


 pmcov pmusic