ドキュメンテーション

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

modalfit

周波数応答関数からのモーダル パラメーター

構文

fn = modalfit(frf,f,fs,mnum)
fn = modalfit(frf,f,fs,mnum,Name,Value)
[fn,dr,ms] = modalfit(___)
[fn,dr,ms,ofrf] = modalfit(___)
[___] = modalfit(sys,f,mnum,Name,Value)

説明

fn = modalfit(frf,f,fs,mnum) は、システムの mnum モードの固有周波数を推定します。システムは、周波数 f およびサンプルレート fs で定義される周波数応答関数 frf を測定したものです。

fn = modalfit(frf,f,fs,mnum,Name,Value) は、名前と値のペアの引数を使用して追加オプションを指定します。

[fn,dr,ms] = modalfit(___) は、さらに、fn で各固有周波数に対応する減衰比とモード形状ベクトルを返します。前の構文の入力の任意の組み合わせを使用します。

[fn,dr,ms,ofrf] = modalfit(___) は、さらに、推定されたモーダル パラメーターに基づいて再構成された周波数応答関数の配列を返します。

[___] = modalfit(sys,f,mnum,Name,Value) は識別されたモデル sys のモーダル パラメーターを推定します。ssesttfest などの推定コマンドを使用して、測定された周波数応答関数から、または時間領域の入力信号および出力信号から sys を作成します。この構文では、名前と値のペアの引数 'DriveIndex''FreqRange' および 'PhysFreq' を指定できます。通常は、ノンパラメトリック法を使用する構文より必要なデータは少なくなります。この構文を使用するには、System Identification Toolbox™ のライセンスを所有していなければなりません。

すべて折りたたむ

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

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

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

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

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

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

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

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

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

C = [1 0];
D = 0;

質量はランダム入力でまず 2000 秒間駆動され、次に静止状態に戻るまで放置されます。状態空間モデルを使用して、すべてゼロの初期状態からの系の時間発展を計算します。質量の変位を時間の関数としてプロットします。

rng default
u = randn(1,N)/2;
u(2001:end) = 0;

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)

システムのモーダル周波数応答関数を推定します。測定信号の 1/2 の長さでハン ウィンドウを使用します。出力を質量の変位として指定します。

wind = hann(N/2);

[frf,f] = modalfrf(u',y',Fs,wind,'Sensor','dis');

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

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

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

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

振動モードの固有周波数と減衰比を推定します。

[fn,dr] = modalfit(frf,f,Fs,1,'FitMethod','PP')
fn = 0.1593
dr = 0.0043

固有周波数を と比較します。これは非減衰システムの理論値です。

theo = 1/(2*pi)
theo = 0.1592

宇宙ステーション モジュールの周波数応答関数 (FRF) 配列からそのモーダル パラメーターを計算します。

3 入力/3 出力の FRF 配列を含む構造体を読み込みます。そのシステムは、320 Hz でサンプリングされます。

load modaldata SpaceStationFRF

frf = SpaceStationFRF.FRF; 
f = SpaceStationFRF.f; 
fs = SpaceStationFRF.Fs;

最小二乗有理関数法を使用して最小 24 のモードのモーダル パラメーターを抽出します。

[fn,dr,ms,ofrf] = modalfit(frf,f,fs,24,'FitMethod','lsrf');

再構成した FRF 配列を測定した配列と比較します。

for ij = 1:3
 for ji = 1:3
    subplot(3,3,3*(ij-1)+ji)
    loglog(f,abs(frf(:,ji,ij)))
    hold on
    loglog(f,abs(ofrf(:,ji,ij)))
    hold off
    axis tight
    title(sprintf('In%d -> Out%d',ij,ji))
    if ij==3
        xlabel('Frequency (Hz)')
    end
 end
end

シンプルな多入力/多出力システムの周波数応答関数とモーダル パラメーターを推定します。

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

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

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

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

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

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

に設定します。

k = 400;
b = 0.1;
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 = [1 0 0 0;0 0 1 0];
D = zeros(2);

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

rng default
u = randn(2,N);

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

入力および出力データを使用して、システムの伝達関数を周波数の関数として推定します。隣接するセグメント間で 9000 サンプルがオーバーラップする 15000 サンプルのハン ウィンドウを使用します。測定出力を変位として指定します。

wind = hann(15000);
nove = 9000;
[FRF,f] = modalfrf(u',y',Fs,wind,nove,'Sensor','dis');

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

