ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

pwelch

ウェルチのパワー スペクトル密度推定

構文

pxx = pwelch(x)
pxx = pwelch(x,window)
pxx = pwelch(x,window,noverlap)
pxx = pwelch(x,window,noverlap,nfft)
[pxx,w] = pwelch(___)
[pxx,f] = pwelch(___,fs)
[pxx,w] = pwelch(x,window,noverlap,w)
[pxx,f] = pwelch(x,window,noverlap,f,fs)
[___] = pwelch(x,window,___,freqrange)
[___] = pwelch(x,window,___,spectrumtype)
[___] = pwelch(x,window,___,trace)
[___,pxxc] = pwelch(___,'ConfidenceLevel',probability)
pwelch(___)

説明

pxx = pwelch(x) では、ウェルチのオーバーラップ セグメント平均推定器を使用して検出された入力信号 x のパワー スペクトル密度 (PSD) 推定 pxx が返されます。x がベクトルの場合、単一チャネルとして取り扱われます。x が行列の場合、PSD は列ごとに個別に計算され、pxx の対応する列に保存されます。x が実数値の場合、pxx は片側 PSD 推定です。x が複素数値の場合、pxx は両側 PSD 推定です。既定では、x は可能な限り長いセグメントに分割され、8 にできるだけ近い (ただし 8 を超えない) 数のセグメントが 50% のオーバーラップで取得されます。各セグメントにはハミング ウィンドウが適用されます。修正ピリオドグラムは PSD 推定を求めるために平均化されます。x の長さを 50% のオーバーラップをもたせて整数値数のセグメントに分割できない場合、x はそれに応じた長さで打ち切られます。

pxx = pwelch(x,window) では、入力ベクトルまたは整数 window を使用して信号がセグメントに分割されます。window がベクトルの場合、pwelch は信号をその長さが window の長さに等しいセグメントに分割します。修正ピリオドグラムは、ベクトル window で乗算された信号のセグメントを使用して計算されます。window が整数の場合、信号は window の長さのセグメントに分割されます。修正ピリオドグラムは長さ window のハミング ウィンドウを使用して計算されます。

pxx = pwelch(x,window,noverlap) は、セグメント間で noverlap 個のサンプルのオーバーラップを使用します。window が整数の場合、noverlapwindow より小さい正の整数でなければなりません。window がベクトルの場合、noverlapwindow の長さより小さい正の整数でなければなりません。noverlap を指定しない場合、または noverlap を空に指定する場合、オーバーラップするサンプル数の既定値はウィンドウの長さの 50% です。

pxx = pwelch(x,window,noverlap,nfft) では、PSD 推定で使用する離散フーリエ変換 (DFT) 点の数が指定されます。既定の nfft は 256 を超える数、またはセグメント長を超える最小の 2 のべき乗です。

[pxx,w] = pwelch(___) では、正規化周波数ベクトル w が返されます。pxx が片側 PSD 推定の場合、w は偶数の nfft に対しては区間 [0,π] をカバーし、奇数の nfft に対しては [0,π) をカバーします。pxx が両側 PSD 推定の場合、w は区間 [0,2π) をカバーします。

[pxx,f] = pwelch(___,fs) は、周波数ベクトル f を単位時間あたりのサイクル数で返します。サンプルレート fs は単位時間あたりのサンプル数です。時間の単位が秒の場合、f の単位はサイクル/秒 (Hz) です。実数値の信号の場合、f は偶数の nfft に対しては区間 [0,fs/2] をカバーし、奇数の nfft に対しては [0,fs/2) をカバーします。複素数値の信号 f は区間 [0,fs) をカバーします。fspwelch の 5 番目の入力でなければなりません。サンプルレートを入力した場合でも、前のオプション引数の既定値を使用するには、これらの引数を空 [] として指定します。

