tfestimate

伝達関数の推定

ページ内をすべて折りたたむ

構文

txy = tfestimate(x,y)

txy = tfestimate(x,y,window)

txy = tfestimate(x,y,window,noverlap)

txy = tfestimate(x,y,window,noverlap,nfft)

txy = tfestimate(___,'mimo')

[txy,w] = tfestimate(___)

[txy,f] = tfestimate(___,fs)

[txy,w] = tfestimate(x,y,window,noverlap,w)

[txy,f] = tfestimate(x,y,window,noverlap,f,fs)

[___] = tfestimate(x,y,___,freqrange)

[___] = tfestimate(___,'Estimator',est)

tfestimate(___)

説明

txy = tfestimate(x,y) では、一連の周波数で評価された入力信号 x と出力信号 y の間で伝達関数の推定を求めます。

x と y の両方がベクトルの場合、それらは同じ長さでなければなりません。
信号の 1 つが行列でもう 1 つがベクトルの場合は、ベクトルの長さが行列の行数と等しくなければなりません。関数はベクトルを展開し、列ごとの伝達関数推定の行列を返します。
x と y が同じ行数だが異なる列数をもつ行列の場合、txy は多入力/多出力 (MIMO) 伝達関数になり、すべての入力信号と出力信号を結合します。txy は 3 次元配列です。x が m 列をもち、y が n 列をもつ場合、txy は n 列および m ページをもちます。詳細については、伝達関数を参照してください。
x と y が等しいサイズの行列の場合、tfestimate は列方向に動作します (txy(:,n) = tfestimate(x(:,n),y(:,n)))。MIMO 推定を取得するには、引数リストに 'mimo' を追加します。

例

txy = tfestimate(x,y,window) は、window を使用して x と y をセグメントに分割し、ウィンドウ処理を実行します。

txy = tfestimate(x,y,window,noverlap) は、隣り合ったセグメント間で noverlap 個のサンプルのオーバーラップを使用します。

txy = tfestimate(x,y,window,noverlap,nfft) は、nfft サンプリング点を使用して離散フーリエ変換を計算します。

txy = tfestimate(___,'mimo') は、行列入力の MIMO 伝達関数を計算します。この構文には、前の構文の入力引数を任意に組み合わせて含めることができます。

[txy,w] = tfestimate(___) は、正規化周波数のベクトル w を返し、これにより伝達関数が推定されます。

例

[txy,f] = tfestimate(___,fs) は、サンプルレート fs で表される周波数のベクトル f を返し、これにより、伝達関数が推定されます。fs は tfestimate に対する 6 番目の数値入力でなければなりません。サンプルレートを入力した場合でも前のオプション引数の既定値を使用するには、これらの引数を空 [] として指定します。

[txy,w] = tfestimate(x,y,window,noverlap,w) は、w で指定される正規化周波数における伝達関数推定を返します。

[txy,f] = tfestimate(x,y,window,noverlap,f,fs) は、f で指定される周波数における伝達関数推定を返します。

[___] = tfestimate(x,y,___,freqrange) では、freqrange で指定される周波数範囲にわたる伝達関数推定が返されます。freqrange の有効なオプションは、'onesided'、'twosided' および 'centered' です。

例

[___] = tfestimate(___,'Estimator',est) は、推定器 est を使用して伝達関数を推定します。est の有効なオプションは、'H1' および 'H2' です。

出力引数を設定せずに tfestimate(___) を使用すると、現在の Figure ウィンドウに伝達関数推定がプロットされます。

例

すべて折りたたむ

2 つのシーケンス間の伝達関数

ライブスクリプトを開く

2 つのシーケンス x および y 間の伝達関数推定を計算して、プロットします。シーケンス x はホワイトガウスノイズから構成されます。y は、x を正規化されたカットオフ周波数 $0.2 π$ ラジアン/サンプルをもつ 30 次ローパスフィルターでフィルター処理した結果です。フィルターの設計には箱型ウィンドウを使用します。伝達関数の推定のために、500 Hz のサンプルレートと長さ 1024 のハミングウィンドウを指定します。

h = fir1(30,0.2,rectwin(31));
x = randn(16384,1);
y = filter(h,1,x);

