パラメトリックモデリング

パラメトリックモデリングについて

パラメトリックモデリング手法では、信号、システム、または変動過程を記述する数学モデルに関するパラメーターが求められます。これらの手法では、システムに関する既知の情報を使用してモデルが決定されます。パラメトリックモデリングの応用例としては、音声や音楽の合成、データ圧縮、高分解スペクトル推定、通信、製造、シミュレーションなどがあります。

使用可能なパラメトリックモデリング関数

このツールボックスのパラメトリックモデリング関数では、有理伝達関数モデルが使用されます。未知のシステムに関する適切な情報 (インパルス、または、周波数応答データ、あるいは、入力シーケンスや出力シーケンス) が与えられると、これらの関数によって、システムをモデリングする線形システムの係数が求められます。

パラメトリックモデリング関数の重要な応用例の1 つは、時間応答、または、周波数応答からフィルターを設計するものです。

ここでは、このツールボックスのパラメトリックモデリング関数について概要を説明します。

領域	関数	説明
時間	`arburg`	レビンソン・ダービンアルゴリズムを使用して、入力データシーケンスをモデル化する全極フィルター係数を作成
	`arcov`	前方予測誤差を最小にすることにより、入力データシーケンスをモデル化する全極フィルター係数を作成
	`armcov`	前方/後方予測誤差を最小にすることにより、入力データシーケンスをモデル化する全極フィルター係数を作成
	`aryule`	自己相関関数の推定を使用して、入力データシーケンスをモデル化する全極フィルター係数を作成
	`lpc`, `levinson`	線形予測符号化。全極再帰型フィルターのインパルス応答が設定したシーケンスに一致するようにフィルターを作成
	`prony`	設定したシーケンスに一致するような IIR フィルターのインパルス応答を作成
	`stmcb`	入力を与えて、設定した出力シーケンスに一致するような IIR フィルター出力を作成
周波数	`invfreqz`, `invfreqs`	複素周波数応答データを与えられ、デジタル、または、アナログフィルターの係数を生成

時間領域ベースのモデリング

関数 lpc、prony、および stmcb では、与えられた時間領域インパルス応答に近似するデジタル有理伝達関数の係数が求められます。アルゴリズムは、結果のモデルの複雑度や精度により異なります。

線形予測

線形予測モデリングでは、信号の各出力サンプル x(k) が、過去の n 個の出力の線形結合である (すなわち、これらの出力から線形的な予測が可能である) ということと、係数がサンプル間で一定であるということを仮定しています。

$x (k) = - a (2) x (k - 1) - a (3) x (k - 2) - \dots - a (n + 1) x (k - n) .$

信号 x の n 次の全極モデルは、以下のステートメントにより得られます。

a = lpc(x,n)

lpc を説明するために、加法性ホワイトノイズをもつ全極フィルターのインパルス応答であるサンプル信号を作成します。

x = impz(1,[1 0.1 0.1 0.1 0.1],10) + randn(10,1)/10;

システムをモデリングする 4 次の全極フィルターの係数は、次のようになります。

a = lpc(x,4)

lpc では、まず xcorr を使用して x の相関関数のバイアス付きの推定値が求められ、次に関数 levinson で実装されているレビンソン・ダービン再帰法を使用してモデル係数 a が求められます。レビンソン・ダービン再帰法は、対称テプリッツ線形方程式を解くための高速のアルゴリズムです。n = 4 に対する lpc のアルゴリズムは以下のようになります。

r = xcorr(x);
r(1:length(x)-1) = [];      % Remove corr. at negative lags
a = levinson(r,4)

以下のように、バイアス付きの相関推定値などの様々な相関推定値を levinson に渡すことにより、ほかの仮定で線形予測係数を作成できます。

r = xcorr(x,'biased');
r(1:length(x)-1) = [];      % Remove corr. at negative lags
a = levinson(r,4)

Prony 法 (ARMA モデリング)

関数 prony では、設定された数の極および零点を使用して信号がモデリングされます。シーケンス x と、分子および分母の次数としてそれぞれ n および m が与えられると、ステートメント

[b,a] = prony(x,n,m)

によって、インパルス応答がシーケンス x を近似する IIR フィルターの分子および分母の係数が求められます。

関数 prony では、Parks and Burrus による[4]に記述されている方法が実装されます。この方法では、AR モデリングの共分散法を修正したものを使用して分母の係数 a を求めた後、生じるフィルターのインパルス応答が x の最初の n + 1 個のサンプルに厳密に一致するように分子の係数 b が求められます。このフィルターは必ずしも安定ではありませんが、データシーケンスが、正しい次数の自己回帰移動平均 (ARMA) 過程の場合、係数を正確に回復できる可能性があります。

メモ

関数 prony および stmcb (次節で説明) は、システム同定用語で言えば、ARX モデルと言ったほうが正確だと言えるでしょう。ARMA モデリングでは、入力時のノイズのみが仮定されていますが、ARX では外部入力が仮定されます。prony と stmcb の入力信号は既知のものです。prony の場合はインパルスで、stmcb の場合には任意のデータ列です。

3 次の IIR フィルターを使用したテストシーケンス x (先の lpc の例) のモデルは、以下のようになります。

[b,a] = prony(x,3,3)

コマンド impz は、このフィルターのインパスル応答がどの程度、元のシーケンスと一致するかを示します。

format long
[x impz(b,a,10)]

最初の 4 つのサンプルは正確に一致していることに注目してください。正確に回復する例として、このインパルス応答からバタワースフィルターの係数を回復してみましょう。

