pmtm

マルチテーパーパワースペクトル密度推定

ページ内をすべて折りたたむ

構文

pxx = pmtm(x)

pxx = pmtm(x,'Tapers',tapertype)

pxx = pmtm(x,nw)

pxx = pmtm(x,m,'Tapers','sine')

pxx = pmtm(___,nfft)

[pxx,w] = pmtm(___)

[pxx,f] = pmtm(___,fs)

[pxx,w] = pmtm(x,nw,w)

[pxx,w] = pmtm(x,m,'Tapers','sine',w)

[pxx,f] = pmtm(___,f,fs)

[___] = pmtm(___,freqrange)

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

[___] = pmtm(___,'DropLastTaper',dropflag)

[___] = pmtm(___,method)

[___] = pmtm(x,e,v,___)

[___] = pmtm(x,dpss_params,___)

pmtm(___)

説明

例

pxx= pmtm(x)は、離散扁長回転楕円体 (スレピアン) 列をテーパーとして使用し、入力信号xのトムソンのマルチテーパーパワースペクトル密度 (PSD) 推定pxxを返します。

例

pxx= pmtm(x,'Tapers',tapertype)は,マルチテーパーPSD推定を計算するときに使用するテーパーのタイプを指定します。'Tapers'とtapertypeの名前と値のペアは、関数呼び出し内のxの後の任意の場所で指定できます。

例

pxx= pmtm(x,nw)は、時間と半帯域幅の積nwを使用し、スレピアンテーパーを使って PSD 推定を計算するときの周波数分解能を制御します。

例

pxx= pmtm(x,m,'Tapers','sine')は、正弦テーパーを使用して PSD 推定値を計算するときに適用するテーパー数または平均化重みを指定します。

例

pxx= pmtm(___,nfft)は、nfft離散フーリエ変換 (DFT) ポイントを前述の任意の構文と組み合わせて使用します。nfftが信号長よりも大きい場合、xには長さnfftまでゼロが追加されます。nfftが信号長よりも小さい場合、信号はnfftを法としてラップされます。

[pxx,w] = pmtm(___)は、pxxが計算される正規化周波数のベクトルを返します。

例

[pxx,f] = pmtm(___,fs)は、周波数ベクトルfを単位時間あたりのサイクル数で返します。fsは、関数呼び出し内のx、nw(正弦テーパーの場合はm)、およびnfftの後で指定しなければなりません。サンプルレートを入力した場合でも、前の引数の既定値を使用するには、これらの引数を空[]として指定します。

[pxx,w] = pmtm(x,nw,w)は、wで指定された正規化周波数でスレピアン列を使用して計算されたマルチテーパー PSD 推定を返します。ベクトルwには少なくとも 2 つの要素が含まれていなければなりません。

例

[pxx,w] = pmtm(x,m,'Tapers','sine',w)は、wで指定された正規化周波数で正弦テーパーを使用して計算されたマルチテーパー PSD 推定を返します。ベクトルwには少なくとも 2 つの要素が含まれていなければなりません。

例

[pxx,f] = pmtm(___,f,fs)は、fで指定された周波数でマルチテーパー PSD 推定を計算します。ベクトルfには、サンプルレートfsと同じ単位で少なくとも 2 つの要素が含まれていなければなりません。

例

[___] = pmtm(___,freqrange)では、freqrangeで指定される周波数範囲全体のマルチテーパー PSD 推定が返されます。

例

[___,pxxc] = pmtm(___,'ConfidenceLevel',probability)は、PSD 推定のprobability× 100% 信頼区間をpxxcに返します。

例

[___] = pmtm(___,'DropLastTaper',dropflag)は、マルチテーパー PSD 推定を計算するときにpmtmが最後のスレピアンテーパーを無視するかどうかを指定します。

例

[___] = pmtm(___,method)は、methodで指定されたメソッドを使用して、個々のテーパー PSD 推定を結合します。この構文は、スレピアンテーパーにのみ適用されます。

例

[___] = pmtm(x,e,v,___)は、eのスレピアンテーパーとvの固有値を使用して PSD を計算します。dpssを使用し、eおよびvを取得します。

例

[___] = pmtm(x,dpss_params,___)は、cell 配列dpss_paramsを使用して、入力引数をdpssに渡します。この構文は、スレピアンテーパーにのみ適用されます。

