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) 加法性ホワイト ノイズを伴う角周波数 $$\pi /4$$ ラジアン/サンプルの離散時間正弦波で構成される入力信号のマルチテーパー PSD 推定を求めます。

N(0,1) 加法性ホワイト ノイズを伴う角周波数 $$\pi /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)

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

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)

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

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

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

nw = 7.5;

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

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

fs = 2e3;

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

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

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

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

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

信号長と等しい DFT 長

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

N(0,1) 加法性ホワイト ノイズを伴う角周波数 $$\pi /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))

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

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)

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

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

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

hold on
plot([1 lv],0.99*[1 1])
hold off

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)

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

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

plot(f,10*log10(pxx))
hold on
plot(f,10*log10(pxxc),'r-.')
xlim([85 175])
xlabel('Hz')
ylabel('dB')
title('Multitaper PSD Estimate with 95%-Confidence Bounds')

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

マルチテーパー法によるマルチチャネル信号の 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);

pmtm(x)

入力引数

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 | single

f — 周波数

ベクトル

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

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

データ型: double | single

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

  この関数は 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 を指定した場合、pmtm は freqrange を無視します。

データ型: char | string

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

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

$$S_k(f)=\Delta t{\left|{\displaystyle \sum _{n=0}^{N-1}{g}_k\left(n\right)}x\left(n\right)e^{-j2\pi fn\Delta t}\right|}^2$$

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

$$S^{(\text{MT})}(f)=\frac{1}{K}{\displaystyle \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] におけるこのシーケンスのエネルギー集中度の測度となります。

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

$${\displaystyle \sum _{m=0}^{N-1}\frac{\sin \left(2\pi W\left(n-m\right)\right)}{\pi \left(n-m\right)}}{g}_k\left(m\right)={\lambda }_k\left(N,W\right)g_k\left(n\right),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{\pi kn}{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 つのセットのテーパーをプロットします。

for kj = 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.

拡張機能

C/C++ コード生成

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

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

• 入力引数 nfft および nw は、コード生成時にスカラーでなければなりません。

• nfft と nw が可変サイズの場合、関数は nfft を周波数ベクトル f、nw をスレピアン テーパー e の行列と見なします。

バージョン履歴

R2006a より前に導入

R2023b: C および C++ コードの生成

関数 pmtm の C および C++ コードを生成します。この機能を使用するには、MATLAB^® Coder™ が必要です。

参考

dpss | periodogram | pwelch

トピック

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