[pxx,w] = pwelch(x,window,noverlap,w) では、ベクトル w で指定される正規化周波数での両側ウェルチ PSD 推定が返されます。ベクトル w は少なくとも 2 つの要素をもたなければなりません。

[pxx,f] = pwelch(x,window,noverlap,f,fs) では、ベクトル f で指定される周波数での両側ウェルチ PSD 推定が返されます。ベクトル f は少なくとも 2 つの要素をもたなければなりません。f の周波数の単位は単位時間あたりのサイクルです。サンプルレート fs は単位時間あたりのサンプル数です。時間の単位が秒の場合、f の単位はサイクル/秒 (Hz) です。

[___] = pwelch(x,window,___,freqrange) では、freqrange で指定される周波数範囲にわたるウェルチ PSD 推定が返されます。freqrange に指定できる有効なオプションは、'onesided''twosided' または 'centered' です。

[___] = pwelch(x,window,___,spectrumtype) では、spectrumtype'psd' に指定されている場合は PSD 推定が返され、spectrumtype'power' に指定されている場合はパワー スペクトルが返されます。

[___] = pwelch(x,window,___,trace) は、trace'maxhold' に指定されている場合は最大ホールド スペクトル推定を返し、trace'minhold' に指定されている場合は最小ホールド スペクトル推定を返します。

[___,pxxc] = pwelch(___,'ConfidenceLevel',probability) では、PSD 推定の probability × 100% 信頼区間が pxxc で返されます。

出力引数を設定せずに pwelch(___) を使用すると、現在の Figure ウィンドウにウェルチ PSD 推定がプロットされます。

すべて折りたたむ

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの離散時間正弦波で構成される入力信号のウェルチ PSD 推定を求めます。

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの正弦波を作成します。再現可能な結果が必要な場合は、乱数発生器をリセットします。信号長は 320 サンプルです。

rng default

n = 0:319;
x = cos(pi/4*n)+randn(size(n));

既定のハミング ウィンドウと DFT 長を使用してウェルチ PSD 推定を求めます。既定のセグメント長は 71 サンプル、DFT 長は 256 点で、周波数分解能は ラジアン/サンプルとなっています。信号は実数値のため、ピリオドグラムは片側で、256/2+1 の点が含まれます。

pxx = pwelch(x);
plot(10*log10(pxx))

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの離散時間正弦波で構成される入力信号のウェルチ PSD 推定を求めます。

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの正弦波を作成します。再現可能な結果が必要な場合は、乱数発生器をリセットします。信号のサンプル数は 512 です。

rng default

n = 0:511;
x = cos(pi/3*n)+randn(size(n));

信号を 132 サンプルの長さのセグメントに分割するウェルチ PSD 推定を求めます。信号のセグメントは 132 サンプルの長さのハミング ウィンドウで乗算されます。オーバーラップされたサンプルの数は指定されていないため、132/2 = 66 に設定されます。DFT 長は 256 点で、周波数分解能は ラジアン/サンプルとなっています。信号は実数値のため、PSD 推定は片側で、256/2+1 = 129 の点が含まれます。正規化周波数の関数として PSD をプロットします。

segmentLength = 132;
[pxx,w] = pwelch(x,segmentLength);

plot(w/pi,10*log10(pxx))
xlabel('\omega / \pi')

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの離散時間正弦波で構成される入力信号のウェルチ PSD 推定を求めます。

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの正弦波を作成します。再現可能な結果が必要な場合は、乱数発生器をリセットします。信号長は 320 サンプルです。

rng default

n = 0:319;
x = cos(pi/4*n)+randn(size(n));

信号を 100 サンプルの長さのセグメントに分割するウェルチ PSD 推定を求めます。信号のセグメントは 100 サンプルの長さのハミング ウィンドウで乗算されます。オーバーラップされたサンプルの数は 25 です。DFT 長は 256 点で、周波数分解能は ラジアン/サンプルとなっています。信号は実数値のため、PSD 推定は片側で、256/2+1 の点が含まれます。

