Statistical Signal Processing (Signal Processing Toolbox)

Signal Processing Toolbox

ノンパラメトリック法

つぎの節では、ノンパラメトリック推定法のピリオドグラム、修正ピリオドグラム、Welch 法、マルチテーパ(multitaper)法を、関連した CSD 関数、伝達関数推定、コヒーレンス関数と共に説明します。

ピリオドグラム

プロセスのパワースペクトル推定の一つの方法は、プロセスのあるサンプルの離散時間フーリエ変換を求め、結果の大きさを二乗するものです。この推定を、ピリオドグラム と言います。

長さ L の信号 x_L[n] の PSD のピリオドグラム推定は、つぎのように表わせます。

ここで、つぎのように表わせます。

X_L(f) の実際の計算は、有限個の周波数点 N でのみ実行します。この計算は、通常、FFT を使います。特に、ピリオドグラム法で行う作業の大部分は、N-点の PSD 推定を計算することです。

ここで、つぎのように表わせます。

L よりも大きい2のベキ乗数の中で最小の値になるように、N を選択します。 X_L[f_k] を推定するには、x_L[n] にゼロを付加して、長さを N にします。L > N の場合、X_L[f_k] を計算する前に、x_L[n] を N に分けます。

例として、つぎの 1001 要素の信号 xn を考えます。この中には、2つの正弦波と雑音が含まれています。

randn('state',0);
fs = 1000;           % サンプリング周波数
t = (0:fs)/fs;       % 1秒の間にサンプルを分布
A = [1 2];           % 正弦波の振幅 (行ベクトル)
f = [150;140];       % 正弦波の周波数 (列ベクトル)
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));

注意最後の3行は、正弦波の和を表現するもので、

   xn 
      = sin(2*pi*150*t) + 2*sin(2*pi*140*t) + 0.1*randn(size(t));

と等価です。

PSD のピリオドグラム推定は、つぎのステートメントで計算できます。

Pxx = periodogram(xn,[],'twosided',1024,fs);

出力引数を設定しない場合、つぎのように、推定結果をプロット表示します。

periodogram(xn,[],'twosided',1024,fs);

平均パワーは、つぎの和で積分を近似することで、計算できます。

Pow = (fs/length(Pxx)) * sum(Pxx)
Pow =
    2.5028

片側 PSD 推定から平均パワーを計算することもできます。

Pxxo = periodogram(xn,[],1024,fs);
Pow = (fs/(2*length(Pxxo))) * sum(Pxxo)
Pow =
    2.4979

ピリオドグラムの性能

つぎの節は、leakage, resolution, bias, variance による問題に関するピリオドグラムの性能を議論します。

スペクトルリーケージ：ピリオドグラムで議論した、有限長信号 x_L[n] の PSD 、または、パワースペクトルを考えましょう。無限長の信号 x[n] に、有限長の箱型ウインドウ w_R[n] を乗算した結果として、 x_L[n] を解釈することは、しばしば有益なことになります。

時間領域での乗算は、周波数領域でのコンボリューションになるので、上の表現のフーリエ変換は、つぎのようになります。

ピリオドグラム用に前に使われた表現

は、ピリオドグラムに現れるこのコンボリューションの影響を示しています。

コンボリューションの影響は、正弦波データに対して、一番理解し易いものです。x[n] は、M 個の複素正弦波の和から構成されていると仮定します。

このスペクトルは、つぎのようになります。

有限長のデータ列については、つぎのように表わせます。

有限長の信号のスペクトルで、Dirac デルタ関数が、の型の項で置き換えられていました。これは、周波数 f_k を中心とする箱型の周波数応答に対応します。

箱型ウインドウの周波数応答は、つぎに示す Sinc 信号の型をしています。

プロットには、メインロブといくつかのサイドロブが表示されています。サイドロブの中の最大のものは、メインロブより約 13.5 dB 小さくなっています。これらのロブは、スペクトルリーケージ(spectral leakage) として知られています。無限長の信号は、離散周波数点 f_k に厳密な意味で、集中したパワーをもっている一方、ウインドウを適用された(または、打ち切られた)信号は、離散周波数点 f_k の近傍で、"leaked"パワーを示しています。

短い箱型ウインドウの周波数応答は、より長いウインドウの周波数応答よりも Dirac デルタ関数への近似が良くないので、スペクトルリーケージは、データ長が短い場合に非常に重要になります。100サンプルのつぎのデータ列を考えましょう。

