このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

modalfrf

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

構文

frf = modalfrf(x,y,fs,window)

frf = modalfrf(x,y,fs,window,noverlap)

frf = modalfrf(___,Name,Value)

[frf,f,coh] = modalfrf(___)

[frf,f] = modalfrf(sys)

frf = modalfrf(sys,f)

modalfrf(___)

説明

frf = modalfrf(x,y,fs,window) は、周波数応答関数 frf の行列を励起信号 x および応答信号 y から推定します。これらはすべてレート fs でサンプリングされます。出力 frf は H₁ の推定で、ウェルチ法を使用して、window で信号をウィンドウ処理して計算されます。x と y は同じ行数でなければなりません。x または y が行列の場合、各列が信号を表します。

システム応答 y には加速度測定値が含まれるものとします。周波数応答関数を変位測定値または速度測定値から計算するには、引数 'Sensor' を使用します。modalfrf は、センサーの種類に関係なく、周波数応答関数を常に動的柔軟性 (レセプタンス) 形式で出力します。

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

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

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

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

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

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

例

ハンマー励起の周波数応答関数

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

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

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

• Yhammer $$\text{—}$$ 入力に対するシステムの応答。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')

MIMO 周波数応答関数

ランダム ノイズで励起される 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')

SISO システムの周波数応答関数

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

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

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

Fs = 1;
dt = 1/Fs;
N = 3000;
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{r}$ は測定出力です。状態空間行列は次のようになります。

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

$\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(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\pi$ と比較します。これは非減衰システムの理論値です。

theo = 1/(2*pi)

theo = 0.1592

2 つの剛体の振動子のモーダル パラメーター

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

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

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

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

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

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

ここで、$$x={\left[\begin{array}{cccc}{r}_1& v_1& r_2& v_2\end{array}\right]}^T$$ は状態ベクトル、$$r_i$$ および $$v_i$$ はそれぞれ $$i$$ 番目の質量の位置と速度、$$u={\left[\begin{array}{cc}{u}_1& u_2\end{array}\right]}^T$$ は入力駆動力のベクトル、および $$y={\left[\begin{array}{cc}{r}_1& r_2\end{array}\right]}^T$$ は出力ベクトルです。状態空間行列は次のようになります。

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

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

$$A_c=\left[\begin{array}{cccc}0& 1& 0& 0\\ -2k& -2b& k& b\\ 0& 0& 0& 1\\ k/\mu & b/\mu & -2k/\mu & -2b/\mu \end{array}\right],B_c=\left[\begin{array}{cc}0& 0\\ 1& 0\\ 0& 0\\ 0& 1/\mu \end{array}\right].$$

$\mathit{k}=400$、$\mathit{b}=0.1$、および $\mu =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)

入力引数

x — 励起信号

ベクトル | 行列

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

データ型: single | double

y — 応答信号

ベクトル | 行列

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

データ型: single | double

fs — サンプル レート

正のスカラー

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

データ型: single | double

window — ウィンドウ

整数 | ベクトル

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

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

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

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

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

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

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

データ型: single | double

noverlap — オーバーラップするサンプル数

0 (既定値) | 正の整数

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

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

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

データ型: double | single

sys — 識別されたシステム

識別されたパラメーターをもつモデル

識別されたシステム。識別されたパラメーターをもつモデルとして指定します。ssest (System Identification Toolbox) や n4sid (System Identification Toolbox)、tfest (System Identification Toolbox) などの推定コマンドを使用して、時間領域の入力信号と出力信号から 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 でサンプリングされます。

f — 周波数

ベクトル

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

データ型: single | double

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで、Name は引数名で、Value は対応する値です。名前と値の引数は他の引数の後に指定しなければなりませんが、ペアの順序は重要ではありません。

R2021a より前では、コンマを使用して名前と値をそれぞれ区切り、Name を引用符で囲みます。

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

Estimator — 推定器

'H1' (既定値) | 'H2' | 'Hv' | 'subspace'

推定器。'H1'、'H2'、'Hv'、または 'subspace' で指定します。H₁ と H₂ の推定器の詳細については、伝達関数を参照してください。

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

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

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

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

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

Feedthrough — 状態空間モデルの直達の存在

false (既定値) | true