segmentLength = 100;
noverlap = 25;
pxx = pwelch(x,segmentLength,noverlap);

plot(10*log10(pxx))

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの離散時間正弦波で構成される入力信号のウェルチ PSD 推定を求めます。

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの正弦波を作成します。再現可能な結果が必要な場合は、乱数発生器をリセットします。信号長は 320 サンプルです。

rng default

n = 0:319;
x = cos(pi/4*n) + randn(size(n));

信号を 100 サンプルの長さのセグメントに分割するウェルチ PSD 推定を求めます。50% の既定オーバーラップを使用します。DFT 長が 640 点となるように指定し、 ラジアン/サンプルの周波数が DFT ビン (ビン 81) に対応するようにします。信号は実数値のため、PSD 推定は片側で、640/2+1 の点が含まれます。

segmentLength = 100;
nfft = 640;
pxx = pwelch(x,segmentLength,[],nfft);

plot(10*log10(pxx))
xlabel('rad/sample')
ylabel('dB / (rad/sample)')

N(0,1) 加法性ホワイト ノイズを伴う 100 Hz の正弦波で構成される信号を生成します。再現可能な結果が必要な場合は、乱数発生器をリセットします。サンプルレートは 1 kHz で、信号の持続時間は 5 秒です。

rng default

fs = 1000;
t = 0:1/fs:5-1/fs;
x = cos(2*pi*100*t) + randn(size(t));

前の信号についてウェルチのオーバーラップ セグメント平均 PSD 推定を求めます。300 個のサンプルがオーバーラップする 500 サンプル長のセグメント長を使用します。500 DFT 点を使用して 100 Hz が直接 DFT ビンにあたるようにします。サンプルレートを入力して周波数ベクトルを Hz で出力します。結果をプロットします。

[pxx,f] = pwelch(x,500,300,500,fs);

plot(f,10*log10(pxx))

xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')

加法性ホワイト ノイズを伴う 100 Hz の正弦波で構成される信号を生成します。再現可能な結果が必要な場合は、乱数発生器をリセットします。サンプルレートは 1 kHz で、信号の持続時間は 5 秒です。

rng default

fs = 1000;
t = 0:1/fs:5-1/fs;

noisevar = 1/4;
x = cos(2*pi*100*t)+sqrt(noisevar)*randn(size(t));

ウェルチ法を使用して DC を中央に揃えたパワー スペクトルを求めます。300 個のサンプルがオーバーラップし、DFT 長が 500 点の 500 サンプル長のセグメント長を使用します。結果をプロットします。

[pxx,f] = pwelch(x,500,300,500,fs,'centered','power');

plot(f,10*log10(pxx))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
grid

-100 Hz と 100 Hz におけるパワーが、振幅 1 の実数値正弦波に対し 1/4 の、期待したパワーに近いことがわかります。1/4 からの偏差は加法性ノイズの影響によるものです。

200 kHz で 0.1 秒間サンプリングされた、ノイズを含む 3 つの正弦波とチャープで構成される信号を作成します。正弦波の周波数は 1 kHz、10 kHz および 20 kHz です。正弦波の振幅とノイズ レベルはさまざまです。ノイズのないチャープの周波数は 20 kHz で始まり、サンプリング中に 30 kHz まで線形に増加します。

Fs = 200e3; 
Fc = [1 10 20]'*1e3; 
Ns = 0.1*Fs;

t = (0:Ns-1)/Fs;
x = [1 1/10 10]*sin(2*pi*Fc*t)+[1/200 1/2000 1/20]*randn(3,Ns);
x = x+chirp(t,20e3,t(end),30e3);

信号のウェルチ PSD 推定と最大ホールドおよび最小ホールド スペクトルを計算します。結果をプロットします。

[pxx,f] = pwelch(x,[],[],[],Fs);
pmax = pwelch(x,[],[],[],Fs,'maxhold');
pmin = pwelch(x,[],[],[],Fs,'minhold');

