pwelch
ウェルチのパワー スペクトル密度推定
構文
説明
[ は、信号 pxx,f] = pwelch(x,win,nOverlap,freqSpec)x のウェルチ PSD 推定と周波数 f (ラジアン/サンプル) を返します。このとき、pwelch 関数は次のように動作します。
winを使用して信号をセグメントに分割し、各セグメントにウィンドウを適用する。nOverlapで指定したオーバーラップ長を使用して、隣接するセグメント間でサンプルをオーバーラップさせる。ウィンドウが適用された各セグメントに対し、
freqSpecで指定した離散フーリエ変換 (DFT) 点の数に基づいて、または指定した周波数において DFT を計算する。
これらの入力引数のいずれかに既定値を使用するには、それらを空 [] として指定します。
[___] = pwelch(___, は、上記の任意の構文に対して、周波数範囲、スペクトル タイプ、およびトレース モードも追加で指定します。これらの入力引数は、任意の組み合わせで指定できます。freqRange,spectrumType,trace)
出力引数を設定せずに pwelch(___) を使用すると、現在の Figure ウィンドウにウェルチ PSD 推定またはパワー スペクトルがプロットされます。
例
加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの離散時間正弦波で構成される入力信号のウェルチ PSD 推定を求めます。
加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの正弦波を作成します。再現可能な結果が必要な場合は、乱数発生器をリセットします。信号長は サンプルです。
rng("default")
n = 0:319;
x = cos(pi/4*n) + randn(size(n));既定のハミング ウィンドウと DFT 長を使用してウェルチ PSD 推定を求めます。既定のセグメント長は 71 サンプル、DFT 長は 256 点で、周波数分解能は ラジアン/サンプルとなっています。信号は実数値のため、ピリオドグラムは片側で、256/2+1 の点が含まれます。ウェルチ PSD 推定をプロットします。
pxx = pwelch(x); pwelch(x)
計算を繰り返します。
長さ のセクションに信号を分割します。このアクションは、可能な限り長いセグメントに分割される信号と等価で、8 にできるだけ近い (ただし 8 を超えない) 数のセグメントが 50% のオーバーラップで取得されます。
ハミング ウィンドウを使用して、セクションにウィンドウを適用します。
隣接するセクション間で 50% のオーバーラップを指定します。
FFT を計算するには、 点を使用します。ここで、 です。
2 つの方法の結果が同じになることを確認します。
Nx = length(x); nsc = floor(Nx/4.5); nov = floor(nsc/2); nff = max(256,2^nextpow2(nsc)); pxxt = pwelch(x,hamming(nsc),nov,nff); maxerr = max(abs(abs(pxxt(:)) - abs(pxx(:))))
maxerr = 0
長さが等しく、セクション間のオーバーラップが 50% である 8 つのセクションに信号を分割します。前のステップと同じ FFT 長を指定します。ウェルチ PSD 推定を計算し、前の 2 つの手順と結果が同じであることを確認します。
ns = 8; ov = 0.5; lsc = floor(Nx/(ns-(ns-1)*ov)); pxxt8 = pwelch(x,lsc,floor(ov*lsc),nff); maxerr8 = max(abs(abs(pxxt8(:)) - abs(pxx(:))))
maxerr8 = 0
加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの離散時間正弦波で構成される入力信号のウェルチ 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,pow2db(pxx)) xlabel("Normalized Frequency (\times \pi rad/sample)") title("Welch Power Spectral Density Estimate")
加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの離散時間正弦波で構成される入力信号のウェルチ 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(pow2db(pxx)) xlabel("Frequency Samples") title("Welch Power Spectral Density Estimate")
加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの離散時間正弦波で構成される入力信号のウェルチ 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(pow2db(pxx)) xlabel("Frequency Points") title("Welch Power Spectral Density Estimate")
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,pow2db(pxx)) xlabel("Frequency (Hz)") ylabel("PSD (dB/Hz)")
200 kHz で 0.1 秒間サンプリングされた、ノイズを含む 3 つの正弦波とチャープで構成される信号を作成します。正弦波の周波数は 1 kHz、10 kHz および 20 kHz です。正弦波の振幅とノイズ レベルはさまざまです。ノイズのないチャープの周波数は 20 kHz で始まり、サンプリング中に 30 kHz まで線形に増加します。
rng("default")
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 サンプル) をオーバーラップさせます。
L = 200; win = hamming(L); nOverlap = 100; [pxx,f,pxxc] = pwelch(x,win,nOverlap,200,Fs,ConfidenceLevel=0.95);
信頼区間と共に WOSA PSD 推定をプロットし、100 Hz と 150 Hz 付近の周波数関心領域を拡大します。100 および 150 Hz のごく近傍における信頼限界の下限は、100 および 150 Hz の近傍外における信頼限界の上限より大幅に上になります。
plot(f,pow2db(pxx)) hold on plot(f,pow2db(pxxc),"-.",Color=[0.866 0.329 0]) hold off xlim([25 250]) xlabel("Frequency (Hz)") ylabel("PSD (dB/Hz)") title("Welch Estimate with 95%-Confidence Bounds")
加法性ホワイト ノイズを伴う 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");
結果をプロットします。-100 Hz と 100 Hz におけるパワーが振幅 1 の実数値正弦波に対して 1/4 となり、期待したパワーに近くなります。1/4 からの偏差は加法性ノイズの影響によるものです。
plot(f,pow2db(pxx)) xlabel("Frequency (Hz)") ylabel("Magnitude (dB)") grid
加法性ホワイト ガウス ノイズを伴う 3 つの正弦波から構成されるマルチチャネル信号の 1024 サンプルを生成します。正弦波の周波数は、、、および ラジアン/サンプルです。ウェルチ法を使用して信号の PSD を推定し、プロットします。
N = 1024;
n = 0:N-1;
w = pi./[2;3;4];
rng("default")
x = cos(w*n)' + randn(length(n),3);
pwelch(x)
R2026a 以降
指定したターゲット座標軸とパネル コンテナーに、4 つの信号のウェルチ パワー スペクトル密度 (PSD) 推定とウェルチ パワー スペクトルをプロットします。
4 つの振動信号 (サンプル レート 10 kHz、3 秒間) を作成します。
Fs = 10e3;
t = 0:1/Fs:3;
x1 = sinc(Fs/2.5*(t-mean(t)));
x2 = sum(cos(2*pi*600*[1 3 5 7]'.*t),1) + randn(size(t))/1e4;
x3 = exp(1j*pi*sin(4*t)*Fs/10);
x4 = chirp(t,Fs/10,t(end),Fs/2.5,"quadratic");ターゲット座標軸へのウェルチ PSD 推定とパワー スペクトルのプロット
新しい Figure ウィンドウの南西隅と北東隅に 2 つの座標軸を作成します。
fig = figure; ax1 = axes(fig,Position=[0.09 0.1 0.52 0.45]); ax2 = axes(fig,Position=[0.55 0.7 0.42 0.25]);
信号 x1 のウェルチ PSD 推定と信号 x2 のウェルチ パワー スペクトルを、Figure の南西軸と北東軸にそれぞれプロットします。256 サンプルのカイザー ウィンドウ、220 サンプルのオーバーラップ長、および 512 個の DFT 点を使用します。
g = kaiser(256,5);
ol = 220;
nfft = 512;
pwelch(x1,g,ol,nfft,Fs,Parent=ax1)
pwelch(x2,g,ol,nfft,Fs,"power",Parent=ax2)
ターゲット UI 座標軸へのウェルチ PSD 推定のプロット
新しい UI Figure ウィンドウの北西隅に座標軸を作成します。
uif = uifigure(Position=[100 100 720 540]); ax3 = uiaxes(uif,Position=[5 305 300 200]);
信号 x3 のウェルチ PSD 推定を Figure の座標軸にプロットします。周波数は 0 kHz を中心として表示します。
pwelch(x3,g,ol,nfft,Fs,"centered",Parent=ax3) title(ax3,"Welch PSD in UI Axes")
ターゲット パネル コンテナーへのウェルチ PSD 推定のプロット
UI Figure ウィンドウの南東隅にパネル コンテナーを追加します。
p = uipanel(uif,Position=[300 5 400 325], ... Title="Welch PSD in Panel Container", ... BackgroundColor="white");
信号 x4 のウェルチ PSD 推定をパネル コンテナーにプロットします。信頼水準を 90% に設定します。
pwelch(x4,g,ol,nfft,Fs,ConfidenceLevel=0.9,Parent=p)
入力引数
出力引数
詳細
参照
[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.
