線形予測と自己回帰モデリング

ライブスクリプトを開く

この例は、自己回帰モデリングと線形予測の関係を比較する方法を示します。線形予測と自己回帰モデリングという 2 つの問題は、同じ数値結果を得ることが可能です。両方の問題の最終的な目的は、線形フィルターのパラメーターを確定することです。しかし、それぞれの問題で使われるフィルターは異なります。

はじめに

線形予測の目的は、過去のサンプルの線形結合に基づいて、自己回帰過程の今後のサンプルを的確に予測できる FIR フィルターを決定することです。実際の自己回帰信号と予測信号の差は、予測誤差と呼ばれます。理想的には、この誤差がホワイトノイズです。

自己回帰モデリングの目的は、ホワイトノイズによって励起されたとき、モデル化する自己回帰過程と同じ統計値をもつ信号を生成する全極 IIR フィルターを決定することです。

入力がホワイトノイズである全極フィルターを使用した AR 信号の生成

LPC 関数と FIR フィルターを使用して、目的の自己回帰信号を作成するために使用するパラメーターを生成します。ここでは FIR1 および LPC の使用は重要ではありません。たとえば、d を [1 1/2 1/3 1/4 1/5 1/6 1/7 1/8] などの単純な数値で置き換え、p0 を 1e-6 などに置き換えることができますが、フィルターの形の方が扱いやすいため、これを使用します。

b = fir1(1024, .5);
[d,p0] = lpc(b,7);

自己回帰信号を生成するには、分散 p0 のホワイトガウスノイズを使用して全極フィルターを励起します。分散 p0 を取得するには、ノイズ発生器で 'gain' 項として SQRT(p0) を使用しなければなりません。

rng(0,'twister'); % Allow reproduction of exact experiment
u = sqrt(p0)*randn(8192,1); % White Gaussian noise with variance p0

ここで、AR 信号を生成するため、ホワイトガウスノイズ信号と全極フィルターを使用します。

x = filter(1,d,u);

ユール・ウォーカー法を使用した信号からの AR モデルの検出

ユール・ウォーカー式を解くことで、指定された信号 x の統計値と一致する統計値をもつ AR 信号をホワイトノイズによって励起されたときに生成する全極フィルターのパラメーターを決定できます。これは、自己回帰モデリングと呼ばれます。ユール・ウォーカー式を解くには、自己相関関数 x を推定する必要があります。レビンソンアルゴリズムを使用すると、効率的にユール・ウォーカー式を解くことができます。ここでは、関数 ARYULE を使用します。

[d1,p1] = aryule(x,7);

AR モデルと AR 信号の比較

次に、AR 信号 x をモデル化するために使用した全極フィルターの周波数応答を計算します。フィルターがホワイトガウスノイズによって励起されたとき、周波数応答の振幅の 2 乗にホワイトノイズ入力の分散を乗算すると、このフィルターの出力のパワースペクトル密度が得られることはよく知られています。この出力パワースペクトル密度を計算する 1 つの方法として、関数 freqz を使用します。

[H1,w1] = freqz(sqrt(p1),d1);

自己回帰信号 x のモデル化がどの程度適切かを確認するため、関数 freqz を使用して計算したモデルの出力のパワースペクトル密度と、PERIODOGRAM を使用して計算した x のパワースペクトル密度の推定値を重ね合わせます。ピリオドグラムは 2*pi でスケーリングされ、片側であることに注意してください。比較するにはピリオドグラムを調整する必要があります。

periodogram(x)
hold on
hp = plot(w1/pi,20*log10(2*abs(H1)/(2*pi)),'r'); % Scale to make one-sided PSD
hp.LineWidth = 2;
xlabel('Normalized frequency (\times \pi rad/sample)')
ylabel('One-sided PSD (dB/rad/sample)')
legend('PSD estimate of x','PSD of model output')

Figure contains an axes object. The axes object with title Periodogram Power Spectral Density Estimate contains 2 objects of type line. These objects represent PSD estimate of x, PSD of model output.

LPC を使用した線形予測の実行

線形予測の問題に戻ります。まず FIR 予測フィルターを決定します。このためには LPC を使用しますが、LPC の結果は若干の解釈が必要です。LPC はホワイトニングフィルター A(z) 全体の係数を返します。このフィルターは、自己回帰信号 x を入力として受け取り、予測誤差を出力として返します。ただし、A(z) の予測フィルターは形式 B(z) = 1- A(z) で組み込まれています (ここで、B(z) が予測フィルター)。LPC で計算した係数と誤差分散は、実質的には関数 ARYULE で計算した係数と誤差分散と同じですが、その解釈が異なります。

[d2,p2] = lpc(x,7);
[d1.',d2.']

ans = 8×2

    1.0000    1.0000
   -3.5245   -3.5245
    6.9470    6.9470
   -9.2899   -9.2899
    8.9224    8.9224
   -6.1349   -6.1349
    2.8299    2.8299
   -0.6997   -0.6997

上記の説明のように A(z) から B(z) を抽出し、FIR 線形予測子フィルターを使用して、過去の値の線形結合に基づき、自己回帰信号の今後の値の推定値を取得します。

xh = filter(-d2(2:end),1,x);

実際の信号と予測信号の比較

7 タップ FIR 予測フィルターを使用して実行した内容を把握するため、元の自己回帰信号と、線形予測子から得られた信号推定値をプロットします (200 サンプル)。ここで、予測フィルターの 1 サンプル遅延を考慮します。

cla
stem([x(2:end),xh(1:end-1)])
xlabel('Sample time')
ylabel('Signal value')
legend('Original autoregressive signal','Signal estimate from linear predictor')
axis([0 200 -0.08 0.1])

Figure contains an axes object. The axes object with title Periodogram Power Spectral Density Estimate contains 2 objects of type stem. These objects represent Original autoregressive signal, Signal estimate from linear predictor.

予測誤差の比較

予測誤差パワー (分散) は、LPC の 2 番目の出力として返されます。その値は、理論的には AR モデリングの問題で全極フィルターを励起するホワイトノイズの分散 (p1) と同じです。予測誤差自体からも、この分散を推定できます。

p3 = norm(x(2:end)-xh(1:end-1),2)^2/(length(x)-1);

以下のすべての値は理論的には同じです。違いは、これに含まれるさまざまな計算および近似誤差が原因です。

[p0 p1 p2 p3]

ans = 1×4
10^-5 ×

    0.5127    0.5305    0.5305    0.5068

参考

aryule | lpc