[b,a] = butter(4,.2);
h = impz(b,a,26);
[bb,aa] = prony(h,4,4);

この例を実行してみてください。bb および aa は、10^-13 の許容誤差で元のフィルター係数と一致します。

スティグリッツ・マクブライド法 (ARMA モデリング)

関数 stmcb は、所定の近似インパルス応答 x および任意の数の零点と極について、システム b(z)/a(z) の係数を求めます。この関数は、システムの挙動を表わす入出力シーケンス、あるいは単にシステムのインパルス応答のどちらかを使用して、未知のシステムを同定します。既定モードでは、stmcb は、prony と同様に働きます。

[b,a] = stmcb(x,3,3)

stmcb は、また、以下のように与えられた入力シーケンスや出力シーケンスと一致するシステムも求めます。

y = filter(1,[1 1],x);     % Create an output signal.
[b,a] = stmcb(y,x,0,1)

この例では、stmcb は x から y の作成に使用されたシステムを正確に同定しています。

スティグリッツ・マクブライド法は、フィルター出力と与えられた出力信号との間の信号の誤差を最小限に抑えようとしながら、分子および分母の係数を同時に解く高速反復アルゴリズムです。このアルゴリズムは、通常急速に収束しますが、モデルの次数が大きすぎる場合には収束しない可能性があります。prony の場合、同様に、stmcb で得られるフィルターは必ずしも安定しません。これは、その厳密なモデリング手法のためです。

stmcb は、いくつかの重要なアルゴリズムを制御するパラメーターの設定を行います。データのモデリングに問題がある場合には、これらのパラメーターを変更してください。既定の 5 回の反復回数を変更し、分母の係数の初期推定値を設定するには、次のようにします。

n = 10;             % Number of iterations
a = lpc(x,3);       % Initial estimates for denominator
[b,a] = stmcb(x,3,3,n,a);

この関数は、ユーザーが初期値を設定しない場合には、prony で作成した全極モデルを初期設定値として使用します。

関数 lpc、prony、および stmcb を比較するために、以下のようにそれぞれの場合の誤差が計算されます。

a1 = lpc(x,3);
[b2,a2] = prony(x,3,3);
[b3,a3] = stmcb(x,3,3);
[x-impz(1,a1,10)  x-impz(b2,a2,10)  x-impz(b3,a3,10)]

IIR モデルの次数を与えて、モデリング機能を比較すると、上の結果から、この例の場合には stmcb が最良で、続いて、prony、lpc の順となることがわかります。この相対的な性能はモデリング関数の特色を示しています。

周波数領域ベースのモデリング

関数 invfreqs および invfreqz では、freqs および freqz の逆演算が実行され、与えられた複素周波数応答に一致する指定の次数のアナログまたはデジタル伝達関数が求められます。以下の例では、invfreqz を示していますが、この考察は invfreqs にも適用されます。

単純なデジタルフィルターの周波数応答から元のフィルター係数を求めるには、次のようにします。

[b,a] = butter(4,0.4)         % Design Butterworth lowpass
[h,w] = freqz(b,a,64);        % Compute frequency response
[b4,a4] = invfreqz(h,w,4,4)   % Model: n = 4, m = 4

周波数ベクトル w はラジアン/サンプル単位ですが、この周波数は、等間隔である必要はありません。invfreqz は、この周波数データに適合する任意の次数のフィルターが求められます。3 次のフィルターの場合には、以下のようになります。

[b4,a4] = invfreqz(h,w,3,3)   % Find third-order IIR

invfreqs と invfreqz は、共に実数係数のフィルターが設計されます。正の周波数 f のデータ点に対して、これらの関数は周波数応答を f と -f の両方で適合させます。

既定では、invfreqz は方程式誤差法を使用して、データから最良のモデルを同定します。これは、以下の式で b と a を求めます。

$\min_{b, a} \sum_{k = 1}^{n} w t (k) {| h (k) A (w (k)) - B (w (k)) |}^{2}$

ここでは、連立線形方程式を作成し、MATLAB^® の \ 演算子で解くことにより、b と a を求めます。ここで、A(w(k)) と B(w(k)) は、それぞれ周波数 w(k) における多項式 a と b のフーリエ変換です。n は周波数点 (h と w の長さ) の数です。wt(k) は、異なる周波数での誤差に対してその誤差に重みを付けます。次の構文

invfreqz(h,w,n,m,wt)

は、重みベクトルを含みます。このモードでは、invfreqz から得られるフィルターが安定であるという保証はありません。

invfreqz は、実際の周波数応答点と希望の応答との間の二乗誤差に重みを付けたものの和を最小にする、直接的な問題を解く優れた ("出力誤差") アルゴリズムを使用しています。

$\min_{b, a} \sum_{k = 1}^{n} w t (k) {| h (k) - \frac{B (w (k))}{A (w (k))} |}^{2}$

このアルゴリズムを使用するには、以下のように重みベクトルパラメーターの後に反復回数を入力引数として設定します。

wt = ones(size(w));   % Create unit weighting vector
[b30,a30] = invfreqz(h,w,3,3,wt,30)  % 30 iterations

得られるフィルターは常に安定しています。

数値的な適合の優位性を検証するには、以下のようにします。

sum(abs(h-freqz(b4,a4,w)).^2)    % Total error, algorithm 1
sum(abs(h-freqz(b30,a30,w)).^2)  % Total error, algorithm 2

パラメトリック モデリング