pwelch (Signal Processing Toolbox)

Welch法を使った信号のパワースペクトル密度(PSD)の推定

表示

[Pxx,w] = pwelch(x)
[Pxx,w] = pwelch(x,nwin)
[Pxx,w] = pwelch(x,nwin,noverlap)
[Pxx,w] = pwelch(x,nwin,noverlap,nfft)
[Pxx,f] = pwelch(x,nwin,noverlap,nfft,fs)
[...] = pwelch(x,nwin,noverlap,...,'range')
pwelch(...)

詳細

[Pxx,w] = pwelch(x) は、スペクトル推定のWelchの平均化修正ピロオドグラム法を使って、入力信号ベクトルxのパワースペクトル密度Pxxを計算します。つぎの書式を使います。

ベクトルxは、同じ長さの8つの部分(お互い50%の重なりをもつ)に分割します。

上の8つの分割を行った残りの部分は削除します。

各部分には、その部分と同じ長さをもつHammingウインドウ(関数 hammingを参照)を適用します。

パワースペクトル密度は、パワー/ラジアン/サンプル単位で計算されます。対応する周波数ベクトルwは、ラジアン/サンプルで計算され、Pxxと同じ長さをもっています。

xが実数値入力ベクトルの場合、片側PSDが計算され、複素数の場合、両側PSDが計算されます。

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

表 7-34: 長さNのFFTのPSD ベクトル特性
実数/複素数入力データ
Pxxの長さ
対応する正規化された周波数の範囲

実数値
(N/2) +1
[0, ]

複素数値
N
[0, 2)

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

[Pxx,w] = pwelch(x,nwin) は、つぎのいずれかを使って、修正ピリオドグラムを計算します。

nwinが正の整数の場合、nwinはHammingウインドウの長さ

nwinがベクトルの場合、ウインドウに重みも設定できます。

この書式では、入力ベクトルxは、お互い50%の重なりをもちながら、ある整数倍で、同じ長さをもつように分割し、残りの部分は、削除します。nwinに空ベクトル[]を設定すると、信号データは、8つに分割され、その各々にHammingウインドウが適用されます。

[Pxx,w] = pwelch(x,nwin,noverlap) は、nwinに従ってxを分割し、整数noverlapを使って、隣同士の分割で共有する部分の信号のサンプル数を設定します。引数noverlapは、ユーザが設定するウインドウの長さよりも短い必要があります。noverlapを空ベクトル[]に設定すると、pwelchは、デフォルトの重ね合わせが50%になるように、xの分割を行います。

[Pxx,w] = pwelch(x,nwin,noverlap,nfft) は、Welch法を使って、PSDを計算しますが、FFTの長さを整数nfftで設定することができます。nfftを空ベクトル[]で設定すると、上記の書式でリストされているNに対するデフォルト値を使います。

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

表 7-35: PSD と周波数ベクトルの特性
実数/複素数入力データ
nfft 偶数/奇数
Pxxの長さ
w の範囲

実数値
偶数
(nfft/2 + 1)
[0, ]

実数値
奇数
(nfft + 1)/2
[0, )

複素数値
偶数または奇数
nfft
[0, 2)

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

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

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

表 7-36: fsを指定して、PSD と
周波数ベクトルの特性
実数/複素数入力データ
nfft 偶数/奇数
fの範囲

実数値
偶数
[0, fs/2]

実数値
奇数
[0, fs/2)

複素数値
偶数または奇数
[0, fs)

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

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

'twosided': 周波数範囲[0,fs)で両側PSDを計算します。これは、xが複素数の場合の周波数範囲のデフォルトの設定です。

fsを空ベクトル[]とする場合、周波数範囲は[0,1)です。

fsを設定していない場合、周波数範囲は[0, 2)です。

'onesided': 実数xに対して設定した周波数範囲上での片側PSDを計算します。これは、xが実数の場合の周波数範囲のデフォルトの設定です。

文字引数'range' は、引数noverlapの後の場合、任意の位置に設定することができます。

出力変数を設定しない pwelch(x,...) は、カレントフィギュアウインドウに周波数に対してdB表示のPSDをプロット表示します。

例題

正弦波とノイズから構成される信号のPSDを計算しましょう。サンプリング周波数は、1000 Hzです。32のサンプルのオーバラップをもった長さ33サンプルのウインドウとデフォルトのFFTの長さを使って、両側PSD計算結果を表示します。

randn('state',0);
Fs = 1000;   t = 0:1/Fs:.3;
x = cos(2*pi*t*200)+randn(size(t));  % 200Hzの正弦波とノイズから構成
pwelch(x,33,32,[],Fs,'twosided')

アルゴリズム

pwelchは、Welch法を使って、パワースペクトル密度を計算します。

入力信号ベクトル x は、nwinとnoverlap(またはデフォルト値)に従って、k個のオーバラップした部分に分割します。
設定した(または、デフォルト)ウインドウは、xを分割したものの各々に適用されます。
nfft点のFFTの計算は、ウインドウを適用したデータに対して実行されます。
各ウインドウを適用された部分の(修正)ピリオドグラムが計算されます(関数periodogramのリファレンスの式 7-22を参照してください)。
修正されたピリオドグラムすべてから平均を計算し、スペクトル推定S(ej)を作ります。
結果求まるスペクトルは、次式のようなパワースペクトル密度を計算するためにスケーリングされます。

、ここで、Fは、
- サンプリング周波数を設定していない場合、2
- サンプリング周波数を設定している場合、fsです。

xを分割する数k は、つぎのように計算されます。

nwinを設定していないか、または、nwinを空ベクトル[]とすると、8です。

nwinを空でないベクトルまたはスカラで設定すると、です。

この方程式の中で、mは信号ベクトルxの長さで、oはオーバラップのサンプル数(noverlap)で、lは各部分の長さ(ウインドウの長さ)です。

参考

pburg
Burg法を使って、パワースペクトル密度を計算

pcov
共分散法を使ったパワースペクトル密度の計算

peig
固有ベクトル法を使った擬似スペクトルの推定

periodogram
ピリオドグラムを使ったパワースペクトル密度の計算

pmcov
修正共分散法を使ったパワースペクトル密度の計算

pmtm
マルチテーパ法を使ったパワースペクトル密度の計算

pmusic
MUSIC アルゴリズムを使った擬似スペクトルの推定

psdplot
パワースペクトル密度データのプロット

pyulear
Yule-Walker AR法を使ったパワースペクトル密度の計算

`pburg`	Burg法を使って、パワースペクトル密度を計算
`pcov`	共分散法を使ったパワースペクトル密度の計算
`peig`	固有ベクトル法を使った擬似スペクトルの推定
`periodogram`	ピリオドグラムを使ったパワースペクトル密度の計算
`pmcov`	修正共分散法を使ったパワースペクトル密度の計算
`pmtm`	マルチテーパ法を使ったパワースペクトル密度の計算
`pmusic`	MUSIC アルゴリズムを使った擬似スペクトルの推定
`psdplot`	パワースペクトル密度データのプロット
`pyulear`	Yule-Walker AR法を使ったパワースペクトル密度の計算

参考文献

[1] Hayes, M., Statistical Digital Signal Processing and Modeling, John Wiley & Sons, 1996.

[2] Stoica, P., and R.L. Moses, Introduction to Spectral Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1997, pp. 52-54.

[3] Welch, P.D, "The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms," IEEE Trans. Audio Electroacoustics, Vol. AU-15 (June 1967), pp. 70-73.

pulstran pyulear