状態空間モデルの直達の存在。logical 値として指定します。この引数は、'Estimator' を 'subspace' に指定した場合にのみ使用可能です。

データ型: logical

Measurement — 測定コンフィギュレーション

'fixed' (既定値) | 'rovinginput' | 'rovingoutput'

同数の励起および応答チャネルの測定設定。'fixed'、'rovinginput'、または 'rovingoutput' として指定します。

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

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

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

Order — 状態空間モデル次数

1:10 (既定値) | 整数 | 整数の行ベクトル

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

データ型: single | double

Sensor — センサー タイプ

'acc' (既定値) | 'dis' | 'vel'

センサー タイプ。'acc'、'vel'、または 'dis' として指定します。

• 'acc' — システムの応答信号が加速度に比例するように指定します。

• 'vel' — システムの応答信号が速度に比例するように指定します。

• 'dis' — システムの応答信号が変位に比例するように指定します。

modalfrf は、センサーの種類に関係なく、周波数応答関数を常に動的柔軟性 (レセプタンス) 形式で出力します。

例:非減衰高調波発振器

${\mathit{f}}_{\mathrm{s}}=1/\Delta \mathit{t}$ のレートでサンプリングされた単位質量と弾性定数の単純な非減衰高調波発振器の運動は、次の伝達関数で記述されます。

$$H(z)=\frac{N_{Sensor}(z)}{1-2z^{-1}\cos \Delta t+z^{-2}}$$,

分子は、測定される振幅によって異なります。

• 変位: $$N_{dis}(z)=(z^{-1}+z^{-2})(1-\cos \Delta t)$$

• 速度: $$N_{vel}(z)=(z^{-1}-z^{-2})\sin \Delta t$$

• 加速度: $$N_{acc}(z)=(1-z^{-1})-(z^{-1}-z^{-2})\cos \Delta t$$

可能性のある 3 種類のシステム応答センサーの周波数応答関数を計算します。2 Hz のサンプル レートと 30,000 のホワイト ノイズのサンプルを入力として使用します。

fs = 2;
dt = 1/fs;
N = 30000;

u = randn(N,1);

ydis = filter((1-cos(dt))*[0 1 1],[1 -2*cos(dt) 1],u);
[frfd,fd] = modalfrf(u,ydis,fs,hann(N/2),Sensor="dis");

yvel = filter(sin(dt)*[0 1 -1],[1 -2*cos(dt) 1],u);
[frfv,fv] = modalfrf(u,yvel,fs,hann(N/2),Sensor="vel");

yacc = filter([1 -(1+cos(dt)) cos(dt)],[1 -2*cos(dt) 1],u);
[frfa,fa] = modalfrf(u,yacc,fs,hann(N/2),Sensor="acc");

loglog(fd,abs(frfd),fv,abs(frfv),fa,abs(frfa))
grid
legend(["dis" "vel" "acc"],Location="best")

いずれの場合も、生成される周波数応答関数は変位に対応する形式になっています。速度と加速度の測定値は、それぞれ変位測定値の 1 次微分と 2 次微分です。周波数応答関数は、システムの固有振動数付近の範囲では等価です。固有振動数から離れたところでは、周波数応答関数は異なります。

出力引数

frf — 周波数応答関数

ベクトル | 行列 | 3 次元配列

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

modalfrf は、センサーの種類に関係なく、周波数応答関数を常に動的柔軟性 (レセプタンス) 形式で出力します。

f — 周波数

ベクトル

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

coh — 多重コヒーレンス行列

行列

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

参照

[1] "Dynamic Stiffness, Compliance, Mobility, and more..." Siemens, last modified 2019, https://community.sw.siemens.com/s/article/dynamic-stiffness-compliance-mobility-and-more.

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

[3] Irvine, Tom. "An Introduction to Frequency Response Functions," Vibrationdata, 2000, https://vibrationdata.com/tutorials2/frf.pdf.

[4] 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 で導入

参考

modalfit | modalsd | n4sid (System Identification Toolbox) | tfestimate

トピック

• 識別したモデルのモード解析
• System Identification Overview (System Identification Toolbox)
• システム同定ワークフロー (System Identification Toolbox)
• Supported Continuous- and Discrete-Time Models (System Identification Toolbox)