rlevinson

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

構文

r = rlevinson(a,eFinal)

[r,u] = rlevinson(___)

[___,k] = rlevinson(___)

[___,eEvolution] = rlevinson(___)

説明

r = rlevinson(a,eFinal) は、a で指定された係数と最終予測誤差パワー e をもつ予測フィルターの自己相関列 r を返します。

[r,u] = rlevinson(___) は、再帰予測フィルターの多項式行列も返します。

[___,k] = rlevinson(___) は、反射係数のリストも返します。

[___,eEvolution] = rlevinson(___) は、逆レビンソン・ダービン再帰法の各反復処理で得られた予測誤差も返します。

例

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

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

50,000 サンプルの信号を生成します。1 kHz のサンプル レートと 50 秒の信号持続時間を指定します。正弦波には 50 Hz と 55 Hz の周波数があります。ノイズの分散は ${\left(0.2\right)}^2$ です。

rng default
Ns = 50e3;
Fs = 1000;
t = (0:Ns-1)'/Fs;

x = sin(2*pi*50*t) + sin(2*pi*55*t) + 0.2*randn(Ns,1);

次数 100 のモデルを使用して自己回帰モデルのパラメーターを推定します。

[ar,eFinal] = arcov(x,100);

逆レビンソン・ダービン再帰法を使用してモデルを解き、自己相関列、再帰予測フィルター、反射係数、および予測誤差を推定します。

[r,u,k,eEvol] = rlevinson(ar,eFinal);

予測誤差をプロットします。最終誤差パワーの水平線を追加します。フィルター次数が増加するにつれて、最終誤差パワーに達するまでエラーが減少します。

plot(eEvol)
hold on
yline(eFinal,"--","eFinal")
hold off
ylim([0 eEvol(1)])
xlabel("Model order")
ylabel("Prediction Error")

次数 1、5、25、50、および 100 について、パワー スペクトル密度 (PSD) を推定します。元の自己回帰パラメーターの PSD を推定します。

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

for idx = 1:numel(N)
    order = N(idx);
    ARpred = flipud(u(:,order));
    psdARpred(:,idx) = 1./abs(fft(ARpred,nFFT)).^2;
end
psdAR = 1./abs(fft(ar,nFFT)).^2;

PSD 推定をプロットします。予測フィルターの PSD 推定値は、次数が増加するにつれて、元の自己回帰パラメーターに徐々に近づきます。

F = (0:1/nFFT:1/2-1/nFFT)*Fs;

figure
plot(F,pow2db(psdARpred(1:length(psdARpred)/2,:)))
hold on
plot(F,pow2db(psdAR(1:length(psdAR)/2)),"--k")
hold off
grid

