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,___,trace)

[___,pxxc] = pwelch(___,'ConfidenceLevel',probability)

[___] = pwelch(___,spectrumtype)

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 が整数の場合、noverlap は window より小さい正の整数でなければなりません。window がベクトルの場合、noverlap は window の長さより小さい正の整数でなければなりません。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) をカバーします。fs は pwelch の 5 番目の入力でなければなりません。サンプル レートを入力した場合でも、前のオプション引数の既定値を使用するには、これらの引数を空 [] として指定します。

[pxx,w] = pwelch(x,window,noverlap,w) では、ベクトル w で指定される正規化周波数での両側ウェルチ PSD 推定が返されます。ベクトル w には少なくとも 2 つの要素が含まれていなければなりません。そうでない場合は、関数が nfft として解釈するためです。

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

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

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

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

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

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

例

既定の入力によるウェルチ推定

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

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

rng default

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

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

pxx = pwelch(x);

pwelch(x)

計算を繰り返します。

• 長さ $$nsc=\lfloor N_x/4.5\rfloor$$ のセクションに信号を分割します。このアクションは、可能な限り長いセグメントに分割される信号と等価で、8 にできるだけ近い (ただし 8 を超えない) 数のセグメントが 50% のオーバーラップで取得されます。

• ハミング ウィンドウを使用して、セクションにウィンドウを適用します。

• 隣接するセクション間で 50% のオーバーラップを指定します。

• FFT を計算するには、$$\max (256,2^p)$$ 点を使用します。ここで、$$p=\lceil {\log }_{2}nsc\rceil$$ です。

2 つの方法の結果が同じになることを確認します。

Nx = length(x);
nsc = floor(Nx/4.5);
nov = floor(nsc/2);
nff = max(256,2^nextpow2(nsc));

t = pwelch(x,hamming(nsc),nov,nff);

maxerr = max(abs(abs(t(:))-abs(pxx(:))))

maxerr = 
0

長さが等しく、セクション間のオーバーラップが 50% である 8 つのセクションに信号を分割します。前のステップと同じ FFT 長を指定します。ウェルチ PSD 推定を計算し、前の 2 つの手順と結果が同じであることを確認します。

ns = 8;
ov = 0.5;
lsc = floor(Nx/(ns-(ns-1)*ov));

t = pwelch(x,lsc,floor(ov*lsc),nff);

maxerr = max(abs(abs(t(:))-abs(pxx(:))))

maxerr = 
0

指定セグメント長によるウェルチ推定

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

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

rng default

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

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

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

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

セグメントのオーバーラップ指定によるウェルチ推定

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

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

rng default

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

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

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

plot(10*log10(pxx))

指定 DFT 長によるウェルチ推定

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

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

rng default

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

信号を 100 サンプルの長さのセグメントに分割するウェルチ PSD 推定を求めます。50% の既定オーバーラップを使用します。DFT 長が 640 点となるように指定し、$$\pi /4$$ ラジアン/サンプルの周波数が 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)')

ウェルチ法による Hz 単位周波数信号の PSD 推定

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)')

最大ホールド スペクトルと最小ホールド スペクトル

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')

95% 信頼限界の上限および下限

この例は、ウェルチのオーバーラップ セグメント平均 (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 の近傍外における信頼限界の上限より大幅に上になります。

DC を中央に揃えたパワー スペクトル

$$N(0,1/4)$$ 加法性ホワイト ノイズを伴う 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 からの偏差は加法性ノイズの影響によるものです。

ウェルチ法によるマルチチャネル信号の PSD 推定

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

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

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

pwelch(x)

入力引数

x — 入力信号

ベクトル | 行列

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

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

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

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

window — ウィンドウ

整数 | ベクトル | []

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

メモ

ハミング ウィンドウの既定のサイドローブ減衰量は 42.5 dB です。この減衰は、この値 (ピーク スペクトル コンテンツを基準とする値) を下回るスペクトル コンテンツをマスクすることがあります。選択するウィンドウを変えることにより、分解能 (たとえば箱型ウィンドウを使用) とサイドローブ減衰量 (たとえばハン ウィンドウを使用) のバランスをとることができます。詳細については、ウィンドウ デザイナーを参照してください。

データ型: single | double

noverlap — オーバーラップするサンプルの数

正の整数 | []

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

nfft — DFT 点の数

max(256,2^nextpow2(length(window))) (既定値) | 整数 | []

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

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

データ型: single | double

fs — サンプル レート

正のスカラー

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

w — 正規化周波数

ベクトル

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

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

データ型: double | single

f — 周波数

ベクトル

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

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

データ型: double | single

freqrange — PSD 推定の周波数範囲

'onesided' | 'twosided' | 'centered'

'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) サイクル/単位時間になります。

  この関数は 0 とナイキスト周波数以外のすべての周波数でパワーを 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) サイクル/単位時間になります。

入力として周波数のベクトル w または f と freqrange を指定した場合、pwelch は freqrange を無視します。

データ型: char | string

spectrumtype — パワー スペクトルのスケーリング

'psd' (既定値) | 'power'

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

trace — トレース モード

'mean' (既定値) | 'maxhold' | 'minhold'

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

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

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

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

probability — PSD 推定の信頼区間

0.95 (既定値) | (0,1) の範囲のスカラー

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

出力引数

pxx — PSD 推定

ベクトル | 行列

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

w — 正規化周波数

ベクトル

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

f — 巡回周波数

ベクトル

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

pxxc — 信頼限界

行列

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

データ型: single | double

詳細

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

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

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

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

参照

[1] Hayes, Monson H. Statistical Digital Signal Processing and Modeling. New York: John Wiley & Sons, 1996.

[2] Stoica, Petre, and Randolph Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.

拡張機能

tall 配列

メモリの許容量を超えるような多数の行を含む配列を計算します。

pwelch 関数は tall 配列をサポートしますが、次の使用上の注意と制限事項が伴います。

• 入力 x は tall 行ベクトルにしないでください。

• 引数 window は常に指定しなければなりません。

詳細については、tall 配列を参照してください。

C/C++ コード生成
MATLAB® Coder™ を使用して C および C++ コードを生成します。

GPU コード生成
GPU Coder™ を使用して NVIDIA® GPU のための CUDA® コードを生成します。

スレッドベースの環境

MATLAB® の backgroundPool を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool を使用してコードを高速化します。

使用上の注意および制限:

• 出力引数のない構文はサポートされません。

詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

この関数は、GPU 配列を完全にサポートします。詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

R2024b: pwelch は GPU コード生成をサポートする

関数 pwelch は、グラフィックス処理装置 (GPU) 用のコード生成をサポートします。CUDA^® コードを生成するには、MATLAB^® Coder™ および GPU Coder™ が必要です。

R2024a: 単精度可変サイズ ウィンドウ入力のためのコード生成サポート

関数 pwelch は、コード生成用の単精度可変サイズ ウィンドウ入力をサポートします。

R2023b: 単精度データの使用

関数 pwelch は単精度入力をサポートします。

R2023a: [プロットの作成] ライブ エディター タスクを使用した関数出力の可視化

[プロットの作成] ライブ エディター タスクを使用して、pwelch の出力を対話的に可視化することができるようになりました。さまざまなグラフ タイプを選択し、オプションのパラメーターを設定できます。タスクは、ライブ スクリプトの一部となるコードも自動的に生成します。

参考

アプリ

• 信号アナライザー

関数

• periodogram | pmtm | pspectrum

トピック

• ピリオドグラムにおけるバイアスとばらつき
• スペクトル解析