ドキュメンテーション

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

periodogram

ピリオドグラム パワー スペクトル密度推定

構文

pxx = periodogram(x)
pxx = periodogram(x,window)
pxx = periodogram(x,window,nfft)
[pxx,w] = periodogram(___)
[pxx,f] = periodogram(___,fs)
[pxx,w] = periodogram(x,window,w)
[pxx,f] = periodogram(x,window,f,fs)
[___] = periodogram(x,window,___,freqrange)
[___] = periodogram(x,window,___,spectrumtype)
[___,pxxc] = periodogram(___,'ConfidenceLevel',probability)
[rpxx,f] = periodogram(___,'reassigned')
[rpxx,f,pxx,fc] = periodogram(___,'reassigned')
periodogram(___)

説明

pxx = periodogram(x) では、箱型ウィンドウを使用して検出された入力信号 x のピリオドグラム パワー スペクトル密度 (PSD) 推定 pxx が返されます。x がベクトルの場合、単一チャネルとして取り扱われます。x が行列の場合、PSD は列ごとに個別に計算され、pxx の対応する列に保存されます。x が実数値の場合、pxx は片側 PSD 推定です。x が複素数値の場合、pxx は両側 PSD 推定です。離散フーリエ変換 (DFT) の点の数 nfft の最大値は 256 または 信号長よりも大きい最小の 2 のべき乗です。

pxx = periodogram(x,window) では、ウィンドウ window を使用して、修正ピリオドグラム PSD 推定が返されます。window は、x と同じ長さのベクトルです。

pxx = periodogram(x,window,nfft) は、離散型フーリエ変換 (DFT) で nfft 個の点を使用します。nfft が信号長よりも大きい場合、x には長さ nfft までゼロが追加されます。nfft が信号長よりも小さい場合、信号は nfft を法としてラップされ、datawrap を使用して合計されます。たとえば、nfft が 4 の入力信号 [1 2 3 4 5 6 7 8] は、sum([1 5; 2 6; 3 7; 4 8],2) のピリオドグラムになります。

[pxx,w] = periodogram(___) では、正規化周波数ベクトル w が返されます。pxx が片側ピリオドグラムの場合、wnfft が偶数であれば区間 [0,π] をカバーし、nfft が奇数であれば [0,π) をカバーします。pxx が両側ピリオドグラムの場合、w は区間 [0,2π) をカバーします。

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

[pxx,w] = periodogram(x,window,w) では、ベクトル w で指定される正規化周波数での両側ピリオドグラム推定が返されます。w には少なくとも 2 つの要素が含まれていなければなりません。

[pxx,f] = periodogram(x,window,f,fs) では、ベクトル f で指定される周波数での両側ピリオドグラム推定が返されます。f には少なくとも 2 つの要素が含まれていなければなりません。f の周波数の単位は単位時間あたりのサイクルです。サンプルレート fs は単位時間あたりのサンプル数です。時間の単位が秒の場合、f の単位はサイクル/秒 (Hz) です。

[___] = periodogram(x,window,___,freqrange) では、freqrange で指定される周波数範囲にわたるピリオドグラムが返されます。freqrange に指定できる有効なオプションは、'onesided''twosided' または 'centered' です。

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

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

[rpxx,f] = periodogram(___,'reassigned') では、各 PSD 推定が、そのエネルギーの中心に最も近い周波数に再代入されます。rpxx には f の各要素に再代入された推定の和が含まれます。

[rpxx,f,pxx,fc] = periodogram(___,'reassigned') では、再代入されない PSD 推定 pxx とエネルギー中心の周波数 fc も返されます。'reassigned' フラグを使用する場合、probability 信頼区間は指定できません。

出力引数なしで periodogram(___) を使用すると、現在の Figure ウィンドウに、ピリオドグラム PSD 推定が単位周波数あたり dB 単位でプロットされます。