randn('state',0)
fs = 1000;           % サンプル周波数
t = (0:fs/10)/fs;    % 1秒の間にサンプルを分布
A = [1 2];           % 正弦波の振幅
f = [150;140];       % 正弦波の周波数
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
periodogram(xn,[],1024,fs);

スペクトルリーケージの影響が、データの長さにのみに依存することは、特筆すべきことです。ピリオドグラムは、周波数サンプルの有限個で計算されると言うことから起因するものではありません。

解像度 解像度 は、個々のスペクトルを区別できるかと言うことで、考えられ、スペクトル推定の性能解析の基本的な考えになっています。

周波数的に比較的に近い2つの正弦波を分離するには、これらの正弦波のどちらか一方に対して、leaked スペクトルのメインロブの幅よりも、2つの周波数間の違いが大きいことが必要です。メインロブ幅は、パワーがピークのメインロブの半分になる点の幅(すなわち、3 dB 幅)で定義されます。この幅は、f_s / L と近似的に等価です。

言い換えれば、周波数 f₁ と f₂ の正弦波に対して、分解可能な条件は、

を満たすことです。

上の例題で、2つの正弦波は、たった10 Hz しか差がありません。ピリオドグラムを使って、2つを明確に分離するために必要な解像度を得るには、100サンプル以上を記録データは必要とします。

この基準が満たされない、すなわち、67サンプルの場合を考えてみましょう。

randn('state',0)
fs = 1000;               % サンプリング周波数
t = (0:fs/15)./fs;       % 67 サンプル
A = [1 2];               % 正弦波の振幅
f = [150;140];           % 正弦波の周波数
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
periodogram(xn,[],1024,fs);

解像度に関する上の議論は、S/N比(SNR)が、比較的高いので、雑音の影響を考えていません。SNR が低い場合、真のスペクトルは、非常に分離が困難で、雑音による影響が、ピリオドグラムをベースにしたスペクトル推定に現れます。つぎに、この例を示します。

randn('state',0)
fs = 1000;           % サンプリング周波数
t = (0:fs/10)./fs;   % 1秒の間にサンプルを分布
A = [1 2];           % 正弦波の振幅
f = [150;140];       % 正弦波の周波数
xn = A*sin(2*pi*f*t) + 2*randn(size(t));
periodogram(xn,[],1024,fs);

ピリオドグラムのバイアス ピリオドグラムは、PSD のバイアス付き推定子です。その期待値は、つぎのように示されます。

これは、スペクトルリーケージの中の X_L(f) に対する最初の表現に似ています。但し、ここでの表現は、大きさよりも、平均パワーの項を使っています。これは、ピリオドグラムで作成された推定は、真の PSD よりむしろ、リーケージ の影響をもつ PSD に対応します。

つぎのことに注意しましょう。

は、本質的に、三角 Bartlett ウインドウ(2 つの箱型パスルのコンボリューションから生じた三角パスル)を与えます。これは、リーケージの影響を受けたパワースペクトルの最大サイドロブが、メインロブのピークの約 27 dB 低い高さを示すことになります。すなわち、正方でない箱型ウインドウに関して、周波数の約2倍の分離を示します。

ピリオドグラムは、漸近的なバイアスの適用されていないもので、これは、データ長が無限大に近付くに連れ、箱型ウインドウの周波数応答は、Dirac デルタ関数により近付く(Bartlett ウインドウでも真)と言う、初期の観測からわかります。しかし、いくつかの場合、ピリオドグラムは、データが長い場合でも、良くない推定子になります。これは、つぎに説明する、ピリオドグラムの分散によるものです。

ピリオドグラムの分散 ピリオドグラムの分散は、近似的に、つぎのように示されます。

これは、データ長 L を無限大にするに連れ、分散は、ゼロ方向に向かないことを示しています。統計的に、ピリオドグラムは、PSD の整合性の取れた推定子ではありません。けれども、ピリオドグラムは、SNR が高い場合のスペクトル推定には、有効なツールです。特に、データが長い場合です。

修正ピリオドグラム

修正ピリオドグラム は、信号の端をスムーズにするために、FFT を計算する前に、時間領域で、信号にウインドウを適用します。これは、サイドロブの高さ、または、スペクトルリーケージのいずれかを小さくする効果があります。この現象は、箱型ウインドウが使われたときに生じる急激な打ち切りにより、信号の中に取り込まれる擬似的な周波数として、サイドロブを解釈する可能性をもたせます。非箱型ウインドウに対して、打ち切られた信号の端の点は、スムーズに減衰していきます。そして、取り込まれた擬似的な周波数は、ほとんど生き残りません。一方、非箱型ウインドウは、メインロブを広げることになり、その結果、解像度を全体的に減少させます。