fs = 500;
tfestimate(x,y,1024,[],[],fs)

Figure contains an axes object. The axes object with title Transfer Function Estimate via Welch contains an object of type line.

fvtool を使用して、伝達関数がフィルターの周波数応答を近似していることを確認します。

fvtool(h,1,'Fs',fs)

{"String":"Figure Figure 1: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB) contains an object of type line.","Tex":"Magnitude Response (dB)","LaTex":[]}

伝達関数の推定を変数に返し、絶対値をデシベル単位でプロットして同じ結果を得ます。

[Txy,f] = tfestimate(x,y,1024,[],[],fs);

plot(f,mag2db(abs(Txy)))

Figure contains an axes object. The axes object contains an object of type line.

SISO 伝達関数

ライブスクリプトを開く

シンプルな単入力/単出力システムの伝達関数を推定し、定義と比較します。

単位弾性定数のバネで壁につながれた単位質量 $m$ で構成される 1 次元の離散時間振動システムについて考えます。センサーによりこの質量の加速度 $a$ を $F_{s} = 1$ Hz でサンプリングします。ダンパーは質量の運動を妨げるために、速度に比例する力を適用します。減衰定数は $b = 0.01$ です。

2000 個の時間サンプルを生成します。サンプリング間隔は $Δ t = 1 / F_{s}$ と定義します。

Fs = 1;
dt = 1/Fs;
N = 2000;
t = dt*(0:N-1);
b = 0.01;

このシステムは次の状態空間モデルで表すことができます。

$\begin{array}{c} x (k + 1) = A x (k) + B u (k), \\ y (k) = C x (k) + D u (k), \end{array}$

ここで、 $x = {[\begin{matrix} r & v \end{matrix}]}^{T}$ は状態ベクトル、 $r$ および $v$ はそれぞれ質量の位置と速度、 $u$ は駆動力、 $y = a$ は測定出力です。状態空間行列は次のようになります。

$A = \exp (A_{c} Δ t), B = A_{c}^{- 1} (A - I) B_{c}, C = [\begin{array}{c} - 1 & - b \end{array}], D = 1,$

$I$ は $2 \times 2$ の単位行列で、連続時間状態空間行列は次のようになります。

$A_{c} = [\begin{array}{c} 0 & 1 \\ - 1 & - b \end{array}], B_{c} = [\begin{array}{c} 0 \\ 1 \end{array}] .$

Ac = [0 1;-1 -b];
A = expm(Ac*dt);

Bc = [0;1];
B = Ac\(A-eye(size(A)))*Bc;

C = [-1 -b];
D = 1;

測定間隔の 1/2 で、質量がランダム入力により駆動されます。状態空間モデルを使用して、すべてゼロの初期状態からの系の時間発展を計算します。質量の加速度を時間の関数としてプロットします。

rng default

u = zeros(1,N);
u(1:N/2) = randn(1,N/2);

y = 0;
x = [0;0];
for k = 1:N
    y(k) = C*x + D*u(k);
    x = A*x + B*u(k);
end

plot(t,y)

Figure contains an axes object. The axes object contains an object of type line.

システムの伝達関数を周波数の関数として推定します。DFT 点 2048 個を使用し、形状係数 15 のカイザーウィンドウを指定します。隣接するセグメント間のオーバーラップに既定値を使用します。

nfs = 2048;
wind = kaiser(N,15);

[txy,ft] = tfestimate(u,y,wind,[],nfs,Fs);

離散時間システムの周波数応答関数は、単位円で評価されたシステムの時間領域伝達関数の Z 変換として表すことができます。tfestimate で計算される推定が、この定義と一致することを確認します。

[b,a] = ss2tf(A,B,C,D);

fz = 0:1/nfs:1/2-1/nfs;
z = exp(2j*pi*fz);
frf = polyval(b,z)./polyval(a,z);

plot(ft,20*log10(abs(txy)))
hold on
plot(fz,20*log10(abs(frf)))
hold off
grid
ylim([-60 40])

Figure contains an axes object. The axes object contains 2 objects of type line.

tfestimate の組み込み機能を使用して推定をプロットします。

tfestimate(u,y,wind,[],nfs,Fs)