例

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

例

すべて折りたたむ

既定の入力によるマルチテーパー推定

ライブスクリプトを開く

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

N(0,1) 加法性ホワイトノイズを伴う角周波数 $π / 4$ ラジアン/サンプルの正弦波を作成します。信号長は320サンプルです。既定の時間と半帯域幅の積 (4) と DFT 長を使用してマルチテーパー PSD 推定を求めます。既定の DFT の点の数は 512 です。信号が実数値なので、PSD 推定は片側で、PSD 推定には 512/2+1 点が含まれます。

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

マルチテーパーPSD推定をプロットします。

pmtm(x)

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains an object of type line.

正弦テーパーとスレピアンテーパー

ライブスクリプトを開く

N(0,1) 加法性ホワイトガウスノイズに組み込まれた 2 チャネル信号の 2048 サンプルを生成します。

最初のチャネルは、π/3 およびπ/5 ラジアン/サンプルの正規化周波数をもつ 2 つの正弦波で構成されます。最初の正弦波は 2 番目の振幅の 2 倍です。
2 番目のチャネルの正規化周波数は、π/4 ラジアン/サンプルです。

マルチテーパー法を使用して、0.1πラジアン/サンプルから 0.4πラジアン/サンプルまでの 1024 サンプル区間で信号の PSD を推定します。均等に重み付けされた 13 の正弦テーパーを使用します。

n = (0:2047)'; x = [sin(pi./[3 5].*n)*[2 1]' sin(pi/4*n)] + randn(length(n),2); w = linspace(0.1,0.4,1024); ntp = 13; pmtm(x,ntp,'Tapers','sine',w*pi)

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains 2 objects of type line.

計算を繰り返しますが、今度は 13 のテーパーを線形の降順で重み付けします。'Tapers','sine'の名前と値のペアは、関数呼び出し内のxの後の任意の場所に配置できます。

pmtm(x,(ntp:-1:1)/sum(1:ntp),w*pi,'Tapers','sine')

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains 2 objects of type line.

計算を繰り返しますが、今度は 13 のスレピアンテーパーを使用し、時間と半帯域幅の積を 7.5 に指定します。

nw = 7.5; pmtm(x,{nw,ntp},w*pi)

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains 2 objects of type line.

計算を繰り返しますが、今度はサンプルレートを 2 kHz に指定します。

fs = 2e3; pmtm(x,{nw,ntp},w*(fs/2),fs)

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains 2 objects of type line.

時間と半帯域幅の積の指定

ライブスクリプトを開く

指定した時間と半帯域幅の積でマルチテーパー PSD 推定を求めます。

N(0,1) 加法性ホワイトノイズを伴う角周波数 $π / 4$ ラジアン/サンプルの正弦波を作成します。信号長は320サンプルです。時間と半帯域幅の積を 2.5 としてマルチテーパー PSD 推定を求めます。分解能帯域幅は $[- 2.5 π / 320, 2.5 π / 320]$ ラジアン/サンプルです。既定の DFT の点の数は 512 です。信号が実数値なので、PSD 推定は片側で、PSD 推定には 512/2+1 点が含まれます。

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

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains an object of type line.

信号長と等しい DFT 長

ライブスクリプトを開く

N(0,1) 加法性ホワイトノイズを伴う角周波数 $π / 4$ ラジアン/サンプルの離散時間正弦波で構成される入力信号のマルチテーパー PSD 推定を求めます。信号長と同じ DFT 長を使用します。

N(0,1) 加法性ホワイトノイズを伴う角周波数 $π / 4$ ラジアン/サンプルの正弦波を作成します。信号長は320サンプルです。時間と半帯域幅の積を 3、DFT 長を信号長と同じにしてマルチテーパー PSD 推定を求めます。信号は実数値のため、既定では長さが 320/2+1 に等しい片側 PSD 推定が返されます。

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

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains an object of type line.

サンプルレートによるマルチテーパー推定

ライブスクリプトを開く

1 kHz でサンプリングされた信号のマルチテーパー PSD 推定を求めます。信号は、N(0,1) 加法性ホワイトノイズを伴う 100 Hz の正弦波です。信号の持続時間は 2 秒です。時間と半帯域幅の積 3 および信号長と同じ DFT 長を使用します。

