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\pi$$ ラジアン/サンプルをもつ 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)

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

freqz(h,1,[],fs)

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

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

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

SISO 伝達関数

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

単位弾性定数のバネで壁につながれた単位質量 $\mathit{m}$ (kg 単位) の物体で構成される 1 次元の離散時間振動系について考えます。センサーによりこの質量の加速度 $\mathit{a}$ を ${\mathit{F}}_{\mathit{s}}=1$ Hz でサンプリングします。ダンパーによって速度に比例する力が物体に加わるため、物体の運動が妨げられます。減衰定数は $\mathit{b}=0.01$ kg/s です。

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

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

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

$$\begin{array}{c}x(k+1)=Ax(k)+Bu(k),\\ y(k)=Cx(k)+Du(k),\end{array}$$

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

$$A=\exp (A_c\Delta t),B=A_c^{-1}(A-I)B_c,C=\left[\begin{array}{cc}-1& -b\end{array}\right],D=1,$$

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

$$A_c=\left[\begin{array}{cc}0& 1\\ -1& -b\end{array}\right],B_c=\left[\begin{array}{c}0\\ 1\end{array}\right].$$

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)

システムの伝達関数を周波数の関数として推定します。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,mag2db(abs(txy)))
hold on
plot(fz,mag2db(abs(frf)))
hold off
grid
ylim([-60 40])

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

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

MIMO システムの伝達関数の推定

多入力/多出力 (MIMO) システムの伝達関数を推定します。

両側にバネとダンパーが接続された 2 つの質量は、理想的な 1 次元離散時間振動システムを形成します。システム入力配列 u は、質量に適用されるランダムな駆動力で構成されます。システム出力配列 y には、初期の基準位置からの質量の観測変位が格納されます。システムは 40 Hz のレート Fs でサンプリングされます。

MIMO システム入力、システム出力、およびサンプル レートを含むデータ ファイルを読み込みます。この例で使用するデータを生成したシステムは、Frequency-Response Analysis of MIMO Systemの例で解析しています。

load MIMOdata

システム データと関数 tfestimate を使用して、システムの周波数領域の伝達関数を推定し、プロットします。"mimo" オプションを選択して、4 つすべての伝達関数を生成します。$2^{14}$ 個のサンプリング点を使用して離散フーリエ変換を計算し、信号を 5000 サンプルのセグメントに分割し、各セグメントにハン ウィンドウを適用します。隣接するセグメント間に 2500 サンプルのオーバーラップを指定します。

wind = hann(5000);
nfs = 2^14;
nov = 2500;

[tXY,ft] = tfestimate(u,y,wind,nov,nfs,Fs,"mimo");

tiledlayout flow
for jk = 1:2
    for kj = 1:2
        nexttile
        plot(ft,mag2db(abs(tXY(:,jk,kj))))
        grid on
        ylim([-120 0])
        title("Input "+jk+", Output "+kj)
        xlabel("Frequency (Hz)")
        ylabel("Magnitude (dB)")
    end
end

入力引数

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 | single

f — 周波数

ベクトル

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

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

データ型: double | single

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) サイクル/単位時間です。

入力として周波数のベクトル w または f と freqrange を指定した場合、tfestimate は freqrange を無視します。

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_{yx}\left(f\right)}{P_{xx}\left(f\right)},$$

  ここで、P_yx は x と y のクロス パワー スペクトル密度、P_xx は x のパワー スペクトル密度です。この推定は、ノイズがシステム入力と相関していないと仮定しています。

  多入力/多出力 (MIMO) システムの場合、H₁ 推定器は以下のようになります。

  $$H_1(f)=P_{YX}(f)P_{XX}^{-1}(f)=\left[\begin{array}{cccc}{P}_{y_1{x}_1}(f)& P_{y_1{x}_2}(f)& \cdots & P_{y_1{x}_m}(f)\\ P_{y_2{x}_1}(f)& P_{y_2{x}_2}(f)& \cdots & P_{y_2{x}_m}(f)\\ ⋮& ⋮& \ddots & ⋮\\ P_{y_n{x}_1}(f)& P_{y_n{x}_2}(f)& \cdots & P_{y_n{x}_m}(f)\end{array}\right]{\left[\begin{array}{cccc}{P}_{x_1{x}_1}(f)& P_{x_1{x}_2}(f)& \cdots & P_{x_1{x}_m}(f)\\ P_{x_2{x}_1}(f)& P_{x_2{x}_2}(f)& \cdots & P_{x_2{x}_m}(f)\\ ⋮& ⋮& \ddots & ⋮\\ P_{x_m{x}_1}(f)& P_{x_m{x}_2}(f)& \cdots & P_{x_m{x}_m}(f)\end{array}\right]}^{-1}$$

  これは m 入力と n 出力の場合です。ここで、

  • P_yixk は k 番目の入力と i 番目の出力のクロス パワー スペクトル密度です。

  • P_xixk は k 番目および i 番目の入力のクロス パワー スペクトル密度です。

  2 入力と 2 出力の場合、推定器は以下のような行列になります。

  $$H_1(f)=\frac{\left[\begin{array}{cc}{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{array}\right]}{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\left(f\right)=\frac{P_{yy}\left(f\right)}{P_{xy}\left(f\right)},$$

  ここで、P_yy は y のパワー スペクトル密度、P_xy = P^*_yx は x と y のクロス パワー スペクトル密度の複素共役です。この推定は、ノイズがシステム出力と相関していないと仮定しています。

  MIMO システムの場合、H₂ 推定器が良好に定義されるのは、入力と出力が同数 (n = m) の場合のみです。推定器は以下のようになります。

  $$H_2(f)=P_{YY}(f)P_{XY}^{-1}(f)=\left[\begin{array}{cccc}{P}_{y_1{y}_1}(f)& P_{y_1{y}_2}(f)& \cdots & P_{y_1{y}_n}(f)\\ P_{y_2{y}_1}(f)& P_{y_2{y}_2}(f)& \cdots & P_{y_2{y}_n}(f)\\ ⋮& ⋮& \ddots & ⋮\\ P_{y_n{y}_1}(f)& P_{y_n{y}_2}(f)& \cdots & P_{y_n{y}_n}(f)\end{array}\right]{\left[\begin{array}{cccc}{P}_{x_1{y}_1}(f)& P_{x_1{y}_2}(f)& \cdots & P_{x_1{y}_n}(f)\\ P_{x_2{y}_1}(f)& P_{x_2{y}_2}(f)& \cdots & P_{x_2{y}_n}(f)\\ ⋮& ⋮& \ddots & ⋮\\ P_{x_n{y}_1}(f)& P_{x_n{y}_2}(f)& \cdots & P_{x_n{y}_n}(f)\end{array}\right]}^{-1},$$

  ここで、

  • P_yiyk は k 番目および i 番目の出力のクロス パワー スペクトル密度です。

  • P_xiyk は 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++ コードを生成します。

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

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

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

この関数は、GPU 配列を完全にサポートします。詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

R2024a: 単精度可変サイズ ウィンドウ入力のためのコード生成サポート

関数 tfestimate は、コード生成用の単精度可変サイズ ウィンドウ入力をサポートします。

R2023b: 単精度データと gpuArray オブジェクトの使用

関数 tfestimate は、単精度入力と gpuArray オブジェクトをサポートします。gpuArray オブジェクトを使用するには、Parallel Computing Toolbox™ が必要です。

参考

cpsd | mscohere | periodogram | pwelch