すべて折りたたむ

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの離散時間正弦波で構成される入力信号のピリオドグラムを求めます。

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの正弦波を作成します。信号長は 320 サンプルです。既定の箱型ウィンドウと DFT 長を使用してピリオドグラムを求めます。DFT 長は信号長よりも大きい最小の 2 のべき乗または 512 点です。信号は実数値で偶数の長さをもつため、ピリオドグラムは片側で 512/2+1 点を含みます。

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
[pxx,w] = periodogram(x);
plot(w,10*log10(pxx))

出力なしで periodogram を使用してプロットを繰り返します。

periodogram(x)

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの離散時間正弦波で構成される入力信号の修正ピリオドグラムを求めます。

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの正弦波を作成します。信号長は 320 サンプルです。ハミング ウィンドウと既定の DFT 長を使用して修正ピリオドグラムを求めます。DFT 長は信号長よりも大きい最小の 2 のべき乗または 512 点です。信号は実数値で偶数の長さをもつため、ピリオドグラムは片側で 512/2+1 点を含みます。

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

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの離散時間正弦波で構成される入力信号のピリオドグラムを求めます。信号長と同じ DFT 長を使用します。

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルの正弦波を作成します。信号長は 320 サンプルです。既定の箱型ウィンドウと信号長と同じ DFT 長を使用してピリオドグラムを求めます。信号は実数値のため、既定では長さが 320/2+1 に等しい片側ピリオドグラムが返されます。

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
nfft = length(x);
periodogram(x,[],nfft)

1700 年から 1987 年にかけて毎年サンプリングされたウォルフ (相対黒点) 数データのピリオドグラムを求めます。

相対黒点数のデータを読み込みます。既定の箱型ウィンドウと DFT の点の数 (この例では 512) を使用してピリオドグラムを求めます。これらのデータのサンプルレートは 1 サンプル/年です。ピリオドグラムをプロットします。

load sunspot.dat
relNums=sunspot(:,2);

[pxx,f] = periodogram(relNums,[],[],1);

plot(f,10*log10(pxx))
xlabel('Cycles/Year')
ylabel('dB / (Cycles/Year)')
title('Periodogram of Relative Sunspot Number Data')

上の図から、およそ 0.1 サイクル/年にピリオドグラムのピークがあり、つまりおよそ 10 年が周期であることがわかります。

加法性ホワイト ノイズを伴う角周波数 ラジアン/サンプルと ラジアン/サンプルの 2 つの離散時間正弦波で構成される入力信号のピリオドグラムを求めます。 および ラジアン/サンプルにおける両側ピリオドグラム推定を求めます。結果を片側ピリオドグラムと比較します。

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

[pxx,w] = periodogram(x,[],[pi/4 pi/2]);
pxx
pxx = 1×2

   14.0589    2.8872

[pxx1,w1] = periodogram(x);
plot(w1/pi,pxx1,w/pi,2*pxx,'o')
legend('pxx1','2 * pxx')
xlabel('\omega / \pi')

得られたピリオドグラムの値は片側ピリオドグラムの値の 1/2 です。特定の周波数のセットでピリオドグラムを評価する場合、出力は両側推定になります。

N(0,1) 加法性ホワイト ノイズを伴う周波数 100 Hz および 200 Hz の 2 つの正弦波で構成される信号を作成します。サンプリング周波数は 1 kHz です。100 Hz と 200 Hz における両側ピリオドグラムを求めます。

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

freq = [100 200];
pxx = periodogram(x,[],freq,fs)
pxx = 1×2

    0.2647    0.2313

次の例は、ピリオドグラムでの信頼限界の使い方を示します。統計的有意性の必要条件ではありませんが、周囲の PSD 推定の信頼下限が信頼上限を超えるピリオドグラムの周波数は、時系列において有意な振幅を明確に示します。

N(0,1) 加法性ホワイト ノイズを伴う 100 Hz と 150 Hz の正弦波を重ね合わせて構成される信号を作成します。2 つの正弦波の振幅は 1 で、サンプリング周波数は 1 kHz です。

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

