ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

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) の複素共役を表します。

[r(1)r(2)r(p)r(2)r(1)r(p1)r(p)r(2)r(1)][a(2)a(3)a(p+1)]=[r(2)r(3)r(p+1)]

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

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

A(z)=1+a(2)z1++a(n+1)zp

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

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

R1=UE1U

ここで、

R=[r(1)r(2)r(p)r(2)r(1)r(p1)r(p)r(2)r(1)]

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

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

U=[a1(1)a2(2)ap+1(p+1)0a2(1)ap+1(p)00ap+1(p1)00ap+1(1)]

ここで、ai(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−1 の UDU* 分解内の行列 E の対角要素から構成されます。

R1=UE1U

すべて折りたたむ

自己回帰モデルを使用して、ノイズ内の 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])

参考文献

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

拡張機能

参考

| | |

R2006a より前に導入