nfs = 2048;
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(f,20*log10(abs(FRF(:,jk,kj))))
        hold on
        plot(fz*Fs,20*log10(abs(frf(jk,:,kj))))
        hold off
        axis([0 Fs/2 -100 0])
        title(sprintf('Input %d, Output %d',jk,kj))
    end
end

出力引数なしで modalfrf の構文を使用して推定をプロットします。

figure
modalfrf(u',y',Fs,wind,nove,'Sensor','dis')

システムの固有周波数、減衰比、モード形状を推定します。計算にはピーク選択法を使用します。

[fn,dr,ms] = modalfit(FRF,f,Fs,2,'FitMethod','pp');
fn
fn = 
fn(:,:,1) =

    3.8466    3.8466
    3.8495    3.8495


fn(:,:,2) =

    3.8492    3.8490
    3.8552   14.4684

固有周波数を非減衰システムの理論予測と比較します。

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

    3.8470
   14.4259

いくつかのバーストのランダム ノイズで励起される 2 入力/3 出力システムの固有周波数、減衰比、モード形状を計算します。各バーストは 1 秒間持続し、各バーストの終了と次の開始の間は 2 秒間です。データは 4 kHz でサンプリングされています。

データ ファイルを読み込みます。入力信号と出力信号をプロットします。

load modaldata

subplot(2,1,1)
plot(Xburst)
title('Input Signals')
subplot(2,1,2)
plot(Yburst)
title('Output Signals')

周波数応答関数を計算します。長さがバースト周期に等しく、隣接するセグメント間のオーバーラップがない箱型ウィンドウを指定します。

burstLen = 12000;
[frf,f] = modalfrf(Xburst,Yburst,fs,burstLen);

安定化ダイアグラムを可視化して、安定した固有周波数を返します。最大モデル次数として 30 モードを指定します。

figure
modalsd(frf,f,fs,'MaxModes',30);

プロットを拡大します。平均応答関数は 373 Hz、852 Hz、1371 Hz で最大値をもち、これはシステムの物理周波数に対応します。最大値を変数に保存します。

phfr = [373 852 1371];

最小二乗複素指数 (LSCE) アルゴリズムを使用してモーダル パラメーターを計算します。モデル次数として 6 モードを指定し、安定化ダイアグラムから決定された物理周波数として 3 モードを指定します。関数により、各入力基準に対して固有周波数と減衰比の 1 つのセットが生成されます。

[fn,dr,ms,ofrf] = modalfit(frf,f,fs,6,'PhysFreq',phfr);

再構成された周波数応答関数をプロットし、元のものと比較します。