Figure contains an axes object. The axes object with title Transfer Function Estimate via Welch contains an object of type line.

MIMO システムの伝達関数

ライブスクリプトを開く

シンプルな多入力/多出力システムの伝達関数を推定します。

2 つの壁の間に拘束された 2 つの質量 $m_{1}$ および $m_{2}$ で構成される理想的な 1 次元の振動システムについて考えます。物体は $m_{1} = 1$ および $m_{2} = μ$ を満たします。各質量は弾性定数 $k$ をもつバネでそれぞれ近い方の壁につながれています。同じバネで 2 つの質量がつながれています。3 つのダンパーは、速度に比例する力を適用することによって質量の運動を妨げます。減衰定数は $b$ です。センサーによりこれらの質量の加速度 $a_{1}$ および $a_{2}$ を $F_{s} = 50$ Hz でサンプリングします。

30000 回のサンプルを生成しますが、これは 600 秒に相当します。サンプリング間隔は $Δ t = 1 / F_{s}$ と定義します。

Fs = 50;
dt = 1/Fs;
N = 30000;
t = dt*(0:N-1);

このシステムは次の状態空間モデルで表すことができます。

$\begin{array}{c} x (k + 1) = A x (k) + B u (k), \\ y (k) = C x (k) + D u (k), \end{array}$

ここで、 $x = {[\begin{array}{c} r_{1} & v_{1} & r_{2} & v_{2} \end{array}]}^{T}$ は状態ベクトル、 $r_{i}$ および $v_{i}$ はそれぞれ $i$ 番目の質量の位置と速度、 $u = {[\begin{array}{c} u_{1} & u_{2} \end{array}]}^{T}$ は入力駆動力のベクトル、および $y = {[\begin{array}{c} a_{1} & a_{2} \end{array}]}^{T}$ は出力ベクトルです。状態空間行列は次のようになります。

$A = \exp (A_{c} Δ t), B = A_{c}^{- 1} (A - I) B_{c}, C = [\begin{array}{c} - 2 k & - 2 b & k & b \\ k / μ & b / μ & - 2 k / μ & - 2 b / μ \end{array}], D = [\begin{array}{c} 1 & 0 \\ 0 & 1 / μ \end{array}],$

$I$ は $4 \times 4$ の単位行列で、連続時間状態空間行列は次のようになります。

$A_{c} = [\begin{array}{c} 0 & 1 & 0 & 0 \\ - 2 k & - 2 b & k & b \\ 0 & 0 & 0 & 1 \\ k / μ & b / μ & - 2 k / μ & - 2 b / μ \end{array}], B_{c} = [\begin{array}{c} 0 & 0 \\ 1 & 0 \\ 0 & 0 \\ 0 & 1 / μ \end{array}] .$

$k = 400$ 、 $b = 0$ 、および $μ = 1 / 10$ を設定します。

k = 400;
b = 0;
m = 1/10;

Ac = [0 1 0 0;-2*k -2*b k b;0 0 0 1;k/m b/m -2*k/m -2*b/m];
A = expm(Ac*dt);
Bc = [0 0;1 0;0 0;0 1/m];
B = Ac\(A-eye(4))*Bc;
C = [-2*k -2*b k b;k/m b/m -2*k/m -2*b/m];
D = [1 0;0 1/m];

測定全体をとおして質量はランダム入力により駆動されます。状態空間モデルを使用して、すべてゼロの初期状態からの系の時間発展を計算します。

rng default
u = randn(2,N);

x = [0;0;0;0];
for kk = 1:N
    y(:,kk) = C*x + D*u(:,kk);
    x = A*x + B*u(:,kk);
end

入力および出力データを使用して、システムの伝達関数を周波数の関数として推定します。'mimo' オプションを指定して、4 つすべての伝達関数を生成します。5000 サンプルのハンウィンドウを使用し、信号をセグメントに分割します。隣り合ったセグメント間のオーバーラップを 2500 サンプル、DFT 点を $2^{14}$ に指定します。推定をプロットします。

wind = hann(5000);
nov = 2500;

[q,fq] = tfestimate(u',y',wind,nov,2^14,Fs,'mimo');

