levinson

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

構文

a = levinson(r,n)

[a,e,k] = levinson(___)

説明

a = levinson(r,n) は、自己相関列として r をもつ n 次の自己回帰線形過程の係数を返します。

また、[a,e,k] = levinson(___) は予測誤差 e および反射係数 k を返します。

例

自己回帰過程の係数

次の式で表される自己回帰過程の係数を推定します。

$$x(n)=0.1x(n-1)-0.8x(n-2)-0.27x(n-3)+w(n).$$

a = [1 0.1 -0.8 -0.27];

分散 0.4 のホワイト ノイズをフィルター処理することによって過程の実現を生成します。

v = 0.4;
w = sqrt(v)*randn(15000,1);
x = filter(1,a,w);

相関関数を推定します。負のラグの相関値を破棄します。レビンソン・ダービン再帰法を使用して、モデル係数を推定します。予測誤差が入力の分散に対応していることを確認します。

[r,lg] = xcorr(x,'biased');
r(lg<0) = [];

[ar,e] = levinson(r,numel(a)-1)

ar = 1×4

    1.0000    0.0772   -0.7954   -0.2493

e = 
0.3909

16 次モデルの反射係数を推定します。95% の信頼限界の外側に位置している反射係数のみが正しいモデル次数に対応していることを確認します。詳細については、偏自己相関列による AR 次数選択を参照してください。

[~,~,k] = levinson(r,16);
stem(k,'filled')

conf = sqrt(2)*erfinv(0.95)/sqrt(15000);
hold on
[X,Y] = ndgrid(xlim,conf*[-1 1]);
plot(X,Y,'--r')
hold off

複数の実現の予測誤差

次の式で表される自己回帰過程の係数を生成します。

$$x(n)=0.1x(n-1)-0.8x(n-2)-0.27x(n-3)+w(n).$$

a = [1 0.1 -0.8 -0.27];

分散が異なるホワイト ノイズをフィルター処理することにより、過程の実現を 5 件作成します。

nr = 5;
v = rand(1,nr)

v = 1×5

    0.8147    0.9058    0.1270    0.9134    0.6324

w = sqrt(v).*randn(15000,nr);
x = filter(1,a,w);

相関関数を推定します。相互相関の項と負のラグの相関値を破棄します。レビンソン・ダービン再帰法を使用して、正しいモデル次数の予測誤差を推定し、その予測誤差が入力ノイズ信号の分散に対応することを確認します。

[r,lg] = xcorr(x,'biased');

[~,e] = levinson(r(lg>=0,1:nr+1:end),numel(a)-1)

e = 5×1

    0.7957
    0.9045
    0.1255
    0.9290
    0.6291

入力引数

r — 自己相関列

ベクトル | 行列

自己相関列。ベクトルまたは行列として指定します。r が行列の場合、関数は r の各列の係数を求め、それらを a の行に返します。

例: [r,lg] = xcorr(randn(1000,1),'biased'); r(lg<0) = [] では、正のラグに関する 1000 サンプルのランダム信号の自己相関列を推定します。

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

n — モデル次数

length(r) – 1 (既定値) | 正の整数スカラー

モデル次数。正の整数スカラーとして指定します。

データ型: single | double

出力引数

a — 自己回帰線形過程の係数

行ベクトル | 行列

自己回帰線形過程の係数。行ベクトルまたは行列として返されます。フィルター係数は、以下のように z^-1 の降べきの順に並べられます。

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

r が行列の場合、a の各行は r の列に対応します。

e — 予測誤差

スカラー | 列ベクトル

予測誤差。スカラーまたは列ベクトルとして返されます。r が行列の場合、e の各要素は r の列に対応します。

k — 反射係数

列ベクトル | 行列

反射係数。長さ n の列ベクトルとして返されます。r が行列の場合、k の各列は r の列に対応します。

メモ

k は、a 係数の計算中に内部的に計算されます。このため、k を同時に返す方が、tf2latc を使用して a を k に変換するよりも効率的です。

アルゴリズム

レビンソン・ダービン再帰法は、設定された確定的な自己相関列をもつ全極 IIR フィルターを求めるためのアルゴリズムです。この再帰法は、フィルター設計、コード化、スペクトル推定などで使用されています。levinson で生成されるフィルターは、最小位相となります。

levinson では、対称テプリッツ連立線形方程式を解けます。

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

ここで、r = [r(1) ... r(n + 1)] は入力自己相関ベクトルであり、r(i)^* は r(i) の複素共役を表します。入力 r は、自己相関係数のベクトルで、ラグ 0 は最初の要素 r(1) です。

メモ

r が有効な自己相関列でない場合、解が存在する場合であっても、関数 levinson は NaN を返すことがあります。

このアルゴリズムは O(n²) フロップスを必要とし、n が大きい場合は、MATLAB^® のバックスラッシュ コマンドよりはるかに効率的です。ただし、関数 levinson ではできるだけ高速に演算を行うために、低次の場合には \ を使用しています。

参照

[1] Ljung, Lennart. System Identification: Theory for the User. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

拡張機能

C/C++ コード生成

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

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

• この関数のコード生成には DSP System Toolbox™ ソフトウェアが必要です。

• 指定する場合、入力 n はスカラーでなければなりません。

バージョン履歴

R2006a より前に導入

参考

lpc | prony | rlevinson | schurrc | stmcb