ドキュメンテーション

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

modalfrf

モード解析の周波数応答関数

説明

frf = modalfrf(x,y,fs,window) は、周波数応答関数 frf の行列を励起信号 x および応答信号 y から推定します。これらはすべてレート fs でサンプリングされます。出力 frfH1 の推定で、ウェルチ法を使用して、window で信号をウィンドウ処理して計算されます。xy は同じ行数でなければなりません。x または y が行列の場合、各列が信号を表します。周波数応答関数の行列 frf は動的柔軟性に関して計算され、システム応答 y には加速度の測定値が含まれます。

frf = modalfrf(x,y,fs,window,noverlap) は、隣り合ったセグメント間のオーバーラップを noverlap サンプルに指定します。

frf = modalfrf(___,Name,Value) は、前の構文の入力の任意の組み合わせを使用して、名前と値のペアの引数によりオプションを指定します。オプションには推定器、測定設定、およびシステム応答を測定するセンサーのタイプが含まれます。

[frf,f,coh] = modalfrf(___) はさらに、各周波数応答関数に対応する周波数ベクトル、および多重コヒーレンス行列を返します。

[frf,f] = modalfrf(sys) は、識別されたモデル sys の周波数応答関数を計算します。ssestn4sidtfest などの推定コマンドを使用して、時間領域の入力信号と出力信号から sys を作成します。この構文では、名前と値のペアの引数 'Sensor' 以外は使用できません。この構文を使用するには、System Identification Toolbox™ のライセンスを所有していなければなりません。

frf = modalfrf(sys,f) では、計算する周波数 frf を指定します。この構文では、名前と値のペアの引数 'Sensor' 以外は使用できません。この構文を使用するには、System Identification Toolbox のライセンスを所有していなければなりません。

出力引数なしで modalfrf(___) を使用すると、周波数応答関数が現在の Figure にプロットされます。プロットは最初の 4 つの励起と 4 つの応答に制限されます。

すべて折りたたむ

単入力/単出力のハンマー励起の周波数応答関数を可視化します。

以下を含むデータ ファイルを読み込みます。

  • Xhammer 周期的に与えられる 5 回のハンマー打撃から成る入力励起信号。

  • Yhammer 入力に対するシステムの応答。Yhammer は変位として測定されます。

この信号は 4 kHz でサンプリングされています。励起および出力信号をプロットします。

load modaldata

subplot(2,1,1)
plot(thammer,Xhammer(:))
ylabel('Force (N)')
subplot(2,1,2)
plot(thammer,Yhammer(:))
ylabel('Displacement (m)')
xlabel('Time (s)')

周波数応答関数を計算して表示します。箱型ウィンドウを使用して信号にウィンドウを適用します。ウィンドウがハンマー打撃間の区間をカバーするように指定します。

clf
winlen = size(Xhammer,1);
modalfrf(Xhammer(:),Yhammer(:),fs,winlen,'Sensor','dis')

ランダム ノイズで励起される 2 入力/2 出力システムの周波数応答関数を計算します。

入力励起信号 Xrand、システム応答 Yrand を含むデータ ファイルを読み込みます。5000 サンプルのハン ウィンドウおよび隣接するデータ セグメント間に 50% のオーバーラップを使用して、周波数応答関数を計算します。出力測定値を変位として指定します。

load modaldata
winlen = 5000;

frf = modalfrf(Xrand,Yrand,fs,hann(winlen),0.5*winlen,'Sensor','dis');

modalfrf のプロット機能を使用して、応答を可視化します。

modalfrf(Xrand,Yrand,fs,hann(winlen),0.5*winlen,'Sensor','dis')

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

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

3000 個の時間サンプルを生成します。サンプリング間隔は Δt=1/Fs と定義します。

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

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

x(k+1)=Ax(k)+Bu(k),y(k)=Cx(k)+Du(k),

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

A=exp(AcΔt),B=Ac-1(A-I)Bc,C=[10],D=0,

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