関数 periodogram は、データに適用するウインドウを指定することにより、修正されたピリオドグラムを計算することができます。たとえば、箱型ウインドウと Hamming ウインドウを比べてみましょう。

randn('state',0)
fs = 1000;           % サンプリング周波数
t = (0:fs/10)./fs;   % 1秒の間にサンプルを分布
A = [1 2];           % 正弦波の振幅
f = [150;140];       % 正弦波の周波数
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
periodogram(xn,boxcar(length(xn)),1024,fs);

periodogram(xn,hamming(length(xn)),1024,fs);

Hamming ウインドウを適用したピリオドグラムの中で、サイドロブが小さくはなっていますが、2つのメインピークはより広がっていることはわかります。実際に、Hamming ウインドウに対応するメインロブの 3 dB 幅は、箱型ウインドウのものの約2倍になっています。そのために、固定したデータ長に対して、PSD 解像度が Hamming ウインドウにより達成できるものは、箱型ウインドウで達成できるものの約半分になります。メインロブの幅とサイドロブの高さとの関係は、Kaiser ウインドウのような可変ウインドウを使うことにより、拡大することができます。

非箱型ウインドウは、信号の平均パワーに影響を与えます。これは、時間サンプルは、ウインドウを適用することにより、その一部が減衰するためです。これに対する補償を行うには、関数 periodogram を使って、ウインドウの平均パワーが1になるように正規化することです。ウインドウの選択は、信号の平均パワーに影響を与えません。

PSD の修正ピリオドグラム推定は、つぎのようになります。

ここで、U は、ウインドウの正規化定数です。

これは、選択するウインドウには無関係です。正規化定数のような U の加算は、修正ピリオドグラムが漸近的にバイアスが適用されていないことを補償します。

Welch 法

PSD の改良した推定は、Welch [8] により提唱されました。この方法は、時系列データを(重ね合わせ可能)セグメント毎に分割し、各セグメントについて修正ピリオドグラムを計算し、その後、PSD 推定の平均を計算するものです。この結果を、Welch の PSD 推定と言います。

Welch 法は、Signal Processing Toolbox の中の関数 pwelch で実現できます。デフォルトでは、データは、50 % の重ね合わせを含んで、8つのセグメントに分割されます。各セグメントの修正ピリオドグラムを計算するときに、Hamming ウインドウを使います。

修正ピリオドグラムの平均化は、データ全体を使った一つのピリオドグラムの場合と比べて、推定の分散を小さくする傾向があります。しかし、セグメント間の重ね合わせは、冗長な情報を導入する傾向があり、この影響は、非箱型ウインドウを使うことにより軽減することができます。この箱型ウインドウは、セグメントの端のサンプル(重ね合わせサンプル)に設定する重み、または、重要性を軽減するものです。

しかし、上述したように、短いデータと非箱型ウインドウとの組み合わせは、推定子の解像度を低下させる結果になります。まとめると、分散の低減化と解像度の間にはトレードオフがあります。ピリオドグラムに関する改良された推定を得るように、Welch 法のパラメータを取り扱うことができます。特に、SNR が低い場合は、より可能です。このことを、つぎの例題の中で示すことができます。

301サンプルのオリジナル信号を考えましょう。

randn('state',1)
fs = 1000;               % サンプリング周波数
t = (0:0.3*fs)./fs;      % 301 サンプル
A = [2 8];               % 正弦波の振幅 (行ベクトル)
f = [150;140];           % 正弦波の周波数 (列ベクトル)
xn = A*sin(2*pi*f*t) + 5*randn(size(t));
periodogram(xn,boxcar(length(xn)),1024,fs);

50 % の重ね合わせ部をもち、3つのセグメントに分割したデータについての Welch スペクトル推定を行います。

pwelch(xn,boxcar(150),75,512,fs);

上のピリオドグラムの中で、雑音とリーケージは、正弦波によるピークと人工的に作られたピーク(雑音)との区別を不可能にしています。対照的に、Welch 法により作成された PSD は、幅の広いピークを示しますが、2つの正弦波の区別は可能です。これらは、雑音レベルとは、かなり異なるものです。

しかし、さらに分散を低減しようとする場合、解像度の低下が、正弦波のどちらかを検出できないことになります。

pwelch(xn,hamming(100),75,512,fs);

PSD 推定の Welch 法のより詳細な議論については、Kay [2]と Welch [8] を参照してください。

Welch 法でのバイアスと正規化