fs = 1000; t = 0:1/fs:2-1/fs; x = cos(2*pi*100*t)+randn(size(t)); [pxx,f] = pmtm(x,3,length(x),fs);

マルチテーパーPSD推定をプロットします。

pmtm(x,3,length(x),fs)

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains an object of type line.

重み 1 の単一テーパー推定平均

ライブスクリプトを開く

個々にテーパーをかけた直接スペクトル推定が平均の等しい重みを与えられている、マルチテーパー PSD 推定を求めます。

1 kHz でサンプリングされた信号のマルチテーパー PSD 推定を求めます。信号は、N(0,1) 加法性ホワイトノイズを伴う 100 Hz の正弦波です。信号の持続時間は 2 秒です。時間と半帯域幅の積 3 および信号長と同じ DFT 長を使用します。'unity'オプションを使用して、テーパーをかけたそれぞれの直接スペクトル推定に等しい平均重みを与えます。

fs = 1000; t = 0:1/fs:2-1/fs; x = cos(2*pi*100*t)+randn(size(t)); [pxx,f] = pmtm(x,3,length(x),fs,'unity');

マルチテーパーPSD推定をプロットします。

pmtm(x,3,length(x),fs,'unity')

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains an object of type line.

DPSS 列およびその周波数領域の集中度

ライブスクリプトを開く

この例では、DPSS 列の周波数領域の集中度を調べます。例では、スレピアン列を事前計算して、分解能帯域幅におけるエネルギーの集中度が 99% を超えるものを選択することで、入力信号のマルチテーパー PSD 推定を求めます。

信号は、N(0,1) 加法性ホワイトノイズを伴う 100 Hz の正弦波です。信号の持続時間は 2 秒です。

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

時間と半帯域幅の積を 3.5 に設定します。信号長が 2000 サンプル、サンプリング間隔が 0.001 秒の場合、分解能帯域幅は [-1.75,1.75] Hz となります。最初の 10 個のスレピアン列を計算し、指定された分解能帯域幅における周波数の集中度を調べます。

[e,v] = dpss(length(x),3.5,10); lv = length(v); stem(1:lv,v,'filled') ylim([0 1.2]) title('Proportion of Energy in [-w,w] of k-th Slepian Sequence')

Figure contains an axes. The axes with title Proportion of Energy in [-w,w] of k-th Slepian Sequence contains an object of type stem.

エネルギーの集中度が 99% を超えるスレピアン列の数を判別します。選択した DPSS 列を使用してマルチテーパー PSD 推定を求めます。'DropLastTaper'をfalseに設定し、選択されたテーパーをすべて使用します。

holdonplot([1 lv],0.99*[1 1]) holdoff

Figure contains an axes. The axes with title Proportion of Energy in [-w,w] of k-th Slepian Sequence contains 2 objects of type stem, line.

idx = find(v>0.99,1,'last')

idx = 5

[pxx,f] = pmtm(x,e(:,1:idx),v(1:idx),length(x),fs,'DropLastTaper',false);

マルチテーパーPSD推定をプロットします。

pmtm(x,e(:,1:idx),v(1:idx),length(x),fs,'DropLastTaper',false)

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains an object of type line.

DC を中央に揃えたマルチテーパー PSD 推定

ライブスクリプトを開く

N(0,1) 加法性ノイズを伴う 100 Hz 正弦波のマルチテーパー PSD 推定を求めます。データは 1 kHz でサンプリングされます。'centered'オプションを使用して DC を中央に揃えた PSD を求めます。

fs = 1000; t = 0:1/fs:2-1/fs; x = cos(2*pi*100*t)+randn(size(t)); [pxx,f] = pmtm(x,3.5,length(x),fs,'centered');

DC を中央に揃えた PSD 推定をプロットします。

pmtm(x,3.5,length(x),fs,'centered')

Figure contains an axes. The axes with title Power Spectral Density contains an object of type line.

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

ライブスクリプトを開く

次の例は、マルチテーパー PSD 推定を伴う信頼限界の使い方を示します。統計的有意性の必要条件ではありませんが、周囲の PSD 推定の信頼下限が信頼上限を超えるマルチテーパー PSD 推定の周波数は、時系列で大きな振幅を明確に示します。

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

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