Ac=[01-1-b],Bc=[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

固有振動数を 1/2π と比較します。これは非減衰システムの理論値です。

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

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

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

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

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

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

x(k+1)=Ax(k)+Bu(k),y(k)=Cx(k)+Du(k),

ここで、x=[r1v1r2v2]T は状態ベクトル、ri および vi はそれぞれ i 番目の質量の位置と速度、u=[u1u2]T は入力駆動力のベクトル、y=[r1r2]T は出力ベクトルです。状態空間行列は次のようになります。

A=exp(AcΔt),B=Ac-1(A-I)Bc,C=[10000010],D=[0000],

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

Ac=[0100-2k-2bkb0001k/μb/μ-2k/μ-2b/μ],Bc=[00100001/μ].

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

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 入力/6 出力データセットの周波数応答関数を計算します。

入力励起と出力加速度計測定値を含む構造体を読み込みます。システムは、1024 Hz で 3.9 秒間サンプリングされます。

load modaldata SteelFrame
X = SteelFrame.Input;
Y = SteelFrame.Output;
fs = SteelFrame.Fs;

部分空間法を使用して周波数応答関数を計算します。入力信号と出力信号をオーバーラップのない 1000 サンプルのセグメントに分割します。箱型ウィンドウを使用して各セグメントにウィンドウを適用します。モデル次数を 36 に指定します。

[frf,f] = modalfrf(X,Y,fs,1000,'Estimator','subspace','Order',36);

システムの安定化ダイアグラムを可視化します。最大 15 の物理モードを特定します。

modalsd(frf,f,fs,'MaxModes',15)

入力引数

すべて折りたたむ

励起信号。ベクトルまたは行列で指定します。

データ型: single | double

応答信号。ベクトルまたは行列で指定します。

データ型: single | double

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

データ型: single | double

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

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

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

  • 'Estimator''subspace' を指定した場合、modalfrfwindow の形状は無視し、その長さを使用して返される周波数応答関数の周波数点数を決定します。

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

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

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

データ型: single | double

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

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

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

データ型: double | single

識別されたシステム。識別されたパラメーターをもつモデルとして指定します。ssestn4sidtfest などの推定コマンドを使用して、時間領域の入力信号と出力信号から sys を作成します。例は、識別したモデルのモード解析を参照してください。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 でサンプリングされます。

周波数。Hz 表記のベクトルとして指定します。

データ型: single | double

名前と値のペアの引数

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

例: 'Sensor','vel','Est','H1' を指定すると、応答信号が速度の測定値から構成され、推定器として H1 が選択されます。

推定器。'Estimator' と、'H1''H2''Hv''subspace' のいずれかから成るコンマ区切りのペアとして指定します。H1H2 の推定器の詳細は、伝達関数を参照してください。

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

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

  • 'Hv' を使用すると、誤りを含む行列のトレースを最小化して、モデル化された応答データと推定された応答データ間の相違を最小化します。HvH1H2 の幾何平均です。Hv = (H1H2)1/2

    測定は、単入力/単出力 (SISO) でなければなりません。

  • 'subspace' を使用して状態空間モデルを使用する周波数応答関数を計算します。この場合には、noverlap 引数は無視されます。この手法では、通常は、ノンパラメトリックによる方法より必要なデータは少なくなります。詳細は、n4sid を参照してください。

状態空間モデルの直達の存在。'Feedthrough' と logical で構成されるコンマ区切りのペアとして指定します。この引数は、'Estimator''subspace' に指定した場合にのみ使用可能です。

データ型: logical

同数の励起および応答チャネルの測定設定。'Measurement' と、'fixed''rovinginput' または 'rovingoutput' から成るコンマ区切りのペアとして指定します。

  • 'fixed' は、励起ソースとセンサーがシステムの固定位置にあるときに使用します。各励起が各応答に寄与します。

  • 'rovinginput' は、測定値がロービング励起 (または "ロービング ハンマー") テストの結果であるときに使用します。単一のセンサーがシステムの固定位置に保持されます。単一の励起ソースが複数の位置に配置され、位置ごとに 1 つのセンサー応答を生成します。関数は frf(:,:,i) = modalfrf(x(:,i),y(:,i)) を出力します。

  • 'rovingoutput' は、測定値が "ロービング センサー" テストの結果であるときに使用します。単一の励起ソースがシステムの固定位置に保持されます。単一のセンサーが複数の位置に配置され、位置ごとに 1 つの励起に応答します。関数は frf(:,i) = modalfrf(x(:,i),y(:,i)) を出力します。

状態空間モデル次数。'Order' と、整数または整数の行ベクトルで構成されるコンマ区切りのペアとして指定します。整数から成るベクトルを指定すると、関数は指定された範囲から最適な次数値を選択します。この引数は、'Estimator''subspace' に指定した場合にのみ使用可能です。

データ型: single | double

センサー タイプ。'Sensor' と、'acc''dis' または 'vel' から成るコンマ区切りのペアとして指定します。

  • 'acc' — 応答信号電圧が加速度に比例します。

  • 'dis' — 応答信号電圧が変位に比例します。

  • 'vel' — 応答信号電圧が速度に比例します。

出力引数

すべて折りたたむ

周波数応答関数。ベクトル、行列または 3 次元配列で返されます。frf のサイズは p × m × n です。ここで、p は周波数ビンの数、m は応答の数、n は励起信号の数です。

ベクトルとして返される周波数。

多重コヒーレンス行列。行列として返されます。coh は応答信号ごとに 1 つの列をもちます。

参照

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

[2] 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.

参考

| | |

トピック

R2017a で導入