Welch 法は、PSD のバイアス付きの推定子を提供します。期待値は、つぎのようになります。

ここで、L_s は、データセグメントの長さで、U は、修正ピリオドグラムの定義で使った正規化定数です。すべてのピリオドグラムの場合と同じく、Welch 推定子は、漸近的にバイアスのないものになります。固定された長さのデータに対して、Welch 推定のバイアスは、L_s < L であるために、ピリオドグラムのバイアスよりも大きくなります。

Welch 推定子の分散は、使用するウインドウとセグメント間の重ね合わせ量の二つに依存するので、計算することが困難です。基本的に、分散は、平均化に使用する修正ピリオドグラムのセグメント数に反比例します。

Multitaper 法

ピリオドグラムは、長さ L の信号 x_L[n] をL 個の FIR バンドパスフィルタのフィルタバンク(並列配置のフィルタ群)を通したものと解釈できます。これらのバンドパスフィルタの各々の 3 dB バンド幅は、近似的に、f_s / L と等価であることが示されます。これらのバンドパスフィルタの各々の応答は、スペクトルリーケージで議論した箱型ウインドウの応答に似ています。ピリオドグラムは、フィルタリングされた(バンドパスフィルタの出力)信号の1つのサンプルだけを使って、個々のパワーの計算したものと見なすことができます。そして、 x_L[n] の PSD は、各バンドパスフィルタのバンド幅に渡り定数であると仮定しています。

信号の長さを長くすればするほど、各バンドパスフィルタのバンド幅は小さくなり、より望ましいフィルタになり、このフィルタのバンド幅に渡って、定数の PSD の近似を改良することができます。これは、信号の長さを長くするに連れ、ピリオドグラムの PSD 推定が良くなる理由の別な解釈を与えます。しかし、ピリオドグラムの推定の精度を良くする観点において、2つのファクタが明らかになります。まず、箱型ウインドウは、質の悪いバンドパスフィルタです。2つ目は、各バンドパスフィルタの出力でのパワーの計算は、出力信号の単一サンプルをもとにしていて、その結果、非常に、粗い近似になります。

Welch 法は、フィルタバンクを使って、同様な解釈ができます。Welch 法の実現で、いくつかのサンプルが、出力パワーを計算するために使われ、その結果、推定の分散を低減します。一方、各バンドパスフィルタのバンド幅は、ピリオドグラム法に対応するものよりも広くなり、結果として、解像度が低下します。フィルタバンクモデルは、分散と解像度の関係を向上させる新しい解釈を提供するものです。

Thompson による multitaper 法 (MTM) は、これらの結果を使って、改良された PSD 推定を提供するものです。本質的に(ピリオドグラム法で使用した)箱型ウインドウであるバンドパスフィルタを使用する代わりに、MTM 法は、最適なバンドパスフィルタのバンクを使って、推定を計算します。これらの最適 FIR フィルタは、 discrete prolate spheroidal sequences (DPSSs, Slepian sequencesとしても知られています) として、知られた数列群から導出されたものです。

加えて、MTM 法は、分散と解像度とのバランスを考慮し、時間-バンド幅パラメータを使います。このパラメータは、時間-バンド幅の積 NW で与えられ、スペクトルを計算するために使用するテーパの数に直接関連します。推定を行う場合に使用するものは、常に、2*NW-1 です。このことは、NW が増加するに連れ、パワースペクトルの推定が良くなり、推定の分散が減少することを意味しています。しかし、各テーパのバンド幅は、NW に比例します。それで、NW が増加すると、各推定は、スペクトルのリーケージをより鮮明にし、スペクトル推定全体に、バイアスがより適用されます。各データセットに対して、バイスと分散との最適トレードオフを可能にする値が、NW に対して存在します。

MTM 法を実現する Signal Processing Toolbox の関数は、pmtm です。pmtm を使って、前の例題の xn のPSD を計算しましょう。

randn('state',0)
fs = 1000;               % サンプリング周波数
t = (0:fs)/fs;           % 1秒の間にサンプルを分布
A = [1 2];               % 正弦波の振幅
f = [150;140];           % 正弦波の周波数
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
[P,F] = pmtm(xn,4,1024,fs);
plot(F,10*log10(P))          % dB/Hz 単位でプロット
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');

時間-バンド幅の積を小さくすることにより、より大きな分散を押さえて、解像度を増すことができます。

[P1,f] = pmtm(xn,3/2,1024,fs);
plot(f,10*log10(P1))         % dB/Hz 単位でプロット
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');