95% 信頼限界でマルチテーパー PSD 推定を求めます。信頼区間と共に PSD 推定をプロットし、100 Hz と 150 Hz 付近の周波数関心領域を拡大します。

[pxx,f,pxxc] = pmtm(x,3.5,length(x),fs,'ConfidenceLevel', 0.95);情节(f, 10 * log10 (pxx)) holdon情节(f, 10 * log10 (pxxc),'r-.') xlim([85 175]) xlabel('Hz') ylabel('dB') title('Multitaper PSD Estimate with 95%-Confidence Bounds')

Figure contains an axes. The axes with title Multitaper PSD Estimate with 95%-Confidence Bounds contains 3 objects of type line.

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

マルチテーパー法によるマルチチャネル信号の PSD 推定

ライブスクリプトを開く

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

N = 1024; n = 0:N-1; w = pi./[2;3;4]; x = cos(w*n)' + randn(length(n),3); pmtm(x)

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains 3 objects of type line.

入力引数

すべて折りたたむ

`x`—入力信号
ベクトル|行列

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

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

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

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

`tapertype`—テーパータイプ
`'slepian'`(既定値) |`'sine'`

テーパータイプ。'slepian'または'sine'として指定します。

'slepian'— テーパーとして離散扁長回転楕円体 (スレピアン) 列を使用します。
'sine'—正弦テーパーを使用します。

'Tapers'とtapertypeの名前と値のペアは、関数呼び出し内のxの後の任意の場所で指定できます。

データ型:char|string

`nw`—時間と半帯域幅の積
`4`(既定値) |正のスカラー

時間と半帯域幅の積。正のスカラーとして指定します。pmtmは、PSD 推定で 2 ×nw– 1 のスレピアンテーパーを使用します。nwの一般的な選択肢は、2、5/2、3、7/2、または4です。

マルチテーパースペクトル推定では、マルチテーパー推定[–W,W]の分解能帯域幅を指定します。ここで、小さなk > 1については、W = k/NΔtとなります。つまり、W は DFT の周波数分解能を数倍したものです。時間と半帯域幅の積は、分解能半帯域幅と入力信号のサンプル数 N の積です。フーリエ変換が[–W,W](1 に近い固有値) にはっきり集中しているスレピアンテーパーの数は2NW – 1です。

`m`—正弦テーパー数または平均化重み
`7`(既定値) |整数値スカラー|ベクトル

正弦テーパー数または平均化重み。整数値スカラーまたはベクトルとして指定します。

mがスカラーの場合、PSD 推定を計算するときにデータウィンドウとして使用される正弦テーパー数を示します。正弦テーパーは均一に重み付けされています。
mがベクトルの場合、PSD 推定を計算するときに正弦テーパーを平均化するために使用される重みを示します。mの長さは、使用するテーパー数を示します。mの要素の合計は 1 でなければなりません。

データ型:single|double

`nfft`—DFT 点の数
`max(256,2^nextpow2(length(x)))`(既定値) |整数|`[]`

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

データ型:single|double

`fs`—サンプルレート
正のスカラー

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

`w`—正規化周波数
ベクトル

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

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

データ型:double

`f`—周波数
ベクトル

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

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

データ型:double

`dropflag`—最後の DPSS 列を無視するか残すかを示すフラグ
`true`(既定値) |`false`

最後の DPSS 列を落とすか残すかを指示するフラグ。logical として指定します。既定値はtrueで、pmtmは最後のテーパーを無視します。マルチテーパー推定では、最初の 2NW – 1 DPSS 列は 1 に近い固有値をもちます。2NW – 1 未満のシーケンスを使う場合、すべてのテーパーが 1 に近い固有値をもつ傾向があります。dropflagをfalseとして指定すれば、最後のテーパーを残すことができます。

`method`—個々のテーパー PSD 推定の重み
`'adapt'`(既定値) |`'eigen'`|`'unity'`

'adapt'、'eigen'または'unity'で指定する個々のテーパー PSD 推定の重み。既定値は、トムソンの適応的周波数依存重み'adapt'です。これらの重みの計算の詳細は、[2]の 368 ～ 370 ページで説明されています。'eigen'メソッドは、対応するスレピアンテーパーの固有値 (周波数の集中度) によって、各テーパー PSD 推定を重み付けします。'unity'メソッドは各テーパー PSD 推定を同等に重み付けします。