95% 信頼限界でピリオドグラム PSD 推定を求めます。信頼区間と共にピリオドグラムをプロットし、100 Hz と 150 Hz 付近の周波数関心領域を拡大します。

[pxx,f,pxxc] = periodogram(x,rectwin(length(x)),length(x),fs,...
    'ConfidenceLevel',0.95);

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

xlim([85 175])
xlabel('Hz')
ylabel('dB/Hz')
title('Periodogram with 95%-Confidence Bounds')

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

'power' オプションを使用して特定の周波数で正弦波のパワーを推定します。

1 kHz でサンプリングされた 1 秒間の 100 Hz 正弦波を作成します。正弦波の振幅は 1.8 で、これは 1.8²/2 = 1.62 のパワーと等しくなります。'power' オプションを使用してパワーを推定します。

fs = 1000;
t = 0:1/fs:1-1/fs;
x = 1.8*cos(2*pi*100*t);
[pxx,f] = periodogram(x,hamming(length(x)),length(x),fs,'power');
[pwrest,idx] = max(pxx);
fprintf('The maximum power occurs at %3.1f Hz\n',f(idx))
The maximum power occurs at 100.0 Hz
fprintf('The power estimate is %2.2f\n',pwrest)
The power estimate is 1.62

加法性ノイズを伴う 100 Hz の正弦波のピリオドグラムを求めます。データは 1 kHz でサンプリングされています。'centered' オプションを使用して DC を中央に揃えたピリオドグラムを求め、結果をプロットします。

fs = 1000;
t = 0:0.001:1-0.001;
x = cos(2*pi*100*t)+randn(size(t));
periodogram(x,[],length(x),fs,'centered')

ホワイト ガウス ノイズに含まれる 200 Hz の正弦波で構成される信号を生成します。信号は 1 kHz で 1 秒間サンプリングされます。ノイズの分散は 0.01² です。再現可能な結果が必要な場合は、乱数発生器をリセットします。

rng('default')

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
N = length(t);
x = sin(2*pi*t*200)+0.01*randn(size(t));

FFT を使用して、信号長で正規化された信号のパワー スペクトルを計算します。この正弦波はビン内にあるため、すべてのパワーは単一の周波数サンプルに集中しています。片側スペクトルをプロットします。ピークの近傍にズームインします。

q = fft(x,N);
ff = 0:Fs/N:Fs-Fs/N;

ffts = q*q'/N^2
ffts = 0.4997
ff = ff(1:floor(N/2)+1);
q = q(1:floor(N/2)+1);

stem(ff,abs(q)/N,'*')
axis([190 210 0 0.55])

periodogram を使用して、信号のパワー スペクトルを計算します。ハン ウィンドウと 1024 の FFT 長を指定します。200 Hz における推定パワーと実際の値の間の割合差を求めます。

wind = hann(N);

[pun,fr] = periodogram(x,wind,1024,Fs,'power');

hold on
stem(fr,pun)

periodogErr = abs(max(pun)-ffts)/ffts*100
periodogErr = 4.7349

パワー スペクトルを再度計算しますが、今度は再割り当てを使用します。新しい推定をプロットし、その最大値を FFT 値と比較します。

[pre,ft,pxx,fx] = periodogram(x,wind,1024,Fs,'power','reassigned');

stem(fx,pre)
hold off
legend('Original','Periodogram','Reassigned')

reassignErr = abs(max(pre)-ffts)/ffts*100
reassignErr = 0.0779

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

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

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

periodogram(x)

入力引数

すべて折りたたむ

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

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

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

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

ウィンドウ。入力信号と同じ長さの行ベクトルまたは列ベクトルで指定します。window を空として指定した場合、periodogram は箱型ウィンドウを使用します。'reassigned' フラグと空の window を指定した場合、関数は β = 38 のカイザー ウィンドウを使用します。

