ノンパラメトリック法

本節以降では、ノンパラメトリック推定法のピリオドグラム、修正ピリオドグラム、ウェルチ法、およびマルチテーパー法を、関連する CPSD 関数、伝達関数推定、およびコヒーレンス関数とともに説明します。

ピリオドグラム

一般的に、変動過程の PSD を推定する 1 つの方法は、過程のサンプルの離散フーリエ変換を求め (通常は FFT を使用してグリッド上で実行)、結果の振幅の 2 乗を適宜スケーリングするというものです。この推定は、"ピリオドグラム" と呼ばれます。

長さ L を持つ信号 $x_{L} (n)$ の PSD のピリオドグラム推定は、以下のようになります。

$P_{x x} (f) = \frac{1}{L F_{s}} {| \sum_{n = 0}^{L - 1} x_{L} (n) e^{- j 2 π f n / F_{s}} |}^{2},$

ここで Fs はサンプリング周波数です。

実際には、 $P_{x x} (f)$ の計算は有限個の周波数点のみで実行でき、この計算には通常 FFT が使用されます。ピリオドグラム法を実装する場合は通常、次の周波数での $N$ 点の PSD 推定を計算します。

$f_{k} = \frac{k F_{s}}{N}, k = 0, 1, \dots, N - 1 .$

場合によっては、FFT アルゴリズムによるピリオドグラムの計算は、周波数点の数が 2 のべき乗である場合に効率的です。したがって、入力信号に 0 を加え、長さが 2 のべき乗となるように伸ばす操作はめずらしくありません。

ピリオドグラムの例として、以下の 1001 要素の信号 xn を考えます。この中には、2 つの正弦波とノイズが含まれています。

fs = 1000;                % Sampling frequency
t = (0:fs)/fs;            % One second worth of samples
A = [1 2];                % Sinusoid amplitudes (row vector)
f = [150;140];            % Sinusoid frequencies (column vector)
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
% The three last lines are equivalent to
% xn = sin(2*pi*150*t) + 2*sin(2*pi*140*t) + 0.1*randn(size(t));

PSD のピリオドグラム推定は、periodogram を使用して計算できます。この場合、データベクトルにハミングウィンドウが乗算され、修正されたピリオドグラムが作成されます。

[Pxx,F] = periodogram(xn,hamming(length(xn)),length(xn),fs);
plot(F,10*log10(Pxx))
xlabel('Hz')
ylabel('dB')
title('Modified Periodogram Power Spectral Density Estimate')

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

アルゴリズム

ピリオドグラムでは、以下のように FFT の出力の計算とスケーリングを行い、パワーと周波数プロットが生成されます。

入力信号が実数値の場合、結果として得られる FFT の振幅は、ゼロ周波数 (DC) に対し対称的となります。偶数長の FFT では、最初の (1+nfft/2) 点だけが一意です。一意である値の数を決定し、一意な点だけを保存します。
一意の FFT 値の振幅を 2 乗します。振幅の 2 乗を (DC の場合を除いて) $2 / (F_{s} N)$ でスケーリングします。N は、ゼロパディング前の信号の長さです。DC の値は、 $1 / (F_{s} N)$ によりスケーリングします。
一意な点の数、nfft とサンプリング周波数から周波数ベクトルを作成します。
周波数に対する振幅の二乗 FFT をプロットします。

ピリオドグラムの性能

ライブスクリプトを開く

本節以降では、スペクトル漏れ、解像度、バイアスおよび分散を基準にピリオドグラムの性能について説明します。

スペクトル漏れ

有限長 (長さ $L$ ) の信号 $x_{L} (n)$ の PSD を考えます。大半の場合、 $x_{L} (n)$ は、無限長の信号 $x (n)$ に有限長の箱型ウィンドウ $w_{R} (n)$ を乗算した結果と解釈すると便利です。

$x_{L} (n) = x (n) w_{R} (n) .$

時間領域での乗算は、周波数領域での畳み込みになるので、周波数領域でのピリオドグラムの期待値は、次のようになります。