`e`—DPSS (スレピアン) 列
行列

DPSS (スレピアン) 列。行列として指定します。xの長さが N の場合、eは N 行になります。行列eはdpssの出力です。

`v`—DPSS (スレピアン) 列の固有値
ベクトル

列ベクトルとして指定する DPSS (スレピアン) 列の固有値。DPSS 列の固有値は、分解能帯域幅[–W, W]に集中したシーケンスのエネルギーの配分を示します。固有値の範囲は(0, 1)の区間にあり、一般に最初の2NW – 1個の固有値は 1 に近く、その後 0 に向かって減少します。ベクトルvはdpssの出力です。

`dpss_params`—`dpss`の入力引数
cell 配列

cell 配列として指定するdpssの入力引数。dpssへの最初の入力引数は DPSS 列の長さであり、xの長さから取得されるため、dpss_paramsからは省略されます。

例:pmtm(randn(1000,1),{2.5,3})は、時間と半帯域幅の積が 2.5 の最初の 3 つのスレピアン列を使用して、ランダムシーケンスの PSD を計算します。

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

`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 推定はワット/ヘルツ単位となります。

データ型:single|double

`w`— 正規化周波数
ベクトル

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

データ型:double

`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) サイクル/単位時間をカバーします。

データ型:double|single

`pxxc`— 信頼限界
行列

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

データ型:single|double

詳細

すべて折りたたむ

トムソンのマルチテーパースペクトル推定

ピリオドグラムは、広義における定常過程の真のパワースペクトル密度 (PSD) の一致推定器ではありません。ピリオドグラムの変動性を軽減して、PSD の一貫性のある推定を求めるために、マルチテーパー法は相互に直行するウィンドウ、つまり"テーパー"群を使用して求められた修正ピリオドグラムを平均します。互いに直交であることに加え、テーパーは最適な時間および周波数の集中度特性をもちます。テーパーの相互直交性、時間および周波数の集中度のどちらもが、マルチテーパー手法を成功させるうえで欠かせません。トムソンのマルチテーパー法で使用されるスレピアン列の概要については、離散扁長回転楕円体 (スレピアン) 列を参照してください。

マルチテーパー法では、それぞれ異なるスレピアン列をウィンドウとして使って求められた K 個の修正ピリオドグラムを使用します。

$S_{k} (f) = Δ t {| \sum_{n = 0}^{N - 1} g_{k} (n) x (n) e^{- j 2 π f n Δ t} |}^{2}$

これは k 番目のスレピアン列 g_k(n) を使用して求められた修正ピリオドグラムを示します。最も単純な形では、マルチテーパー法は K 個の修正ピリオドグラムをそのまま平均してマルチテーパー PSD 推定を求めます。

$S^{(MT)} (f) = \frac{1}{K} \sum_{k = 0}^{K - 1} S_{k} (f) .$

トムソンのマルチテーパー法 ([4]で説明) は、ウェルチのオーバーラップセグメント平均化の方法に似ています。両方とも、ほぼ無相関な PSD 推定の平均をとります。しかし、2 つのアプローチはこれらの無相関の PSD 推定の求め方に違いがあります。マルチテーパー法は各修正ピリオドグラムにおいて信号全体を使用します。スレピアンテーパーの直交性により、それぞれの修正ピリオドグラムの相関が失われます。ウェルチの方法では、各修正ピリオドグラムにおける信号のセグメントを使用します。セグメント化により、それぞれの修正ピリオドグラムの相関が失われます。

S^(MT)(f)の方程式は、pmtmの'unity'オプションに相当します。しかし、離散扁長回転楕円体 (スレピアン) 列で説明されているように、スレピアン列は対象周波数帯域において等価なエネルギー集中度をもちません。スレピアン列の次数が高くなると、固有値によって与えられる集中度をもつ帯域[–W,W]におけるシーケンスのエネルギーの集中度は低くなります。したがって、平均をとる前に固有値を使用して K 個の修正ピリオドグラムに重み付けすることが有益となります。これはpmtmにおける'eigen'オプションに相当します。