理論上の伝達関数を、単位円で評価された時間領域伝達関数の Z 変換として計算します。

nfs = 2^14;

fz = 0:1/nfs:1/2-1/nfs;
z = exp(2j*pi*fz);

[b1,a1] = ss2tf(A,B,C,D,1);
[b2,a2] = ss2tf(A,B,C,D,2);

frf(1,:,1) = polyval(b1(1,:),z)./polyval(a1,z);
frf(1,:,2) = polyval(b1(2,:),z)./polyval(a1,z);

frf(2,:,1) = polyval(b2(1,:),z)./polyval(a2,z);
frf(2,:,2) = polyval(b2(2,:),z)./polyval(a2,z);

理論上の伝達関数とそれらに対応する推定をプロットします。

for jk = 1:2
    for kj = 1:2
        subplot(2,2,2*(jk-1)+kj)
        plot(fq,20*log10(abs(q(:,jk,kj))))
        hold on
        plot(fz*Fs,20*log10(abs(frf(jk,:,kj))))
        hold off
        grid
        title(['Input ' int2str(kj) ', Output ' int2str(jk)])
        axis([0 Fs/2 -50 100])
    end
end

Figure contains 4 axes objects. Axes object 1 with title Input 1, Output 1 contains 2 objects of type line. Axes object 2 with title Input 2, Output 1 contains 2 objects of type line. Axes object 3 with title Input 1, Output 2 contains 2 objects of type line. Axes object 4 with title Input 2, Output 2 contains 2 objects of type line.

伝達関数は期待値 $ω_{1, 2} / 2 π$ において最大値を持ちます。ここで、 $ω$ はモードの行列の固有値です。

sqrt(eig(k*[2 -1;-1/m 2/m]))/(2*pi)

ans = 2×1

    3.8470
   14.4259

$b = 0.1$ を設定することでシステムに減衰を追加します。減衰システムの時間発展を同じ駆動力で計算します。MIMO 伝達関数の $H_{2}$ の推定を同じウィンドウとオーバーラップを使用して計算します。tfestimate 機能を使用して推定をプロットします。

b = 0.1;

Ac = [0 1 0 0;-2*k -2*b k b;0 0 0 1;k/m b/m -2*k/m -2*b/m];
A = expm(Ac*dt);
B = Ac\(A-eye(4))*Bc;
C = [-2*k -2*b k b;k/m b/m -2*k/m -2*b/m];

x = [0;0;0;0];
for kk = 1:N
    y(:,kk) = C*x + D*u(:,kk);
    x = A*x + B*u(:,kk);
end

clf
tfestimate(u',y',wind,nov,[],Fs,'mimo','Estimator','H2')
legend('I1, O1','I1, O2','I2, O1','I2, O2')

Figure contains an axes object. The axes object with title Transfer Function Estimate via Welch contains 4 objects of type line. These objects represent I1, O1, I1, O2, I2, O1, I2, O2.

yl = ylim;

推定を理論予測と比較します。

[b1,a1] = ss2tf(A,B,C,D,1);
[b2,a2] = ss2tf(A,B,C,D,2);

frf(1,:,1) = polyval(b1(1,:),z)./polyval(a1,z);
frf(1,:,2) = polyval(b1(2,:),z)./polyval(a1,z);

frf(2,:,1) = polyval(b2(1,:),z)./polyval(a2,z);
frf(2,:,2) = polyval(b2(2,:),z)./polyval(a2,z);

plot(fz*Fs,20*log10(abs(reshape(permute(frf,[2 1 3]),[nfs/2 4]))))
legend('I1, O1','I1, O2','I2, O1','I2, O2')
ylim(yl)
grid

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent I1, O1, I1, O2, I2, O1, I2, O2.

入力引数

すべて折りたたむ

`x` — 入力信号
ベクトル | 行列

ベクトルまたは行列として指定される入力信号。

例: cos(pi/4*(0:159))+randn(1,160) は、ホワイトガウスノイズに含まれる正弦波を指定します。

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

`y` — 出力信号
ベクトル | 行列

ベクトルまたは行列として指定される出力信号。

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

`window` — ウィンドウ
整数 | ベクトル | `[]`