plot(f,pow2db(pxx))
hold on
plot(f,pow2db([pmax pmin]),':')
hold off
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
legend('pwelch','maxhold','minhold')

この手順を繰り返し、今度はパワー スペクトル推定が中心になるように計算します。

[pxx,f] = pwelch(x,[],[],[],Fs,'centered','power');
pmax = pwelch(x,[],[],[],Fs,'maxhold','centered','power');
pmin = pwelch(x,[],[],[],Fs,'minhold','centered','power');

plot(f,pow2db(pxx))
hold on
plot(f,pow2db([pmax pmin]),':')
hold off
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
legend('pwelch','maxhold','minhold')

この例は、ウェルチのオーバーラップ セグメント平均 (WOSA) PSD 推定での信頼限界の使い方を示します。統計的有意性の必要条件ではありませんが、周囲の PSD 推定の信頼下限が信頼上限を超えるウェルチ推定の周波数は、時系列で大きな振幅を明確に示します。

N(0,1) 加法性ホワイト ノイズを伴う 100 Hz と 150 Hz の正弦波を重ね合わせて構成される信号を作成します。2 つの正弦波の振幅は 1 で、サンプルレートは 1 kHz です。再現可能な結果が必要な場合は、乱数発生器をリセットします。

rng default
fs = 1000;
t = 0:1/fs:1-1/fs;
x = cos(2*pi*100*t)+sin(2*pi*150*t)+randn(size(t));

95% 信頼限界の WOSA 推定を求めます。セグメント長を 200 に設定し、セグメントの 50% (100 サンプル) をオーバーラップさせます。信頼区間と共に WOSA PSD 推定をプロットし、100 Hz と 150 Hz 付近の周波数関心領域を拡大します。

L = 200;
noverlap = 100;
[pxx,f,pxxc] = pwelch(x,hamming(L),noverlap,200,fs,...
    'ConfidenceLevel',0.95);

plot(f,10*log10(pxx))
hold on
plot(f,10*log10(pxxc),'-.')
hold off

xlim([25 250])
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
title('Welch Estimate with 95%-Confidence Bounds')

100 および 150 Hz のごく近傍における信頼限界の下限は、100 および 150 Hz の近傍外における信頼限界の上限より大幅に上になります。

加法性ホワイト ガウス ノイズを伴う 3 つの正弦波から構成されるマルチチャネル信号の 1024 サンプルを生成します。正弦波の周波数は、 および ラジアン/サンプルです。ウェルチ法を使用して信号の PSD を推定し、プロットします。

N = 1024;
n = 0:N-1;

w = pi./[2;3;4];
x = cos(w*n)' + randn(length(n),3);

pwelch(x)

入力引数

すべて折りたたむ

行または列のベクトル、または行列として指定する入力信号。x が行列の場合、その各列は独立チャネルとして扱われます。

例: cos(pi/4*(0:159))+randn(1,160) は単一チャネルの行ベクトル信号です。

例: cos(pi./[4;2]*(0:159))'+randn(160,2) は 2 チャネル信号です。

データ型: single | double
複素数のサポート: あり

行ベクトル、列ベクトルまたは整数として指定するウィンドウ。window がベクトルの場合、pwelchx を長さが window のオーバーラップのあるセグメントに分割した後、window で指定されるベクトルで各信号セグメントを乗算します。window が整数の場合、pwelch はその整数値の長さのセグメントに分割され、同じ長さのハミング ウィンドウが使用されます。x の長さが noverlap 個のオーバーラップ サンプルをもつ整数個のセグメントに厳密に分割できない場合、x はそれに応じた長さで打ち切られます。window が空に指定されると、既定のハミング ウィンドウにより noverlap 個のオーバーラップ サンプルをもつ x の 8 個のセグメントが求められます。

データ型: single | double