シーケンスの固有値を使用して修正ピリオドグラムの加重平均を求めた場合、スレピアン列における周波数の集中度特性が現れます。しかし、ランダム過程のパワースペクトル密度とスレピアン列の周波数の集中度との相互作用は明らかにはなりません。特に、ランダム過程がパワーをほとんどもたない周波数領域については、修正ピリオドグラムで高次のスレピアン列が使用され、信頼性が低い状態で推定されます。これは周波数依存の適応過程について唱えられており、スレピアン列の周波数の集中度のみならず、時系列におけるパワー分布についての説明にもなっています。この適応重み付けはpmtmにおける'adapt'オプションに相当し、マルチテーパー推定演算の既定値です。

離散扁長回転楕円体 (スレピアン) 列

スレピアン列の偏差は、離散時間または連続的な周波数集中度の問題から生じます。この問題では、インデックスが0, 1, …, N – 1に制限されるすべてのℓ²列について、|W| < 1/2Δtの場合に周波数帯域幅[–W, W]においてエネルギーの集中度が最大になるような列を求めます。

これは、N 行 N 列の半正定値の自己共役作用素について、その固有値と対応する固有ベクトルを求めることを意味しています。したがって、固有値は実数かつ非負であり、各固有値に対応する固有ベクトルは互いに直交になります。特にこの問題では、固有値は 1 以下に制限され、周波数範囲[–W, W]におけるこのシーケンスのエネルギー集中度の測度となります。

固有値問題は次により与えられます。

$\sum_{m = 0}^{N - 1} \frac{\sin (2 π W (n - m))}{π (n - m)} g_{k} (m) = λ_{k} (N, W) g_{k} (n), n, k = 0, 1, 2, \dots, N - 1.$

0 次の DPSS 列g₀は最大の固有値に対応する固有ベクトルです。1 次の DPSS 列g₁は 2 番目に大きい固有値に対応する固有ベクトルであり、0 次のシーケンスと直交します。2 次の DPSS 列g₂は 3 番目に大きい固有値に対応する固有ベクトルであり、2 つの低次の DPSS 列と直交します。演算子は N 行 N 列であるため、N 個の固有ベクトルがあります。しかし、与えられたシーケンスの長さ N と指定の帯域幅[–W, W]について、1 に非常に近い固有値をもつ約2NW – 1の DPSS 列が存在することが確認されます。nwを使用してNWを指定します。

正弦テーパー

[3]で推奨されているスレピアン列の代わりになる正弦テーパーは、次で定義されます。

$g_{k} (n) = \sqrt{\frac{2}{N + 1}} \sin \frac{π k n}{N + 1}, n, k = 1, 2, \dots, N .$

スレピアン列とは異なり、正弦テーパーは直接計算できます。固有値の方程式を設定して解く必要はありません。これにより、正弦テーパーの計算がはるかに高速になります。正弦テーパーのスペクトルの集中度は、スレピアン列のスペクトルの集中度に近いですが、スペクトルの帯域幅を指定するための追加パラメーターは必要ありません。正弦テーパーを使用して計算された PSD 推定の帯域幅は、mを使用してテーパー数を変更することでローカルで調整できます。

スレピアンテーパーと正弦テーパーの比較

時間と半帯域幅の積 3 に対応する最初の 5 つのスレピアンテーパーを生成します。テーパーの長さを 1000 に指定します。

N = 1000; nw = 3; ns = 2*(nw)-1; tprs = dpss(N,nw,ns); lbs ="Slepian";

最初の 5 つの正弦テーパーを生成します。

