ドキュメンテーション

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

invfreqs

周波数応答データからの連続時間フィルター パラメーターの同定

構文

[b,a] = invfreqs(h,w,n,m)
[b,a] = invfreqs(h,w,n,m,wt)
[b,a] = invfreqs(h,w,n,m,wt,iter)
[b,a] = invfreqs(h,w,n,m,wt,iter,tol)
[b,a] = invfreqs(h,w,n,m,wt,iter,tol,'trace')
[b,a] = invfreqs(h,w,'complex',n,m,...)

説明

invfreqs は、freqs 逆演算です。与えられた複素周波数応答に対応する連続時間伝達関数が求められます。実験室の解析の見地から invfreqs は、振幅データや位相データを伝達関数に変換するのに役立ちます。

[b,a] = invfreqs(h,w,n,m) では、伝達関数の実数の分子係数ベクトル b および分母係数ベクトル a が返されます。

H(s)=B(s)A(s)=b(1)sn+b(2)sn1++b(n+1)a(1)sm+a(2)sm1++a(m+1)

この複素周波数応答は、ベクトル w に指定した周波数点において、ベクトル h で与えられます。スカラー n および m は、分子多項式および分母多項式の希望の次数を指定します。

h の長さは w の長さと等しくしなければなりません。invfreqs では、-wconj(h) を使用して、実数フィルターに対して適切な周波数領域の対称性が保証されます。

[b,a] = invfreqs(h,w,n,m,wt) では、周波数に対する近似誤差に重みが付けられます。wt は、w と同じ長さの重み係数のベクトルです。

[b,a] = invfreqs(h,w,n,m,wt,iter) および

[b,a] = invfreqs(h,w,n,m,wt,iter,tol) は優れたアルゴリズムを提供し、これは、結果として得られる線形システムの安定性を保証し、数値を使用した反復手法により最も適合するものを探索します。iter パラメーターでは、解が収束するか、反復回数が iter に達するかのいずれかが先に起こったときに、invfreqs に反復を終了するように指示します。invfreqs では、修正した勾配ベクトルのノルムが tol 未満となるときが収束と定義されています。ここで、tol はオプションのパラメーターで、既定の設定では 0.01 に設定されます。すべてが 1 の重みベクトルを取得するには、次のステートメントを使用します。

invfreqs(h,w,n,m,[],iter,tol)

[b,a] = invfreqs(h,w,n,m,wt,iter,tol,'trace') では、反復の経過報告がテキストで表示されます。

[b,a] = invfreqs(h,w,'complex',n,m,...) では、複素フィルターが作成されます。この場合、対称性は強制されず、周波数は –π ~ π の間のラジアン単位で指定されます。

すべて折りたたむ

簡単な伝達関数を周波数応答データに変換し、その後、元のフィルター係数に戻します。

a = [1 2 3 2 1 4];
b = [1 2 3 2 3];

[h,w] = freqs(b,a,64);
[bb,aa] = invfreqs(h,w,4,5)
bb = 1×5

    1.0000    2.0000    3.0000    2.0000    3.0000

aa = 1×6

    1.0000    2.0000    3.0000    2.0000    1.0000    4.0000

bb および aa は、b および a とそれぞれ等価です。しかし、aa には正の実数部をもつ極があるため、システムが不安定になります。bbaa の極を表示します。

zplane(bb,aa)

invfreqs の反復アルゴリズムを使用して、システムに対する安定近似を求めます。

[bbb,aaa] = invfreqs(h,w,4,5,[],30)
bbb = 1×5

    0.6816    2.1015    2.6694    0.9113   -0.1218

aaa = 1×6

    1.0000    3.4676    7.4060    6.2102    2.5413    0.0001

新しい極をプロットして、システムが安定していることを検証します。

zplane(bbb,aaa)

実験室で収集した振幅データと位相データをシミュレートする 2 つのベクトル magphase を生成します。さらに、周波数のベクトル w も生成します。

rng('default')

fs = 1000;
t = 0:1/fs:2;
mag = periodogram(sin(2*pi*100*t)+randn(size(t))/10,[],[],fs);
phase = randn(size(mag))/10;
w = linspace(0,fs/2,length(mag))';

invfreqs を使用して、このデータを連続時間の伝達関数に変換します。結果をプロットします。

[b,a] = invfreqs(mag.*exp(1j*phase),w,2,2,[],4);

freqs(b,a)

ヒント

高周波数を使用して高次モデルを作成するときには、a および b に良い条件の値を取得するために、w に存在する最も高い周波数の 1/2 などの係数で除算して周波数をスケーリングすることが重要です。これは、時間の再スケーリングと対応します。

アルゴリズム

既定の設定では、invfreqs では方程式誤差法を使用して、データに対する最適なモデルが識別されます。これは、以下の式

minb,ak=1nwt(k)|h(k)A(w(k))B(w(k))|2

で、線形方程式系を作成し、MATLAB®\ 演算子で解くことにより、ba を求めるものです。ここでは、A(w(k)) と B(w(k)) は、それぞれ周波数 w(k) における多項式 ab のフーリエ変換です。n は周波数点 (hw の長さ) の数です。このアルゴリズムは、Levi (参考文献 [1]) に基づいています。この文献には、いくつかの変形が提案されていますが、高周波数に対する重み関数 wt の考慮があまりなされていません。

優れた (出力誤差) アルゴリズムでは、初期推定値として最初のアルゴリズムの出力が使用され、反復探索 ([2] 参照) に対して、減衰付きガウス-ニュートン法が使用されます。この方法では、実際の周波数応答点と希望の周波数応答点の間の二乗誤差の重み付き和を最小化するという、直接的な問題を解けます。

minb,ak=1nwt(k)|h(k)B(w(k))A(w(k))|2

参考文献

[1] Levi, E. C. “Complex-Curve Fitting.” IRE Trans. on Automatic Control. Vol. AC-4, 1959, pp. 37–44.

[2] Dennis, J. E., Jr., and R. B. Schnabel. Numerical Methods for Unconstrained Optimization and Nonlinear Equations.Englewood Cliffs, NJ: Prentice-Hall, 1983.

R2006a より前に導入