両方共に、平均パワーは保存されることに注意してください。

Pow = (fs/1024) * sum(P)
Pow =
    2.4926
Pow1 = (fs/1024) * sum(P1)
Pow1 =
    2.4927

この方法は、離散扁長回転楕円体列(Slepian列)を計算するために、Welch 法より計算時間を必要とします。長いデータ列(10,000点以上)に対して、DPSS を一回計算し、それを MAT-ファイルとして保存しておくことが重要です。M-ファイル dpsssave, dpssload, dpssdir と dpssclear は、MAT-ファイル dpss.mat でセーブした DSPP のデータベースを使用するようにしています。

クロススペクトル密度関数

PSD は、信号 x_n と y_n を使って、つぎのように定義される クロススペクトル密度 (CSD) 関数の特殊なものです。

相関列、共分散列の場合と同じように、ツールボックスは、信号列が有限なので、PSD と CSD を計算します。

2つの同じ長さの信号 x と y のクロススペクトル密度を、Welch 法を使って、推定するには、関数 csd を使い、x の FFT と y の FFT の共役部との積として、ピリオドグラムを作成します。実数値 PSD と異なり、CSD は複素関数です。csd は、関数 pwelch と同様に、x と y の一部分やウインドウと適用した部分を取り扱います。

Sxy = csd(x,y,nfft,fs,window,numoverlap)

信頼区間

付加的な入力引数 p に、信頼区間のパーセントを指定し、引数 numoverlap に0を設定する型で、関数 csd を使うことで、信頼区間を計算することができます。

[Sxy,Sxyc,f] = csd(x,y,nfft,fs,window,0,p)

p は、0 と 1 の間のスカラ値です。この関数は、ウインドウが適用された重ね合わせ部をもたないカイ二乗分布のピリオドグラムと仮定し、信頼区間を計算しています。この仮定は、信号がガウス分布したランダムプロセスの場合に正しいものです。これらの仮定が正しい場合、信頼区間は、つぎのようになります。

[Sxy-Sxyc(:,1) Sxy+Sxyc(:,2)]

は、確率 p をもつ真の CSD をカバーしています。引数 numoverlap に 0 以外の数値を設定した場合、ワーニングが生じ、対象部分に重ね合わせがあり、算出した信頼区間は、信頼性が低いことを示します。

伝達関数の推定

Welch 法の一つのアプリケーションは、ノンパラメトリックシステム同定です。H が線形で、時不変システムで、x(n) と y(n) は、H への入力と、出力と仮定します。そして、x(n) のパワースペクトルは、つぎの式により、x(n) と y(n) の CSD に関連付けられます。

x(n) と y(n) の間の伝達関数の推定は、つぎのようになります。

この方法は、大きさと位相の両方を推定します。関数 tfe は、Welch 法を使って、CSD とパワースペクトルを計算し、伝達関数推定用に、それらの比を形作ります。tfe を関数 csd と同様に使用することができます。

信号 xn に、FIR フィルタを適用し、実際の大きさの応答と推定応答をプロットします。

h = ones(1,10)/10;                % 移動平均フィルタ
yn = filter(h,1,xn);
[HEST,f] = tfe(xn,yn,256,fs,256,128,'none');
H = freqz(h,1,f,fs);
subplot(2,1,1); plot(f,abs(H)); 
title('Actual Transfer Function Magnitude'); 
subplot(2,1,2); plot(f,abs(HEST));
title('Transfer Function Magnitude Estimate'); 
xlabel('Frequency (Hz)');

コヒーレンス関数

2つの信号 x(n) と y(n) の間の大きさの二乗のコヒーレンスは、つぎのように表わせます。

この割合を使った表現は、0と1の間の実数値になり、周波数での x(n) と y(n) の間の相関の尺度になります。

関数 cohere は、データ列 x と y を使って、個々のパワースペクトルと CSD を計算し、CSD の大きさの二乗と個々のパワースペクトルの積との比を出力します。このオプションと演算は、関数 csd と tfe と同じです。

xn とフィルタ出力 yn の周波数に関するコヒーレンス関数は、つぎのようにして求めます。

cohere(xn,yn,256,fs,256,128,'none')

入力データ列の長さ nfft、ウインドウの長さ window とし、cohere 演算では、単一レコードしか演算しないように、一つのウインドウ内での重ね合わせ部分のデータ点数 numoverlap を設定すると、関数は、すべて1からなるものを出力します。これは、線形的な従属性をもつデータに対するコヒーレンス関数は、1になるためです。

スペクトル推定法の概要

パラメトリック法