ウィンドウ。整数、あるいは行ベクトルまたは列ベクトルとして指定します。window は信号をセグメントに分割するために使用します。

window が整数の場合、tfestimate は x と y を長さ window のセグメントに分割し、各セグメントにその長さのハミングウィンドウを適用します。
window がベクトルの場合、tfestimate は x と y をベクトルと同じ長さのセグメントに分割し、window を使用して各セグメントにウィンドウを適用します。

x と y の長さが noverlap 個のオーバーラップサンプルをもつ整数個のセグメントに厳密に分割できない場合、信号はそれに応じた長さで打ち切られます。

window を空として指定した場合、tfestimate はハミングウィンドウを使用して、x と y を noverlap 個のオーバーラップサンプルをもつ 8 個のセグメントに分割します。

利用可能なウィンドウのリストについては、ウィンドウを参照してください。

例: hann(N+1) と (1-cos(2*pi*(0:N)'/N))/2 は、いずれも長さ N + 1 のハンウィンドウを指定します。

データ型: single | double

`noverlap` — オーバーラップするサンプル数
正の整数 | `[]`

オーバーラップするサンプル数。正の整数で指定します。

window がスカラーの場合、noverlap は window より小さくなければなりません。
window がベクトルの場合、noverlap は window の長さより小さくなければなりません。

noverlap を空として指定した場合、tfestimate はセグメント間で 50% のオーバーラップが発生する数を使用します。

データ型: double | single

`nfft` — DFT 点の数
正の整数 | `[]`

正の整数として指定する DFT 点の数。nfft を空として指定した場合、tfestimate によりこの引数が max(256,2^p) に設定されます。ここで p = ⌈log₂ N⌉ は入力信号の長さ N についての値であり、⌈ ⌉ 記号は天井関数を表します。

データ型: single | double

`fs` — サンプルレート
正のスカラー

サンプルレート。正のスカラーで指定します。サンプルレートは単位時間あたりのサンプル数です。時間の単位が秒の場合、サンプルレートの単位は Hz です。

`w` — 正規化周波数
ベクトル

正規化周波数。少なくとも 2 つの要素をもつ行ベクトルまたは列ベクトルとして指定します。正規化周波数の単位はラジアン/サンプルです。

例: w = [pi/4 pi/2]

データ型: double

`f` — 周波数
ベクトル

周波数。少なくとも 2 つの要素をもつ行ベクトルまたは列ベクトルとして指定します。周波数の単位は単位時間あたりのサイクルです。単位時間はサンプルレート fs で指定されます。fs の単位がサンプル/秒であれば、f の単位は Hz です。

例: fs = 1000; f = [100 200]

データ型: double

`freqrange` — 伝達関数推定のための周波数範囲
`'onesided'` | `'twosided'` | `'centered'`

'onesided'、'twosided' または 'centered' で指定する、伝達関数推定の周波数範囲。既定値は、実数値信号の場合は 'onesided'、複素数値信号の場合は 'twosided' です。

'onesided' — 2 つの実数値の入力信号 x と y の間の伝達関数の片側推定を返します。nfft が偶数の場合、txy は nfft/2 + 1 行で、計算区間は [0,π] ラジアン/サンプルです。nfft が奇数の場合、txy は (nfft + 1)/2 行で、計算区間は [0,π) ラジアン/サンプルです。fs を指定すると、nfft が偶数の場合は、対応する計算区間は [0,fs/2] サイクル/単位時間で、nfft が奇数の場合は、[0,fs/2) サイクル/単位時間です。
'twosided' — 実数値または複素数値の 2 つの入力信号 x と y の間の伝達関数の両側推定を返します。この場合、txy は nfft 行で、計算区間は [0,2π) ラジアン/サンプルです。fs を指定した場合、計算区間は [0,fs) サイクル/単位時間となります。
'centered' — 実数値または複素数値の 2 つの入力信号 x と y の間の伝達関数の中央に揃えた両側推定を返します。この場合、txy は nfft 行で、偶数の nfft については区間 (–π,π] ラジアン/サンプルで、奇数の nfft については (–π,π) ラジアン/サンプルで計算されます。fs を指定すると、nfft が偶数の場合は、対応する計算区間は (–fs/2, fs/2] サイクル/単位時間で、nfft が奇数の場合は、(–fs/2, fs/2) サイクル/単位時間です。

`est` — 伝達関数の推定器
`'H1'` (既定値) | `'H2'`

伝達関数の推定器。'H1' または 'H2' で指定します。

ノイズが入力信号と無相関な場合は、'H1' を使用します。
ノイズが出力信号と無相関な場合は、'H2' を使用します。この場合、入力信号の数が出力信号の数と等しくなければなりません。

詳細については、伝達関数を参照してください。

出力引数

すべて折りたたむ

`txy` — 伝達関数の推定
ベクトル | 行列 | 3 次元配列

伝達関数推定。ベクトル、行列または 3 次元配列として返されます。

`w` — 正規化周波数
ベクトル

正規化周波数。実数値の列ベクトルとして返されます。

`f` — 巡回周波数
ベクトル

巡回周波数。実数値の列ベクトルとして返されます。

詳細

すべて折りたたむ

伝達関数

入力 x と出力 y の関係は、線形の時不変 "伝達関数" txy によってモデリングされます。周波数領域では、Y(f) = H(f)X(f) となります。

単入力/単出力システムの場合、伝達関数の H₁ の推定は、以下で与えられます。

$H_{1} (f) = \frac{P_{y x} (f)}{P_{x x} (f)},$
ここで、P_yx は x と y のクロスパワースペクトル密度、P_xx は x のパワースペクトル密度です。この推定は、ノイズがシステム入力と相関していないと仮定しています。
多入力/多出力 (MIMO) システムの場合、H₁ 推定器は以下のようになります。

$H_{1} (f) = P_{Y X} (f) P_{X X}^{- 1} (f) = [\begin{matrix} P_{y_{1} x_{1}} (f) & P_{y_{1} x_{2}} (f) & \dots & P_{y_{1} x_{m}} (f) \\ P_{y_{2} x_{1}} (f) & P_{y_{2} x_{2}} (f) & \dots & P_{y_{2} x_{m}} (f) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ P_{y_{n} x_{1}} (f) & P_{y_{n} x_{2}} (f) & \dots & P_{y_{n} x_{m}} (f) \end{matrix}] {[\begin{matrix} P_{x_{1} x_{1}} (f) & P_{x_{1} x_{2}} (f) & \dots & P_{x_{1} x_{m}} (f) \\ P_{x_{2} x_{1}} (f) & P_{x_{2} x_{2}} (f) & \dots & P_{x_{2} x_{m}} (f) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ P_{x_{m} x_{1}} (f) & P_{x_{m} x_{2}} (f) & \dots & P_{x_{m} x_{m}} (f) \end{matrix}]}^{- 1}$
これは m 入力と n 出力の場合です。ここで、
- P_{y_ix_k} は k 番目の入力と i 番目の出力のクロスパワースペクトル密度です。
- P_{x_ix_k} は k 番目および i 番目の入力のクロスパワースペクトル密度です。
2 入力と 2 出力の場合、推定器は以下のような行列になります。

$H_{1} (f) = \frac{[\begin{matrix} P_{y_{1} x_{1}} (f) P_{x_{2} x_{2}} (f) - P_{y_{1} x_{2}} (f) P_{x_{2} x_{1}} (f) & P_{y_{1} x_{2}} (f) P_{x_{1} x_{1}} (f) - P_{y_{1} x_{1}} (f) P_{x_{1} x_{2}} (f) \\ P_{y_{2} x_{1}} (f) P_{x_{2} x_{2}} (f) - P_{y_{2} x_{2}} (f) P_{x_{2} x_{1}} (f) & P_{y_{2} x_{2}} (f) P_{x_{1} x_{1}} (f) - P_{y_{2} x_{1}} (f) P_{x_{1} x_{2}} (f) \end{matrix}]}{P_{x_{1} x_{1}} (f) P_{x_{2} x_{2}} (f) - P_{x_{1} x_{2}} (f) P_{x_{2} x_{1}} (f)} .$
単入力/単出力システムの場合、伝達関数の H₂ の推定は以下で与えられます。

$H_{2} (f) = \frac{P_{y y} (f)}{P_{x y} (f)},$
ここで、P_yy は y のパワースペクトル密度、P_xy = P^*_yx は x と y のクロスパワースペクトル密度の複素共役です。この推定は、ノイズがシステム出力と相関していないと仮定しています。
MIMO システムの場合、H₂ 推定器が良好に定義されるのは、入力と出力が同数 (n = m) の場合のみです。推定器は以下のようになります。

$H_{2} (f) = P_{Y Y} (f) P_{X Y}^{- 1} (f) = [\begin{matrix} P_{y_{1} y_{1}} (f) & P_{y_{1} y_{2}} (f) & \dots & P_{y_{1} y_{n}} (f) \\ P_{y_{2} y_{1}} (f) & P_{y_{2} y_{2}} (f) & \dots & P_{y_{2} y_{n}} (f) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ P_{y_{n} y_{1}} (f) & P_{y_{n} y_{2}} (f) & \dots & P_{y_{n} y_{n}} (f) \end{matrix}] {[\begin{matrix} P_{x_{1} y_{1}} (f) & P_{x_{1} y_{2}} (f) & \dots & P_{x_{1} y_{n}} (f) \\ P_{x_{2} y_{1}} (f) & P_{x_{2} y_{2}} (f) & \dots & P_{x_{2} y_{n}} (f) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ P_{x_{n} y_{1}} (f) & P_{x_{n} y_{2}} (f) & \dots & P_{x_{n} y_{n}} (f) \end{matrix}]}^{- 1},$
ここで、
- P_{y_iy_k} は k 番目および i 番目の出力のクロスパワースペクトル密度です。
- P_{x_iy_k} は i 番目の入力と k 番目の出力のクロスパワースペクトル密度の複素共役です。

アルゴリズム

tfestimate では、ウェルチの平均ピリオドグラム法が使用されます。詳細は、pwelch を参照してください。

参照

[1] Vold, Håvard, John Crowley, and G. Thomas Rocklin. “New Ways of Estimating Frequency Response Functions.” Sound and Vibration. Vol. 18, November 1984, pp. 34–38.

拡張機能

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

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

名前と値のペアを使用して指定する引数は、コンパイル時の定数でなければなりません。
可変サイズの window は倍精度でなければなりません。

バージョン履歴

R2006a より前に導入

参考

cpsd | mscohere | periodogram | pwelch

tfestimate

構文

説明

例

2 つのシーケンス間の伝達関数

SISO 伝達関数

MIMO システムの伝達関数

入力引数

x — 入力信号 ベクトル | 行列

y — 出力信号 ベクトル | 行列

window — ウィンドウ 整数 | ベクトル | []

noverlap — オーバーラップするサンプル数 正の整数 | []

nfft — DFT 点の数 正の整数 | []

fs — サンプルレート 正のスカラー

w — 正規化周波数 ベクトル

f — 周波数 ベクトル

freqrange — 伝達関数推定のための周波数範囲 'onesided' | 'twosided' | 'centered'

est — 伝達関数の推定器 'H1' (既定値) | 'H2'

出力引数

txy — 伝達関数の推定 ベクトル | 行列 | 3 次元配列

w — 正規化周波数 ベクトル

f — 巡回周波数 ベクトル

詳細

伝達関数

アルゴリズム

参照

拡張機能

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

バージョン履歴

参考

`x` — 入力信号
ベクトル | 行列

`y` — 出力信号
ベクトル | 行列

`window` — ウィンドウ
整数 | ベクトル | `[]`

`noverlap` — オーバーラップするサンプル数
正の整数 | `[]`

`nfft` — DFT 点の数
正の整数 | `[]`

`fs` — サンプルレート
正のスカラー

`w` — 正規化周波数
ベクトル

`f` — 周波数
ベクトル

`freqrange` — 伝達関数推定のための周波数範囲
`'onesided'` | `'twosided'` | `'centered'`

`est` — 伝達関数の推定器
`'H1'` (既定値) | `'H2'`

`txy` — 伝達関数の推定
ベクトル | 行列 | 3 次元配列

`w` — 正規化周波数
ベクトル

`f` — 巡回周波数
ベクトル

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