ドキュメンテーション

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

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 に指定して伝達関数の推定 txy を求めます。

  • xy の両方がベクトルの場合、それらは同じ長さでなければなりません。

  • 信号の 1 つが行列でもう 1 つがベクトルの場合は、ベクトルの長さが行列の行数と等しくなければなりません。関数はベクトルを展開し、列ごとの伝達関数推定の行列を返します。

  • xy が同じ行数だが異なる列数をもつ行列の場合、txy は多入力/多出力 (MIMO) 伝達関数になり、すべての入力信号と出力信号を結合します。txy は 3 次元配列です。x が m 列をもち、y が n 列をもつ場合、txy は n 列および m ページをもちます。詳細は、伝達関数を参照してください。

  • xy が等しいサイズの行列の場合、tfestimate は列方向に動作します (txy(:,n) = tfestimate(x(:,n),y(:,n)))。MIMO 推定を取得するには、引数リストに 'mimo' を追加します。

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

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 を返し、これにより、伝達関数が推定されます。fstfestimate に対する 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 つのシーケンス x および y 間の伝達関数推定を計算して、プロットします。シーケンス x はガウス ノイズから構成されます。y は、x を正規化されたカットオフ周波数 ラジアン/サンプルをもつ 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)

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

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

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

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

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

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

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

2000 個の時間サンプルを生成します。サンプリング間隔は と定義します。

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

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

ここで、 は状態ベクトル、 および はそれぞれ質量の位置と速度、 は駆動力、 は測定出力です。状態空間行列は次のようになります。

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

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,20*log10(abs(txy)))
hold on
plot(fz,20*log10(abs(frf)))
hold off
grid
ylim([-60 40])

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

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

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

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

30000 回のサンプルを生成しますが、これは 600 秒に相当します。サンプリング間隔は と定義します。

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

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

ここで、 は状態ベクトル、 および はそれぞれ 番目の質量の位置と速度、 は入力駆動力のベクトル、 は出力ベクトルです。状態空間行列は次のようになります。

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

に設定します。

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 点 個を指定します。推定をプロットします。

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

伝達関数は期待値 において最大値をもちます。ここで、 はモードの行列の固有値です。

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

    3.8470
   14.4259

を設定することでシステムに減衰を追加します。減衰システムの時間発展を同じ駆動力で計算します。MIMO 伝達関数の の推定を同じウィンドウとオーバーラップを使用して計算します。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')

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

入力引数

すべて折りたたむ

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

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

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

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

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

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

  • window が整数の場合、tfestimatexy を長さ window のセグメントに分割し、各セグメントにその長さのハミング ウィンドウを適用します。

  • window がベクトルの場合、tfestimatexy をベクトルと同じ長さのセグメントに分割し、window を使用して各セグメントにウィンドウを適用します。

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

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

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

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

データ型: single | double

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

  • window がスカラーの場合、noverlapwindow より小さくなければなりません。

  • window がベクトルの場合、noverlapwindow の長さより小さくなければなりません。

noverlap を空として指定した場合、tfestimate はセグメント間で 50% のオーバーラップが発生する数を使用します。セグメントの長さを指定していない場合、関数により noverlap が ⌊N/4.5⌋ に設定されます。ここで、N は入力および出力信号の長さです。

データ型: double | single

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

データ型: single | double

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

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

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

データ型: double

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

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

データ型: double

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

  • 'onesided' — 2 つの実数値の入力信号 xy の間の伝達関数の片側推定を返します。nfft が偶数の場合、txynfft/2 + 1 行で、計算区間は [0,π] ラジアン/サンプルです。nfft が奇数の場合、txy は (nfft + 1)/2 行で、計算区間は [0,π) ラジアン/サンプルです。fs を指定すると、nfft が偶数の場合は、対応する計算区間は [0,fs/2] サイクル/単位時間で、nfft が奇数の場合は、[0,fs/2) サイクル/単位時間です。

  • 'twosided' — 実数値または複素数値の 2 つの入力信号 xy の間の伝達関数の両側推定を返します。この場合、txynfft 行で、計算区間は [0,2π) ラジアン/サンプルです。fs を指定した場合、計算区間は [0,fs) サイクル/単位時間となります。

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

データ型: char

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

  • ノイズが入力信号と無相関な場合は、'H1' を使用します。

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

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

出力引数

すべて折りたたむ

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

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

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

詳細

すべて折りたたむ

伝達関数

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

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

    H1(f)=Pyx(f)Pxx(f),

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

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

    H1(f)=PYX(f)PXX1(f)=[Py1x1(f)Py1x2(f)Py1xm(f)Py2x1(f)Py2x2(f)Py2xm(f)Pynx1(f)Pynx2(f)Pynxm(f)][Px1x1(f)Px1x2(f)Px1xm(f)Px2x1(f)Px2x2(f)Px2xm(f)Pxmx1(f)Pxmx2(f)Pxmxm(f)]1

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

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

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

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

    H1(f)=[Py1x1(f)Px2x2(f)Py1x2(f)Px2x1(f)Py1x2(f)Px1x1(f)Py1x1(f)Px1x2(f)Py2x1(f)Px2x2(f)Py2x2(f)Px2x1(f)Py2x2(f)Px1x1(f)Py2x1(f)Px1x2(f)]Px1x1(f)Px2x2(f)Px1x2(f)Px2x1(f).

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

    H2(f)=Pyy(f)Pxy(f),

    ここで、Pyyy のパワー スペクトル密度、Pxy = P*yxxy のクロス パワー スペクトル密度の複素共役です。この推定は、ノイズがシステム出力と相関をもたないと仮定しています。

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

    H2(f)=PYY(f)PXY1(f)=[Py1y1(f)Py1y2(f)Py1yn(f)Py2y1(f)Py2y2(f)Py2yn(f)Pyny1(f)Pyny2(f)Pynyn(f)][Px1y1(f)Px1y2(f)Px1yn(f)Px2y1(f)Px2y2(f)Px2yn(f)Pxny1(f)Pxny2(f)Pxnyn(f)]1,

    ここで、

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

    • Pxiyk は 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.

R2006a より前に導入