rlevinson

逆レビンソン・ダービン再帰法

構文

r = rlevinson(a,efinal) [r,u] = rlevinson(a,efinal) [r,u,k] = rlevinson(a,efinal) [r,u,k,e] = rlevinson(a,efinal)

説明

逆レビンソン・ダービン再帰法は、r についての次の対称テプリッツ線形方程式系を解く、下降演算アルゴリズムを実行します。ここで、r = [r(1) … r(p + 1)] および r(i)^* は、r(i) の複素共役を表します。

$[\begin{matrix} r (1) & r {(2)}^{*} & \dots & r {(p)}^{*} \\ r (2) & r (1) & \dots & r {(p - 1)}^{*} \\ ⋮ & ⋱ & ⋱ & ⋮ \\ r (p) & \dots & r (2) & r (1) \end{matrix}] [\begin{matrix} a (2) \\ a (3) \\ ⋮ \\ a (p + 1) \end{matrix}] = [\begin{matrix} - r (2) \\ - r (3) \\ ⋮ \\ - r (p + 1) \end{matrix}]$

r = rlevinson(a,efinal) では、ベクトル a が与えられ、r について上記の線形方程式系が解かれます。ここで、a = [1 a(2) … a(p + 1)] です。線形予測のアプリケーションにおいて、r は、予測誤差フィルターへの入力の自己相関列を表します。ここで、r(1) はゼロラグの要素です。以下の図は、このタイプの典型的なフィルターを示しています。ここで、H(z) は最適線形予測子、x(n) は入力信号、 $\hat{x} (n)$ は予測信号、e(n) は予測誤差です。

入力ベクトル a は、z に関して降べきの順に並べられた、この予測誤差フィルターの多項式係数を表します。

$A (z) = 1 + a (2) z^{- 1} + \dots + a (n + 1) z^{- p}$

有効な自己相関列を生成するには、フィルターは最小位相でなければなりません。efinal は、スカラーの予測誤差パワーで、予測誤差信号の分散 σ²(e) と等しくなります。

[r,u] = rlevinson(a,efinal) では、UDU^* 分解から上三角行列 U が返されます。

$R^{- 1} = U E^{- 1} U^{*}$

ここで

$R = [\begin{matrix} r (1) & r {(2)}^{*} & \dots & r {(p)}^{*} \\ r (2) & r (1) & \dots & r {(p - 1)}^{*} \\ ⋮ & ⋱ & ⋱ & ⋮ \\ r (p) & \dots & r (2) & r (1) \end{matrix}]$

で、E は、出力 e に返される対角行列要素です (以下を参照)。この分解により、自己相関行列の逆行列 R⁻¹ を効率的に評価できるようになります。

出力行列 u は、逆レビンソン・ダービン再帰法の各反復からの予測多項式 a を含みます。

$U = [\begin{matrix} a_{1} {(1)}^{*} & a_{2} {(2)}^{*} & \dots & a_{p + 1} {(p + 1)}^{*} \\ 0 & a_{2} {(1)}^{*} & ⋱ & a_{p + 1} {(p)}^{*} \\ 0 & 0 & ⋱ & a_{p + 1} {(p - 1)}^{*} \\ ⋮ & ⋱ & ⋱ & ⋮ \\ 0 & \dots & 0 & a_{p + 1} {(1)}^{*} \end{matrix}]$

ここで、a_i(j) は、i 次の予測フィルター多項式 (たとえば、再帰の i ステップ) の j 番目の係数です。たとえば、5 次の予測フィルター多項式は、次のようになります。

a5 = u(5:-1:1,5)'

u(p+1:-1:1,p+1)' は、入力多項式係数ベクトル a であることに注意してください。

[r,u,k] = rlevinson(a,efinal) では、反射係数を含む長さ p + 1 のベクトル k が返されます。反射係数は、u の最初の行に共役なものを設定しています。

k = conj(u(1,2:end))

[r,u,k,e] = rlevinson(a,efinal) では、逆レビンソン・ダービン再帰法の各反復からの予測誤差を含む長さ p + 1 のベクトルが返されます。e(1) は 1 次モデルからの予測誤差、e(2) は 2 次モデルからの予測誤差というようになります。

これらの予測誤差値は、R⁻¹ の UDU^* 分解内の行列 E の対角要素から構成されます。

$R^{- 1} = U E^{- 1} U^{*}$

例

すべて折りたたむ

最適な自己回帰モデルの次数

ライブスクリプトを開く

自己回帰モデルを使用して、ノイズ内の 2 つの正弦波のスペクトルを推定します。逆レビンソン・ダービン再帰法で返されたモデルのグループから最適なモデル次数を選択します。

信号を生成します。1 kHz のサンプルレートと 50 秒の信号持続時間を指定します。正弦波には 50 Hz と 55 Hz の周波数があります。ノイズの分散は 0.2² です。

Fs = 1000;
t = (0:50e3-1)'/Fs;
x = sin(2*pi*50*t) + sin(2*pi*55*t) + 0.2*randn(50e3,1);

自己回帰モデルパラメーターを推定します。

[a,e] = arcov(x,100);
[r,u,k] = rlevinson(a,e);

次数 1、5、25、50 および 100 について、パワースペクトル密度を推定します。

N = [1 5 25 50 100];
nFFT = 8096;
P = zeros(nFFT,5);

for idx = 1:numel(N)
    order = N(idx);
    ARtest = flipud(u(:,order));
    P(:,idx) = 1./abs(fft(ARtest,nFFT)).^2;
end

PSD 推定をプロットします。

F = (0:1/nFFT:1/2-1/nFFT)*Fs;
plot(F, 10*log10(P(1:length(P)/2,:)))
grid

legend([repmat('Order = ',[5 1]) num2str(N')])
xlabel('Frequency (Hz)')
ylabel('dB')
xlim([35 70])

Figure contains an axes object. The axes object with xlabel Frequency (Hz), ylabel dB contains 5 objects of type line. These objects represent Order = 1, Order = 5, Order = 25, Order = 50, Order = 100.