データ型: single | double

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

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

データ型: char

パワー スペクトルのスケーリング。'psd' または 'power' で指定します。パワー スペクトル密度を返すには、spectrumtype を省略するか、または 'psd' を指定します。各周波数のパワーの推定を求めるには、代わりに 'power' を使用します。'power' を指定すると、'reassigned' フラグを使用している場合を除き、ウィンドウ相当ノイズ帯域幅ごとの PSD 推定をスケーリングします。

データ型: char

(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 推定。実数値の非負の列ベクトルまたは行列として返されます。rpxx の各列は、x の対応する列の再代入された PSD 推定です。

エネルギー中心の周波数。ベクトルまたは行列として指定します。

詳細

すべて折りたたむ

ピリオドグラム

ピリオドグラムは広義定常性ランダム過程のパワー スペクトル密度 (PSD) のノンパラメトリック推定です。ピリオドグラムは自己相関列のバイアスのある推定のフーリエ変換です。fs サンプル/単位時間でサンプリングされた信号 xn の場合、ピリオドグラムは次のように定義されます。

P^(f)=ΔtN|n=0N1xnei2πfn|2,1/2Δt<f1/2Δt,

ここで、Δt はサンプリング間隔です。片側ピリオドグラムの場合、0 とナイキスト周波数 1/2Δt 以外のすべての周波数における値は合計パワーを保全するため 2 倍されます。

周波数の単位がラジアン/サンプルの場合、ピリオドグラムは次のように定義されます。

P^(ω)=12πN|n=0N1xneiωn|2,π<ωπ.

上記の方程式の周波数範囲は freqrange 引数の値によって異なります。入力引数freqrange の説明を参照してください。

1 周期 (巡回周波数の場合は 1/Δt、正規化周波数の場合は 2π) に対する真の PSD の積分 P(f) は、広義定常性ランダム過程の分散に等しくなります。

σ2=1/2Δt1/2ΔtP(f)df.

正規化周波数の場合は、積分範囲を適切に置き換えます。

修正ピリオドグラム

修正ピリオドグラムは、入力時系列をウィンドウ関数で乗算します。適切なウィンドウ関数は非負であり、開始点と終了点でゼロに減衰します。時系列をウィンドウ関数で乗算すると、データが徐々に断続的に "漸減し"、ピリオドグラムでの漏れを軽減します。例は、ピリオドグラムにおけるバイアスとばらつきを参照してください。

hn がウィンドウ関数の場合、修正ピリオドグラムは次のように定義されます。

P^(f)=ΔtN|n=0N1hnxnei2πfn|2,1/2Δt<f1/2Δt,

ここで、Δt はサンプリング間隔です。

周波数の単位がラジアン/サンプルの場合、修正ピリオドグラムは次のように定義されます。

P^(ω)=12πN|n=0N1hnxneiωn|2,π<ωπ.

上記の方程式の周波数範囲は freqrange 引数の値によって異なります。入力引数freqrange の説明を参照してください。

再割り当てされたピリオドグラム

再割り当て手法では、スペクトル推定の局所化が鮮明になり、読み取りと解釈の容易なピリオドグラムが作成されます。この手法では、各 PSD 推定はビンの幾何学的中心から離れ、そのビンのエネルギー中心に再代入されます。これにより、チャープとインパルスの厳密な局所化が行われます。

参照

[1] Fulop, Sean A., and Kelly Fitz. “Algorithms for computing the time-corrected instantaneous frequency (reassigned) spectrogram, with applications.” Journal of the Acoustical Society of America. Vol. 119, January 2006, pp. 360–371.

[2] Auger, François, and Patrick Flandrin. “Improving the Readability of Time-Frequency and Time-Scale Representations by the Reassignment Method.” IEEE® Transactions on Signal Processing. Vol. 43, May 1995, pp. 1068–1089.

R2006a より前に導入