window の長さより小さな正の整数として指定する、オーバーラップするサンプル数。noverlap を省略したり、noverlap を空に指定する場合、セグメント同士が 50% のオーバーラップになるような値が使用されます。

正の整数として指定する DFT 点の数。実数値の入力信号 x では、PSD 推定 pxx は、nfft が偶数の場合、長さ (nfft/2 + 1)、nfft が奇数の場合、長さ (nfft + 1)/2 になります。複素数値の入力信号 x では、PSD 推定は常に長さ nfft になります。nfft が空として指定されている場合、既定の nfft が使用されます。

nfft がセグメント長よりも長い場合、データにゼロが付加されます。nfft がセグメント長よりも短い場合、セグメントは長さが nfft と等しくなるように datawrap によってラップされます。

データ型: single | double

サンプルレート。正のスカラーで指定します。サンプルレートは単位時間あたりのサンプル数です。時間の単位が秒の場合、サンプルレートの単位は Hz です。

正規化周波数。少なくとも 2 つの要素をもつ行ベクトルまたは列ベクトルとして指定します。正規化周波数の単位はラジアン/サンプルです。

例: w = [pi/4 pi/2]

データ型: double

周波数。少なくとも 2 つの要素をもつ行ベクトルまたは列ベクトルとして指定します。周波数の単位は単位時間あたりのサイクルです。単位時間はサンプルレート fs で指定されます。fs の単位がサンプル/秒であれば、f の単位は Hz です。

例: fs = 1000; f = [100 200]

データ型: double

'onesided''twosided' または 'centered' で指定する、PSD 推定の周波数範囲。既定値は、実数値信号の場合は 'onesided'、複素数値信号の場合は 'twosided' です。各オプションに対応する周波数範囲は次のとおりです。

  • 'onesided' — 実数値入力信号 x の片側 PSD を返します。nfft が偶数の場合、pxx の長さは nfft/2 + 1 で、計算区間は [0,π] ラジアン/サンプルです。nfft が奇数の場合、pxx の長さは (nfft + 1)/2、区間は [0,π) ラジアン/サンプルです。fs がオプションで指定されると、対応する区間はそれぞれ、nfft が偶数の場合は [0,fs/2] サイクル/単位時間、奇数の場合は [0,fs/2) サイクル/単位時間になります。

  • 'twosided' — 実数値または複素数値の入力 x について両側 PSD 推定を返します。この場合、pxx の長さは nfft、計算区間は [0,2π) ラジアン/サンプルです。fs がオプションで指定される場合、その区間は [0,fs) サイクル/単位時間になります。

  • 'centered' — 実数値または複素数値の入力 x について中央に揃えた両側 PSD 推定を返します。この場合、pxx の長さは nfft、偶数長の nfft については区間 (–π,π] ラジアン/サンプルで、奇数長の nfft については (–π,π) ラジアン/サンプルで計算されます。fs がオプションで指定されると、対応する区間はそれぞれ、nfft が偶数長の場合は (–fs/2,fs/2] サイクル/単位時間、奇数長の場合は (–fs/2,fs/2) サイクル/単位時間になります。

'psd' または 'power' で指定するパワー スペクトルのスケーリング。spectrumtype を省略するか、'psd' を指定すると、パワー スペクトル密度が返されます。'power' を指定すると、ウィンドウの等価ノイズ帯域幅ごとに PSD 推定をスケーリングします。各周波数のパワーの推定を求めるには 'power' オプションを使用します。

トレース モード。'mean''maxhold' または 'minhold' のいずれかで指定します。既定の設定は 'mean' です。

  • 'mean' — 各入力チャネルのウェルチ スペクトル推定を返します。pwelch は、すべてのセグメントのパワー スペクトル推定を平均化することで各周波数ビンにおけるウェルチ スペクトル推定を計算します。

  • 'maxhold' — 各入力チャネルの最大ホールド スペクトルを返します。pwelch は、すべてのセグメントのパワー スペクトル推定間の最大値を保持することで各周波数ビンにおける最大ホールド スペクトルを計算します。

  • 'minhold' — 各入力チャネルの最小ホールド スペクトルを返します。pwelch は、すべてのセグメントのパワー スペクトル推定間の最小値を保持することで各周波数ビンにおける最小ホールド スペクトルを計算します。