for k = 1:2
    for m = 1:3
        subplot(2,3,m+3*(k-1))
        plot(f/1000,10*log10(abs(frf(:,m,k))))
        hold on
        plot(f/1000,10*log10(abs(ofrf(:,m,k))))
        hold off
        text(1,-50,[['Output ';' Input '] num2str([m k]')])
        ylim([-100 -40])
    end
end
subplot(2,3,2)
title('Frequency-Response Functions')

入力引数

すべて折りたたむ

周波数応答関数。ベクトル、行列または 3 次元配列で指定します。frf のサイズは p × m × n です。ここで、p は周波数ビンの数、m は応答信号の数、n は伝達関数を推定するための励起信号の数です。

例: tfestimate(randn(1,1000),sin(2*pi*(1:1000)/4)+randn(1,1000)/10) は、発振器の周波数応答を近似します。

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

周波数。ベクトルとして指定します。f の要素数は、frf の行数と等しくなければなりません。

データ型: single | double

測定値データのサンプルレート。ヘルツ単位の正のスカラーとして指定します。

データ型: single | double

モード数。正の整数で指定します。

データ型: single | double

識別されたシステム。識別されたパラメーターをもつモデルとして指定します。ssestn4sidtfest などの推定コマンドを使用して、測定された周波数応答関数から、または時間領域の入力信号および出力信号から sys を作成します。例は、識別したモデルのモード解析を参照してください。この入力引数を使用するには、System Identification Toolbox のライセンスを所有していなければなりません。

例: idss([0.5418 0.8373;-0.8373 0.5334],[0.4852;0.8373],[1 0],0,[0;0],[0;0],1) は、識別された状態空間モデルを生成します。生成されるモデルは、単位弾性定数をもつバネと定数 0.01 のダンパーで壁につながれた単位質量に対応します。質量の変位は 1 Hz でサンプリングされます。

例: idtf([0 0.4582 0.4566],[1 -1.0752 0.99],1) は、識別された伝達関数モデルを生成します。生成されるモデルは、単位弾性定数をもつバネと定数 0.01 のダンパーで壁につながれた単位質量に対応します。質量の変位は 1 Hz でサンプリングされます。

名前と値のペアの引数

オプションの引数 Name,Value のコンマ区切りペアを指定します。Name は引数名で、Value は対応する値です。Name は引用符で囲まなければなりません。Name1,Value1,...,NameN,ValueN のように、複数の名前と値のペアの引数を、任意の順番で指定できます。

例: 'FitMethod','pp','FreqRange',[0 500] はピーク選択法を使用して近似を実行し、周波数範囲を 0 ~ 500 Hz に制限します。

推定伝達関数の直達の存在。'Feedthrough' と logical で構成されるコンマ区切りのペアとして指定します。この引数は、'FitMethod''lsrf' に指定した場合にのみ使用可能です。

データ型: logical

近似アルゴリズム。'FitMethod''lsce''lsrf''pp' のいずれかで構成されるコンマ区切りのペアとして指定します。

  • 'lsce'最小二乗複素指数法.'lsce' を指定した場合、fnmnum 個の要素をもつベクトルで、frf のサイズとは関係ありません。

  • 'lsrf' — 最小二乗有理関数推定法。'lsrf' を指定した場合、fnmnum 個の要素をもつベクトルで、frf のサイズとは関係ありません。この手法については、[3]で説明します。詳細は、Continuous-Time Transfer Function Estimation Using Continuous-Time Frequency-Domain Data (System Identification Toolbox)を参照してください。このアルゴリズムは通常、ノンパラメトリック法を使用する方法より必要なデータが少なく、不均一な f に対して有効な唯一のアルゴリズムです。

  • 'pp'ピーク選択法.frf が n 個の励起信号と m 個の応答信号から計算された場合、fnmnum × m × n の配列で、frf ごとに 1 つの fn の推定および 1 つの dr の推定をもちます。

周波数範囲。コンマ区切りペアとして指定します。ペアは 'FreqRange' および f で指定した範囲内の増加する正の値から成る 2 要素ベクトルで構成されます。

データ型: single | double

解析に含める物理モードの固有周波数。コンマ区切りペアとして指定します。ペアは 'PhysFreq' および f でカバーされる範囲内の周波数値のベクトルで構成されます。関数は、ベクトルで指定された値に最も近い固有周波数をもつモードを解析に含めます。ベクトルに m の周波数値が含まれる場合、fndr はそれぞれ m 行をもち、ms は m 列をもちます。この引数を指定しない場合、関数では f の周波数範囲全体が使用されます。

データ型: single | double

駆動点の周波数応答関数のインデックス。'DriveIndex' と正の整数の 2 要素ベクトルから成るコンマ区切りペアとして指定します。ベクトルの最初の要素は、システム応答の数以下でなければなりません。ベクトルの 2 番目の要素は、システム励起の数以下でなければなりません。モード形状は駆動点に基づいて単一モーダルに正規化されます。

例: 'DriveIndex',[2 3] を指定すると、駆動点の周波数応答関数が frf(:,2,3) になります。

データ型: single | double

出力引数

すべて折りたたむ

固有周波数。行列または 3 次元配列として返されます。fn のサイズは、'FitMethod' で指定する近似アルゴリズムの選択に依存します。

  • 'lsce' または 'lsrf' を指定した場合、fnmnum 個の要素をもつベクトルで、frf のサイズとは関係ありません。システムの振動モードが mnum より多い場合、'lsrf' メソッドは固有周波数の昇順で並べ替えた最初の mnum の最小減衰モードを返します。

  • 'pp' を指定すると、fn はサイズが mnum × m × n の配列になり、frf ごとに 1 つの fn の推定および 1 つの dr の推定をもちます。

fn の固有周波数の減衰比。fn と同じサイズの行列または 3 次元配列として返されます。

モード形状ベクトル。行列として返されます。msmnum 列をもち、それぞれに長さ q のモード形状ベクトルが含まれます。ここで、q は励起チャネル数と応答チャネル数のうち大きいほうです。

再構成された周波数応答関数。frf と同じサイズのベクトル、行列または 3 次元配列として返されます。

アルゴリズム

すべて折りたたむ

最小二乗複素指数法

最小二乗複素指数法は、各周波数応答関数に対応するインパルス応答を計算し、Prony 法を使用して、複素減衰正弦波のセットを応答に対して近似します。

サンプリングされた減衰正弦波は以下の形でキャストできます。

si(n)=Aiebin/fscos(2πfin/fs+ϕi)=12Aiejϕiexp((bi/fsj2πfi/fs)n)+12Aiejϕiexp((bi/fs+j2πfi/fs)n)ai+xi+n+aixin,

ここで、

  • fs はサンプルレート。

  • fi は正弦波周波数。

  • bi は減衰係数。

  • Ai および ϕi は正弦波の振幅と位相。

ai は "振幅" と呼ばれ、xi は "極" と呼ばれます。Prony 法はサンプリングされた関数 h(n)N/2 モードの重ね合わせとして表します (したがって、N の振幅と極)。

h(0)=a1x10+a2x20+aNxN0h(1)=a1x11+a2x21++aNxN1h(N1)=a1x1N1+a2x2N1++aNxNN1.

極は多項式の根で係数 c0, c1, …, cN–1 をもちます。

xiN+cN1xiN1++c1xi1+c0xi0=0.

係数を求めるには、h の L = 2N サンプルの自己回帰モデルを使用します。

[h(0)h(1)h(N1)h(1)h(2)h(N)h(LN1)h(LN)h(L2)][c0c1cN1]=[h(N)h(N+1)h(L1)].

極を求めるために、アルゴリズムで関数 roots が使用されます。極を求めたら、周波数と減衰係数を決定するため、極の対数の虚数部と実数部を特定します。最後に、以下を使用して、振幅の解を求めてインパルス応答を再構成します。

[h(0)h(N1)]=[x10xN0x1N1xNN1][a1aN].

以下の単純な MATLAB® 実装でこの手続きの概要を示します。

N = 4;
L = 2*N;
h = rand(L,1);
c = hankel(h(1:N),h(L-N:L-1))\-h(N+1:L);
x = roots([1;c(N:-1:1)]).';
p = log(x);
hrec = x.^((0:L-1)')*(x.^((0:L-1)')\h(1:L));
sum(h-hrec)
ans =

   3.2613e-15 - 1.9297e-16i
複数の周波数応答関数のサンプルを含めるようにシステムを構成し、最小二乗を使用して解を求めることもできます。

ピーク選択法

ピーク選択法では、周波数応答関数のそれぞれの大きなピークが、厳密に 1 つの固有モードに対応すると仮定します。ピーク付近では、システムは自由度 1 の減衰調和振動子のように振る舞うと仮定します。

H(f)=1(2π)21/mf2+j2ζrfrf+fr2H(f)fr2+j2ζrfrfH(f)1(2π)2m=f2H(f),

ここで、H は周波数応答関数、fr は非減衰共振周波数、ζr = b/(4mk)1/2 は相対減衰、b は減衰定数、k は弾性定数、m は質量です。

fp に位置するピークを与え、この手順では、いずれかの側にピークと固定数の点を受け取り、質量項をダミー変数 d に置き換えて、以下の方程式系を解いてモーダル パラメーターを求めます。

[H(fpk)j2fpkH(fpk)1H(fp)j2fpH(fp)1H(fp+k)j2fp+kH(fp+k)1][fr2ζrfrd]=[fpk2H(fpk)fp2H(fp)fp+k2H(fp+k)].

参照

[1] Allemang, Randall J., and David L. Brown. “Experimental Modal Analysis and Dynamic Component Synthesis, Vol. III: Modal Parameter Estimation.” Technical Report AFWAL-TR-87-3069. Air Force Wright Aeronautical Laboratories, Wright-Patterson Air Force Base, OH, December 1987.

[2] Brandt, Anders. Noise and Vibration Analysis: Signal Analysis and Experimental Procedures. Chichester, UK: John Wiley & Sons, 2011.

[3] Ozdemir, Ahmet Arda, and Suat Gumussoy. "Transfer Function Estimation in System Identification Toolbox via Vector Fitting." Proceedings of the 20th World Congress of the International Federation of Automatic Control, Toulouse, France, July 2017.

参考

| | | |

トピック

R2017a で導入