n = 1:N; k = 1:ns; tprs(:,:,2) = sqrt(2/(N+1))*sin(pi*n'*k/(N+1)); lbs(2) ="Sine";

2 つのセットのテーパーをプロットします。

forkj = 1:2 subplot(2,1,kj) plot(tprs(:,:,kj)) title(lbs(kj)) legend(append('k = ',string(k+kj-2)),...'Orientation','horizontal','Location','south') legend('boxoff') ylim([-0.09 0.07])end

参照

[1] McCoy, Emma J., Andrew T. Walden, and Donald B. Percival. "Multitaper Spectral Estimation of Power Law Processes." IEEE^®Transactions on Signal Processing 46, no. 3 (March 1998): 655–68.https://doi.org/10.1109/78.661333.

[2] Percival, Donald B., and Andrew T. Walden. Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques. Cambridge; New York, NY, USA: Cambridge University Press, 1993.

[3] Riedel, Kurt S., and Alexander Sidorenko. “Minimum Bias Multiple Taper Spectral Estimation.” IEEE Transactions on Signal Processing 43, no. 1 (January 1995): 188–95.https://doi.org/10.1109/78.365298.

[4] Thomson, David J. "Spectrum estimation and harmonic analysis." Proceedings of the IEEE 70, no. 9 (1982): 1055–96.https://doi.org/10.1109/PROC.1982.12433.

参考

dpss|periodogram|pwelch

トピック

ピリオドグラムにおけるバイアスとばらつき

R2006a より前に導入

pmtm

構文

説明

例

既定の入力によるマルチテーパー推定

正弦テーパーとスレピアン テーパー

時間と半帯域幅の積の指定

信号長と等しい DFT 長

サンプルレートによるマルチテーパー推定

重み 1 の単一テーパー推定平均

DPSS 列およびその周波数領域の集中度

DC を中央に揃えたマルチテーパー PSD 推定

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

マルチテーパー法によるマルチチャネル信号の PSD 推定

入力引数

x—入力信号ベクトル|行列

tapertype—テーパー タイプ'slepian'(既定値) |'sine'

nw—時間と半帯域幅の積4(既定値) |正のスカラー

m—正弦テーパー数または平均化重み7(既定値) |整数値スカラー|ベクトル

nfft—DFT 点の数max(256,2^nextpow2(length(x)))(既定値) |整数|[]

fs—サンプルレート正のスカラー

w—正規化周波数ベクトル

f—周波数ベクトル

dropflag—最後の DPSS 列を無視するか残すかを示すフラグtrue(既定値) |false

method—個々のテーパー PSD 推定の重み'adapt'(既定値) |'eigen'|'unity'

e—DPSS (スレピアン) 列行列

v—DPSS (スレピアン) 列の固有値ベクトル

dpss_params—dpssの入力引数cell 配列

freqrange—PSD 推定の周波数範囲'onesided'|'twosided'|'centered'

probability—PSD 推定の信頼区間0.95(既定値) |(0,1) の範囲のスカラー

出力引数

pxx— PSD 推定ベクトル | 行列

w— 正規化周波数ベクトル

f— 巡回周波数ベクトル

pxxc— 信頼限界行列

詳細

トムソンのマルチテーパー スペクトル推定

離散扁長回転楕円体 (スレピアン) 列

正弦テーパー

スレピアン テーパーと正弦テーパーの比較

参照

参考

トピック

Signal Processing Toolbox ドキュメンテーション

サポート

MATLAB による 信号処理向けディープラーニング

正弦テーパーとスレピアンテーパー

`x`—入力信号
ベクトル|行列

`tapertype`—テーパータイプ
`'slepian'`(既定値) |`'sine'`

`nw`—時間と半帯域幅の積
`4`(既定値) |正のスカラー

`m`—正弦テーパー数または平均化重み
`7`(既定値) |整数値スカラー|ベクトル

`nfft`—DFT 点の数
`max(256,2^nextpow2(length(x)))`(既定値) |整数|`[]`

`fs`—サンプルレート
正のスカラー

`w`—正規化周波数
ベクトル

`f`—周波数
ベクトル

`dropflag`—最後の DPSS 列を無視するか残すかを示すフラグ
`true`(既定値) |`false`

`method`—個々のテーパー PSD 推定の重み
`'adapt'`(既定値) |`'eigen'`|`'unity'`

`e`—DPSS (スレピアン) 列
行列

`v`—DPSS (スレピアン) 列の固有値
ベクトル

`dpss_params`—`dpss`の入力引数
cell 配列

`freqrange`—PSD 推定の周波数範囲
`'onesided'`|`'twosided'`|`'centered'`

`probability`—PSD 推定の信頼区間
`0.95`(既定値) |(0,1) の範囲のスカラー

`pxx`— PSD 推定
ベクトル | 行列

`w`— 正規化周波数
ベクトル

`f`— 巡回周波数
ベクトル

`pxxc`— 信頼限界
行列

トムソンのマルチテーパースペクトル推定

スレピアンテーパーと正弦テーパーの比較

MATLAB による信号処理向けディープラーニング