$E {{P_{}^{ˆ}}_{x x} (f)} = \frac{1}{F_{s}} \int_{- F_{s} / 2}^{F_{s} / 2} \frac{\sin^{2} (L π (f - f^{'}) / F_{s})}{L \sin^{2} (π (f - f^{'}) / F_{s})} P_{x x} (f^{'}) d f^{'},$

これにより、ピリオドグラムの期待値は、真の PSD とディリクレカーネルの 2 乗との畳み込みであることが示されます。

畳み込みの影響を最も理解しやすいのは、正弦波データの場合です。 $x (n)$ が $M$ 個の複素正弦波の和から構成されていると仮定します。

$x (n) = \sum_{k = 1}^{N} A_{k} e^{j ω_{k} n} .$

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

$X (ω) = \sum_{k = 1}^{N} A_{k} δ (ω - ω_{k}),$

有限長のシーケンスについては、以下のように表せます。

$X (ω) = \int_{- π}^{π} (\sum_{k = 1}^{N} A_{k} δ (ε - ω_{k})) W_{R} (ω - ε) d ε .$

この方程式は以下と等価です。

$X (ω) = \sum_{k = 1}^{N} A_{k} W_{R} (ω - ω_{k}) .$

有限長の信号のスペクトルで、ディラックデルタ関数は $W_{R} (ω - ω_{k})$ の型の項で置き換えられました。これは、周波数 $ω_{k}$ を中心とする箱型ウィンドウの周波数応答に対応します。

箱型ウィンドウの周波数応答は、以下に示す周期的 sinc の形状をもっています。

L = 32;
[h,w] = freqz(rectwin(L)/L,1);
y = diric(w,L);

plot(w/pi,20*log10(abs(h)))
hold on
plot(w/pi,20*log10(abs(y)),'--')
hold off
ylim([-40,0])
legend('Frequency Response','Periodic Sinc')
xlabel('\omega / \pi')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Frequency Response, Periodic Sinc.

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

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

fs = 1000;                 % Sampling frequency
t = (0:fs/10)/fs;          % One-tenth second worth of samples
A = [1 2];                 % Sinusoid amplitudes
f = [150;140];             % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
periodogram(xn,rectwin(length(xn)),1024,fs)

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

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

解像度

"解像度" は、個々のスペクトルを区別できる程度 (能力) を表し、スペクトル推定の性能解析の基本的な考えになっています。

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

言い換えれば、周波数 $f_{1}$ と $f_{2}$ の 2 つの正弦波の分解可能性条件には次が必要です。

$f_{2} - f_{1} > \frac{F_{s}}{L} .$

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

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

fs = 1000;                  % Sampling frequency
t = (0:fs/15)/fs;           % 67 samples
A = [1 2];                  % Sinusoid amplitudes
f = [150;140];              % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
periodogram(xn,rectwin(length(xn)),1024,fs)

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

解像度に関する上の議論は、S/N 比 (SNR) が、相対的にかなり高いためノイズの影響を考えなくても良いです。S/N 比が低い場合、真のスペクトルは、分離が難しく、ノイズにより、ピリオドグラム上でのスペクトルの計算に人工的な邪魔な部分が現れます。例として、次のものを参考としてください。

fs = 1000;                  % Sampling frequency
t = (0:fs/10)/fs;           % One-tenth second worth of samples
A = [1 2];                  % Sinusoid amplitudes
f = [150;140];              % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + 2*randn(size(t));
periodogram(xn,rectwin(length(xn)),1024,fs)

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

ピリオドグラムのバイアス

ピリオドグラムは PSD のバイアス付き推定子です。その期待値は以下のようになることが先に示されました。

$E {{P_{}^{ˆ}}_{x x} (f)} = \frac{1}{F_{s}} \int_{- F_{s} / 2}^{F_{s} / 2} \frac{\sin^{2} (L π (f - f^{'}) / F_{s})}{L \sin^{2} (π (f - f^{'}) / F_{s})} P_{x x} (f^{'}) d f^{'} .$

ピリオドグラムは漸近的に無バイアスで、このことは、データ長が無限大に近付くにつれ、箱型ウィンドウの周波数応答がディラックデルタ関数により近くなるという、前述の観測からも明らかです。しかし場合によっては、ピリオドグラムは、データが長い場合でも PSD の推定法を満たさないものとなります。これは、次に説明する、ピリオドグラムの分散によるものです。

ピリオドグラムの分散

ピリオドグラムの分散は、次のように示されます。

$V a r ({P_{}^{ˆ}}_{x x} (f)) = {\begin{array}{c} P_{x x}^{2} (f), & 0 < f < F_{s} / 2, \\ 2 P_{x x}^{2} (f), & f = 0, F_{s} / 2, \end{array}$

これは、データ長 $L$ を無限大にするにつれ、分散は、ゼロ方向に向かないことを示しています。統計上は、ピリオドグラムは PSD の一致推定器ではありません。しかし、SNR が高い場合のスペクトル推定では、特に長いデータレコードの場合、有効なツールになり得ます。

修正ピリオドグラム

ライブスクリプトを開く

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

periodogram では、データに適用するウィンドウを指定することにより、修正されたピリオドグラムを計算できます。たとえば、既定の箱型ウィンドウとハミングウィンドウを比べます。どちらの場合も、同じ数の DFT 点を指定します。

fs = 1000;                  % Sampling frequency
t = (0:fs/10)/fs;           % One-tenth second worth of samples
A = [1 2];                  % Sinusoid amplitudes
f = [150;140];              % Sinusoid frequencies
nfft = 1024;

xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
periodogram(xn,rectwin(length(xn)),nfft,fs)

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

periodogram(xn,hamming(length(xn)),nfft,fs)

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

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

非箱型ウィンドウ処理は、信号の平均パワーに影響を与えます。これは、時間サンプルは、ウィンドウを適用することにより、その一部が減衰するためです。これを補うため、periodogram と pwelch は、ウィンドウを正規化し、1 の平均パワーをもつようにします。ウィンドウの選択は、信号の平均パワーに影響を与えません。周波数成分が PSD 推定子によって解決されていない場合は、ウィンドウの選択は平均パワーに影響を与えます。

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

${P_{}^{ˆ}}_{x x} (f) = \frac{{| X (f) |}^{2}}{F_{s} L U},$

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

$U = \frac{1}{L} \sum_{n = 0}^{N - 1} {| w (n) |}^{2} .$

L が大きい値の場合、U はウィンドウ長に依存しなくなる傾向があります。U を正規化定数として加算することで、修正ピリオドグラムが漸近的に非バイアスとなるようにします。

ウェルチ法

ライブスクリプトを開く

改善された PSD の推定子としては、ウェルチにより提唱されたものがあります。この方法は、時系列データを (重ね合わせ可能) セグメントごとに分割し、各セグメントについて修正ピリオドグラムを計算し、その後、PSD 推定の平均を計算するものです。この結果を、ウェルチの PSD 推定と言います。ツールボックス関数 pwelch は、ウェルチ法を実装しています。

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

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

301 サンプルからなる信号について考えます。

fs = 1000;             % Sampling frequency
t = (0:0.3*fs)/fs;     % 301 samples
A = [2 8];             % Sinusoid amplitudes (row vector)
f = [150;140];         % Sinusoid frequencies (column vector)

xn = A*sin(2*pi*f*t) + 5*randn(size(t));
periodogram(xn,rectwin(length(xn)),1024,fs)

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

箱型ウィンドウを使用して、3 つのセグメントに分割された 50% のオーバーラップをもつデータに対するウェルチスペクトル推定を行います。

pwelch(xn,rectwin(150),50,512,fs)

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

上のピリオドグラムの中で、ノイズと漏れは、正弦波によるピークと人工的に作られたピーク (ノイズ) との区別を不可能にしています。対照的に、ウェルチ法により作成された PSD は幅の広いピークを示しますが、2 つの正弦波の区別は可能です。これらは、"ノイズフロア" とはかなり異なるものです。

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

pwelch(xn,rectwin(100),75,512,fs)

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

ウェルチ法でのバイアスと正規化

ウェルチ法は PSD のバイアス付き推定子を提供します。PSD 推定の期待値は、次のようになります。

$E {P_{Welch} (f)} = \frac{1}{F_{s} L U} \int_{- F_{s} / 2}^{F_{s} / 2} {| W (f - f^{'}) |}^{2} P_{x x} (f^{'}) d f^{'},$

ここで、L はデータセグメントの長さ、U は修正ピリオドグラムの定義のものと同じ正規化定数、W(f) はウィンドウ関数のフーリエ変換です。すべてのピリオドグラムの場合と同じく、ウェルチ推定子は、漸近的にバイアスのないものになります。固定された長さのデータに対し、ウェルチ推定のバイアスは、セグメントの長さが全データサンプルの長さより短いため、ピリオドグラムのバイアスよりも大きくなります。

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

マルチテーパー法

ライブスクリプトを開く

ピリオドグラムは、長さ $L$ の信号 $x_{L} (n)$ を $L$ 個の FIR バンドパスフィルターから成るフィルターバンク (並列配置のフィルター群) に通したものと解釈できます。これらのバンドパスフィルターのそれぞれの 3 dB 帯域幅は、近似的に $f_{s} / L$ と等価であることが示されます。これらのバンドパスフィルターのそれぞれの振幅応答は箱型ウィンドウの応答に似ています。ピリオドグラムは、フィルター処理された信号の 1 つのサンプルのみを使用して、フィルターが適用された個々の信号のパワー (各バンドパスフィルターの出力) を計算したものと見なすことができます。そして、 $x_{L} (n)$ の PSD は、各バンドパスフィルターの帯域幅にわたり定数であると仮定しています。

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

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

トンプソンの "マルチテーパー法" (MTM) はこれらの結果を使用して、改良された PSD 推定を行います。本質的に箱型ウィンドウであるバンドパスフィルターを (ピリオドグラム法の場合のように) 使用する代わりに、MTM 法は、一群の最適なバンドパスフィルターを使用して、推定を計算します。これらの最適 FIR フィルターは、"離散扁長回転楕円体列" (DPSS、"スレピアン列" とも呼ばれる) として知られるシーケンス群から導出されたものです。

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

MTM 法を実装する Signal Processing Toolbox™ の関数は pmtm です。信号の PSD を計算するには pmtm を使用します。

fs = 1000;                % Sampling frequency
t = (0:fs)/fs;            % One second worth of samples
A = [1 2];                % Sinusoid amplitudes
f = [150;140];            % Sinusoid frequencies

xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
pmtm(xn,4,[],fs)

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

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

pmtm(xn,1.5,[],fs)

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

この方法は、離散扁長回転楕円体列を計算するため、ウェルチ法と比べてより長い計算時間を要します。長いデータ列 (10,000 点以上) の場合、DPSS を 1 回計算し、それらを MAT ファイルとして保存しておくと便利です。dpsssave、dpssload、dpssdir および dpssclear が、保存した DPSS のデータベースを MAT ファイル dpss.mat 内に保持するために用意されています。

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

PSD は、"クロススペクトル密度" (CPSD) 関数の特殊なものであり、2 つの信号 x(n) と y(n) との間で次のように定義されます。

$P_{x y} (ω) = \frac{1}{2 π} \sum_{m = - \infty}^{\infty} R_{x y} (m) e^{- j ω m} .$

相関列および共分散列の場合と同じように、ツールボックスでは、信号列が有限であることから PSD と CPSD を "推定" します。

ウェルチ法を使用して 2 つの同じ長さの信号 x および y のクロススペクトル密度を推定するには、関数 cpsd を使い、x の FFT と y の FFT の共役の積としてピリオドグラムを作成します。実数値 PSD と異なり、CPSD は複素関数です。cpsd は、関数 pwelch と同様に、x と y の分割したものやウィンドウ処理した部分を取り扱います。

Sxy = cpsd(x,y,nwin,noverlap,nfft,fs)

伝達関数の推定

ライブスクリプトを開く

ウェルチ法の使用方法の 1 つとして、ノンパラメトリックシステム同定があります。H は線形の時不変システムで、x(n) および y(n) は、それぞれ H への入力と出力であると仮定します。この場合、x(n) のパワースペクトルは、以下の式により、x(n) および y(n) の CPSD に関連付けられます。

$P_{y x} (ω) = H (ω) P_{x x} (ω) .$

x(n) と y(n) の間の伝達関数の推定値は以下のようになります。

$H_{}^{ˆ} (ω) = \frac{{P_{}^{ˆ}}_{y x} (ω)}{{P_{}^{ˆ}}_{x x} (ω)} .$

この方法は、振幅と位相の両方の情報を推定します。関数 tfestimate は、ウェルチ法を使用して CPSD とパワースペクトルを計算し、伝達関数推定用にそれらの商を生成します。tfestimate を関数 cpsd と同様に使用できます。

ホワイトガウスノイズに含まれる 2 つの正弦波で構成される信号を生成します。

rng('default')

fs = 1000;                % Sampling frequency
t = (0:fs)/fs;            % One second worth of samples
A = [1 2];                % Sinusoid amplitudes
f = [150;140];            % Sinusoid frequencies

xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));

信号 xn を FIR 移動平均フィルターでフィルター処理します。実際の振幅応答と推定応答を計算します。

h = ones(1,10)/10;	            % Moving-average filter
yn = filter(h,1,xn);

[HEST,f] = tfestimate(xn,yn,256,128,256,fs);
H = freqz(h,1,f,fs);

結果をプロットします。

subplot(2,1,1)
plot(f,abs(H))
title('Actual Transfer Function Magnitude')
yl = ylim;
grid
subplot(2,1,2)
plot(f,abs(HEST))
title('Transfer Function Magnitude Estimate')
xlabel('Frequency (Hz)')
ylim(yl)
grid

Figure contains 2 axes objects. Axes object 1 with title Actual Transfer Function Magnitude contains an object of type line. Axes object 2 with title Transfer Function Magnitude Estimate contains an object of type line.

コヒーレンス関数

ライブスクリプトを開く

2 つの信号 x(n) と y(n) の間の振幅二乗コヒーレンスは、以下のようになります。

$C_{x y} (ω) = \frac{{| P_{x y} (ω) |}^{2}}{P_{x x} (ω) P_{y y} (ω)} .$

この商は 0 と 1 の間の実数値になり、周波数 $ω$ での x(n) と y(n) の間の相関度を測定します。

関数 mscohere は、シーケンス xn と yn を使用して個々のパワースペクトルと CPSD を計算し、CPSD の振幅二乗と個々のパワースペクトルの積との比を返します。このオプションと演算は、関数 cpsd および tfestimate と同じです。

ホワイトガウスノイズに含まれる 2 つの正弦波で構成される信号を生成します。信号は 1 kHz で 1 秒間サンプリングされます。

rng('default')

fs = 1000;
t = (0:fs)/fs;
A = [1 2];                % Sinusoid amplitudes
f = [150;140];            % Sinusoid frequencies

xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));

信号 xn を FIR 移動平均フィルターでフィルター処理します。xn のコヒーレンス関数およびフィルター出力 yn を周波数関数として計算してプロットします。

h = ones(1,10)/10;
yn = filter(h,1,xn);

mscohere(xn,yn,256,128,256,fs)

Figure contains an axes object. The axes object with title Coherence Estimate via Welch contains an object of type line.

1 つのウィンドウ内の入力シーケンスの長さ、ウィンドウの長さ、オーバーラップ部分のデータ点数が、mscohere が単一のレコードに対してのみ演算を行うようなものである場合、関数はすべて 1 を返します。これは、線形的な従属性をもつデータに対するコヒーレンス関数は 1 になるためです。