(0,1) の範囲のスカラーとして指定する、真の PSD のカバレッジ確率。出力 pxxc は真の PSD の probability × 100% 区間推定の下限および上限を含みます。

出力引数

すべて折りたたむ

PSD 推定。実数値の非負の列ベクトルまたは行列として返されます。pxx の各列が x の対応する列の PSD 推定です。PSD 推定の単位は、単位周波数あたりの時系列データの 2 乗振幅単位となります。たとえば、入力データがボルト単位の場合、PSD 推定の単位は単位周波数あたりのボルトの 2 乗となります。ボルト単位の時系列では、抵抗が 1 Ωという想定でサンプルレートをヘルツで指定した場合、PSD 推定はワット/ヘルツ単位となります。

データ型: single | double

正規化周波数。実数値の列ベクトルとして返されます。pxx が片側 PSD 推定の場合、w は偶数の nfft に対しては区間 [0,π] をカバーし、奇数の nfft に対しては [0,π) をカバーします。pxx が両側 PSD 推定の場合、w は区間 [0,2π) をカバーします。DC を中央に揃えた PSD 推定の場合、nfft が偶数のときは w は区間 (–π,π] をカバーし、nfft が奇数のときは (–π,π) をカバーします。

データ型: double

巡回周波数。実数値の列ベクトルとして返されます。片側 PSD 推定では、fnfft が偶数の場合は区間 [0,fs/2] をカバーし、nfft が奇数の場合は [0,fs/2) をカバーします。両側 PSD 推定では、f は区間 [0,fs) をカバーします。DC を中央に揃えた PSD 推定では、f は、nfft が偶数長であれば区間 (-fs/2, fs/2] サイクル/単位時間をカバーし、nfft が奇数長であれば (-fs/2, fs/2) サイクル/単位時間をカバーします。

データ型: double

信頼限界。実数値要素をもつ行列として返されます。行列の行のサイズは、PSD 推定 pxx の長さと同じです。pxxc の列数は pxx の 2 倍です。奇数番号の列には信頼区間の下限が、偶数番号の列にはその上限がそれぞれ含まれています。したがって、pxxc(m,2*n-1) は、推定 pxx(m,n) に対応する信頼区間の下限で、pxxc(m,2*n) はその上限となります。信頼区間のカバレッジ確率は、probability の入力値によって決まります。

データ型: single | double

詳細

すべて折りたたむ

ウェルチのオーバーラップ セグメント平均スペクトル推定

ピリオドグラムは、広義における定常プロセスの真のパワー スペクトル密度の一致推定器ではありません。ピリオドグラムの分散を軽減するウェルチの手法では、時系列は通常、オーバーラップのある形でセグメントに分割されます。

ウェルチ法では、セグメントごとに修正ピリオドグラムを計算し、これらの推定を平均してパワー スペクトル密度の推定を求めます。その過程が広義定常性で、ウェルチ法が時系列のさまざまなセグメントの PSD 推定を使用するため、修正ピリオドグラムは真の PSD のほぼ無相関の推定を表現し、平均化によりばらつきが縮小されます。

セグメントは通常ハミング ウィンドウのようなウィンドウ関数で乗算されるため、ウェルチ法は修正ピリオドグラムを平均していることになります。通常セグメントはオーバーラップするため、1 つのセグメント内でウィンドウによってテーパーをかけられたセグメントの始まりと終わりのデータ値は、隣接するセグメントの両端から離れたところにあります。これはウィンドウ処理による情報の損失を防ぎます。

R2006a より前に導入