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

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

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

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

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

パラメトリックモデリング関数の重要な応用例の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:长度(x) 1) = [];% Remove corr. at negative lagsa = levinson(r,4)

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

r = xcorr(x,'biased'); r(1:length(x)-1) = [];% Remove corr. at negative lagsa = 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は、このフィルターのインパスル応答がどの程度、元のシーケンスと一致するかを示します。

formatlong[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 iterationsa = 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 を求めます。

ここでは、線形方程式系を作成し、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は,実際の周波数応答点と希望の応答との間の二乗誤差に重みを付けたものの和を最小にする、直接的な問題を解く優れた ("出力誤差") アルゴリズムを使用しています。

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

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

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

FVTool (振幅と位相応答を選択) を使用して、最初のアルゴリズムと 2 番目のアルゴリズムの結果を元のバタワースフィルターと比較します。

fvtool(b,a,b4,a4,b30,a30)

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

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