legend(["Order = "+N';"Original AR"])
xlabel("Frequency (Hz)")
ylabel("Power (dB/Hz)")
xlim([35 70])

入力引数

a — 予測フィルター係数

ベクトル

予測フィルター係数。ベクトルとして指定します。

メモ

有効な自己相関列 r を生成するには、a に関連付けられたフィルターが最小位相でなければなりません。

a が最小位相でない場合、解が存在する場合でも、rlevinson 関数は NaN または Inf の値を返すことがあります。

データ型: single | double
複素数のサポート: あり

eFinal — 最終予測誤差パワー

スカラー

最終予測誤差パワー。スカラーとして指定します。

データ型: single | double
複素数のサポート: あり

出力引数

r — 自己相関列

列ベクトル

自己相関列。a と同じ数の要素をもつ列ベクトルとして返されます。

u — 再帰予測フィルターの多項式行列

行列

再帰予測フィルターの多項式行列。a の要素数と同じ数の行と列をもつ正方行列として返されます。

k — 反射係数

列ベクトル

反射係数のリスト。長さ p の列ベクトルとして返されます。ここで、p+1 は a の要素数です。

eEvolution — 再帰予測誤差

行ベクトル

逆レビンソン・ダービン再帰法の各反復処理で得られた再帰予測誤差のリスト。行ベクトルとして返されます。ベクトル eEvolution の長さは p です。ここで、p+1 は a の要素数です。

アルゴリズム

逆レビンソン・ダービン再帰法は、対称テプリッツ連立線形方程式を解くためのステップダウン アルゴリズムを実装します。

$$\left[\begin{array}{ccccc}r(1)& r{(2)}^{\ast }& \cdots & r{(p-1)}^{\ast }& r{(p)}^{\ast }\\ r(2)& r(1)& \cdots & r{(p-2)}^{\ast }& r{(p-1)}^{\ast }\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ r(p-1)& r(p-2)& \cdots & r(1)& r{(2)}^{\ast }\\ r(p)& r(p-1)& \cdots & r(2)& r(1)\end{array}\right]\left[\begin{array}{c}a(2)\\ a(3)\\ ⋮\\ a(p)\\ a(p+1)\end{array}\right]=\left[\begin{array}{c}-r(2)\\ -r(3)\\ ⋮\\ -r(p)\\ -r(p+1)\end{array}\right],$$

ここで、r = [r(1) r(2) … r(p) r(p + 1)] および a = [1 a(2) … a(p + 1)] です。^* という表記は複素共役転置演算を表します。

線形予測アプリケーションの場合:

• ベクトル r は、予測誤差フィルターへの入力の自己相関列を表します。ここで、r(1) はゼロラグの要素です。テプリッツ行列 R は、自己相関行列です。

  $$R=\left[\begin{array}{ccccc}r(1)& r{(2)}^{\ast }& \cdots & r{(p-1)}^{\ast }& r{(p)}^{\ast }\\ r(2)& r(1)& \cdots & r{(p-2)}^{\ast }& r{(p-1)}^{\ast }\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ r(p-1)& r(p-2)& \cdots & r(1)& r{(2)}^{\ast }\\ r(p)& r(p-1)& \cdots & r(2)& r(1)\end{array}\right]$$

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

• スカラー p は、予測フィルターの次数を表します。

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

次の図は、このタイプの典型的なフィルターを示しています。ここで、H(z) は最適線形予測子、x(n) は入力信号、$$\widehat{x}(n)$$ は予測信号、e(n) は予測誤差信号です。予測誤差の分散 var(e) は、スカラーの予測誤差パワーです。

rlevinson 関数は、予測フィルター係数をもつベクトル a が与えられた場合、この連立方程式を解いて自己相関列 r を計算します。有効な自己相関列を生成するには、フィルターが最小位相でなければなりません。スカラー eFinal は、最終予測誤差パワー var(e) を表します。rlevinson 関数が実装するレビンソン・ダービン再帰法は、この値を目標推定誤差として目指します。

出力 u、k、または eEvolution を要求した場合、rlevinson は、自己相関行列の逆行列 R⁻¹ の UDU^* 分解も実行します。この分解により行列 U と E が求められ、自己相関行列の逆行列 R⁻¹ を効率的に求めることができるようになります。

$$R^{-1}=U{E}^{-1}{U}^{\ast }$$

出力引数 u は行列 U を返します。この行列には、逆レビンソン・ダービン再帰法の各反復処理で得られた再帰予測フィルターの多項式 (a₁、a₂、…、a_p+1) が格納されます。

$$U=\left[\begin{array}{cccc}{a}_1{(1)}^{\ast }& a_2{(2)}^{\ast }& \cdots & a_{p+1}{(p+1)}^{\ast }\\ 0& a_2{(1)}^{\ast }& \cdots & a_{p+1}{(p)}^{\ast }\\ 0& 0& \cdots & a_{p+1}{(p-1)}^{\ast }\\ ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & 0& a_{p+1}{(1)}^{\ast }\end{array}\right]$$

各ベクトル a_i には再帰多項式の係数がリストされます。ここで、a_i(j) は、(i-1) 次の予測フィルター多項式 (すなわち再帰処理の i 番目のステップ) における j 番目の係数です。たとえば、4 次予測フィルター多項式の係数は a₅ として取得できます。ここで、a5 = u(5:-1:1,5) であり、このベクトルは多項式 $$A_5(z)=a_5(1)+a_5(2)z^{-1}+\mathrm{...}+a_5(5)z^{-4}$$ を表します。したがって、U の (p+1) 番目の列には入力多項式の係数ベクトル a が格納されます。これは、u(p+1:-1:1,p+1)' を使用して取得できます。

出力引数 k は、U の最初の行の非対角要素の共役である反射係数をもつ列ベクトルを返します。つまり、k は u(1,2:end)' として返されます。

$$k=\left[\begin{array}{c}{a}_2{(2)}^{\ast }\\ a_3{(3)}^{\ast }\\ ⋮\\ a_{p+1}{(p+1)}^{\ast }\end{array}\right]$$

出力引数 eEvolution は、行列 E の最初の要素を除く対角要素から得られる再帰予測誤差が格納された行ベクトルを返します。

$$E=\left[\begin{array}{cccc}{e}_1& 0& \cdots & 0\\ 0& e_2& \ddots & 0\\ ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & 0& e_{p+1}\end{array}\right]$$

$$eEvolution=[\begin{array}{cccc}{e}_2& e_3& \cdots & e_{p+1}\end{array}]$$

各要素 e_i には、(i-1) 次の予測フィルター多項式 (すなわち再帰処理の i 番目のステップ) における再帰予測誤差がリストされます。たとえば、4 次予測フィルター多項式の予測誤差は、e5 = eEvolution(4) として取得できます。これは、多項式 $$A_5(z)=a_5(1)+a_5(2)z^{-1}+\mathrm{...}+a_5(5)z^{-4}$$ を予測フィルターとして使用した場合の予測誤差です。

• E の最初の対角要素 e₁ は、自己相関列 r の最初の要素 r(1) と一致します。

• E の (p+1) 番目の対角要素 e₁ (これは eEvolution の最後の要素でもあります) は、最終予測誤差パワー eFinal と一致します。

参照

[1] Kay, Steven M. Modern Spectral Estimation: Theory and Application. Englewood Cliffs, NJ: Prentice-Hall, 1988.

拡張機能

C/C++ コード生成

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

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

既定の次元の選択に関する MATLAB との非互換性 (MATLAB Coder)を参照してください。

バージョン履歴

R2006a より前に導入

参考

levinson | lpc | prony | schurrc | stmcb