【MATLAB】半値幅法の弱点を克服！高校数学でわかる最新のモード同定法（カーブフィット）

はじめに：半値幅法の課題と新しい手法の登場

これまでに、当サイトではモード円適合やプロニー法（ポリレファレンス法）、偏分反復法、直交多項式法（RFP法）などの高度なカーブフィット手法を解説してきました。

しかし、これらの手法は行列演算などの専門的な数学知識が必要であり、自分で一からプログラムを組むのはハードルが高いのが実情です。そのため、振動の専門家ではない現場の設計者たちの間では、直感的に理解しやすい半値幅法がいまだに根強く使われています。

ところが、半値幅法には以下のような致命的なデメリット（精度の低下）が存在します。

• 周波数分解能 $\Delta f$ が粗いと、正確なピーク値と半値点が読めない
• 減衰が非常に小さい場合、ピークが鋭すぎてデータ点が存在しない
• 低周波数帯に共振がある場合、近似式の前提が崩れて誤差が大きくなる

つまり、私のようなモード同定マニアではない一般のエンジニア向けに、「半値幅法のように簡単で、かつ半値幅法の弱点を克服した高精度な同定手法」が求められていました。

そんな中、2023年に中部大学から非常に実用的な新しいモード同定手法（リンクの109番）が発表されました。この手法は半値幅法よりも簡単に理解でき、複雑な行列演算も不要です。今回はこの画期的なモード同定手法の理論と実装について解説します。

 

 

新しいモード同定手法の理論（2023年発表）

本手法は1自由度法に分類されます。つまり、1つの共振ピークに対して、共振周波数とモード減衰比を同定するアプローチです。隣り合う共振が近接していない孤立モードにおいて非常に有効です（近接モードの意味がわからない方はこちらをご覧ください）。

$r$ 次の共振周波数近傍のコンプライアンス $G$（変位／力の周波数応答関数）は、式(1)で表せます。ここで、$\omega$ は角周波数、$m_r$、$c_r$、$k_r$ はそれぞれ $r$ 次のモード質量、モード減衰、モード剛性を意味し、$\phi_{ri}$ と $\phi_{rj}$ は加振点 $i$ および応答点 $j$ の $r$ 次モードベクトルの成分です。

$$G(\omega) = \frac{x_j}{F_i e^{j\omega t}} = \frac{\phi_{ri}\phi_{rj}}{-m_r\omega^2 + j c_r\omega + k_r} \tag{1}$$

この式を、剛性 $k_r$ と減衰比 $\zeta$、および角共振周波数 $\Omega$ を用いて無次元化してまとめると式(2)になります。

$$G(\omega) = \frac{\frac{\phi_{ri}\phi_{rj}}{k_r}}{1 – \beta^2 + 2j\zeta\beta}, \quad \beta=\frac{\omega}{\Omega} \tag{2}$$

一般的に、実験モード解析ではコンプライアンス（変位）ではなく、加速度センサーから得られるアクセレランス（加速度／力の周波数応答関数）を用いることが多いため、式(2)に $-\omega^2$ を掛けてアクセレランス $L(\omega)$ の式(3)に変換します。

$$L(\omega) = \frac{-\omega^2 \frac{\phi_{ri}\phi_{rj}}{k_r}}{1 – \beta^2 + 2j\zeta\beta}, \quad \beta=\frac{\omega}{\Omega} \tag{3}$$

式(3)を実部と虚部に分ける

式(3)の分母を実数化（複素共役を掛ける）して、実部と虚部に分解すると式(4)および式(5)となります。

$$Real[L(\omega)] = \frac{-\omega^2 \frac{\phi_{ri}\phi_{rj}}{k_r}(1 – \beta^2)}{(1 – \beta^2)^2 + (2\zeta\beta)^2} \tag{4}$$

$$Imag[L(\omega)] = \frac{\omega^2 \frac{\phi_{ri}\phi_{rj}}{k_r}(2\zeta\beta)}{(1 – \beta^2)^2 + (2\zeta\beta)^2} \tag{5}$$

ここで、式(5)を式(4)で割り、実部と虚部の比をとると、非常にシンプルな式(6)が得られます。（余計な定数や係数がすべて綺麗に相殺されるのがポイントです）

$$\frac{Imag[L(\omega)]}{Real[L(\omega)]} = \frac{-2\zeta\beta}{1 – \beta^2} \tag{6}$$

式(3)を振幅と位相で表現する

一方で、アクセレランス $L(\omega)$ を極座標形式（振幅と位相）で表現すると式(7)になります。

$$L(\omega) = |\omega^2 G(\omega)| e^{j\theta} \tag{7}$$

オイラーの公式を用いて実部と虚部に分けると式(8)、式(9)となります。

$$Real[L(\omega)] = |\omega^2 G(\omega)| \cos\theta \tag{8}$$
$$Imag[L(\omega)] = |\omega^2 G(\omega)| \sin\theta \tag{9}$$

同様に、実部と虚部の比をとると $\tan\theta$ が導かれます。

$$\frac{Imag[L(\omega)]}{Real[L(\omega)]} = \tan\theta \tag{10}$$

式(6)と式(10)から共振周波数とモード減衰比を導出

導出過程は異なりますが、式(6)と式(10)はどちらも「アクセレランスの虚部と実部の比」を表しているため等価です（式(6) = 式(10)）。

実験データには、アクセレランスの位相データ（実測結果）が存在します。そこで、$r$ 次の共振周波数近傍の適当な2つの角周波数 $\omega_1$ と $\omega_2$ を選ぶと、以下の連立方程式(11)、式(12)が成立します。

$$\frac{-2\zeta \frac{\omega_1}{\Omega}}{1 – \left(\frac{\omega_1}{\Omega}\right)^2} = \tan\theta_1 \tag{11}$$

$$\frac{-2\zeta \frac{\omega_2}{\Omega}}{1 – \left(\frac{\omega_2}{\Omega}\right)^2} = \tan\theta_2 \tag{12}$$

この式(11)および式(12)を、未知数である $\Omega$（角共振周波数）と $\zeta$（モード減衰比）について解くと、見事に式(13)と式(14)が導き出されます。

$$\Omega = \sqrt{\frac{\omega_1 \omega_2 (\omega_1 – \omega_2) \frac{\tan\theta_2}{\tan\theta_1}}{\omega_2 – \omega_1 \frac{\tan\theta_2}{\tan\theta_1}}} \tag{13}$$

$$\zeta = \frac{1}{2} \tan\theta_1 \left(\frac{\Omega^2 – \omega_1^2}{\Omega\omega_1}\right) \tag{14}$$

【この手法の凄さ】
これまで数多くのモード同定手法を勉強してきましたが、行列式や反復計算を使わず、高校の代数計算レベルの式展開だけで共振周波数とモード減衰比を一発で同定できるアルゴリズムは他にありません。しかも必要なデータは、共振付近のたった「2点の周波数と位相データ」だけです。いやー、本当にすごい発想ですね。

 

例題での精度検証

では、図1の4自由度モデルの周波数応答を例題として実際に解いてみましょう。
図1を見ると、1次共振周波数は明らかに 8Hz から 10Hz の間に存在しています。そこで、データ点として 8Hz と 10Hz の2点だけを抽出し、式(13)と式(14)に代入して1次の共振周波数とモード減衰比を計算してみます。

図1：周波数応答関数のグラフ

 

計算結果を表1に示します。
たった2点のデータを代入しただけで、非常に高精度に同定できていることがわかります。半値幅法のように「正確にピーク周波数を探す」必要もありません。

表1：同定結果の比較

	真値	同定結果	誤差 [%]
共振周波数 [Hz]	9.4321	9.5061	-0.7851
モード減衰比	0.01	0.0093	7.0806

 

MATLABプログラム（例題の実装コード）

上記の計算をMATLABで実行するコードです。実行すると以下のような図が出力され、同定された共振周波数に赤い破線が引かれます。

実行メインコード

clear all;clc;close all

freq=1:1:100; %計算周波数の定義
w=2*pi*freq; %角振動数

m_vec=ones(1,4);
k_vec=(1:4)*2*10^4;
[M]=eval_Mmatrix(m_vec); %質量行列
[K]=eval_Kmatrix(k_vec); %剛性行列

[V,D]=eig(K,M); %固有値D、固有ベクトルV
wn=sqrt(diag(D)); %固有角振動数
fn=wn/(2*pi); %固有振動数

mr=diag(V.'*M*V); %モード質量行列
kr=diag(V.'*K*V); %モード剛性行列
c_cr=2*sqrt(mr.*kr);
modal_dampimg=0.01; %モード減衰比（0.01と設定）
cr=c_cr*modal_dampimg; %モード減衰係数

F=zeros(length(m_vec),1); F(1)=1; %外力ベクトル(全周波数帯を1で加振する)
Xj=zeros(length(m_vec),length(freq)); %加速度ベクトル（計算前にベクトルを定義してメモリを確保）

% % % % モード法によるFRF合成
ii=1; %加振点
for ii1=1:length(m_vec) %応答点
    for ii2=1:length(wn)
        Xj(ii1,:)=Xj(ii1,:)+w.^2.*V(ii1,ii2)*V(ii,ii2)./( -mr(ii2).*w.^2 + 1i*cr(ii2)*w +kr(ii2) )*F(ii);
    end
end


% % % % % 手法の適用：ω1とω2の設定（8Hzと10Hz）
id1=8;
id2=10;
w1=w(id1);
w2=w(id2);

% % % % % tanθ1とtanθ2を算出
tan1=tan(angle(Xj(1,id1)));
tan2=tan(angle(Xj(1,id2)));

% % % % % 式(13)(14)による同定
wr=sqrt( (w1*w2*(w1-w2*tan2/tan1))/(w2-w1*tan2/tan1) ); %式(13) 角共振周波数
fn_re=wr/(2*pi);  %周波数に変換
zeta1=-1/2*tan1*((wr^2-w1^2)/(wr*w1));  %式(14) モード減衰比


% % % % % グラフ描画と結果出力
figure(1)
semilogy(freq,abs(Xj(1,:)),'k')
hold on
semilogy([fn_re fn_re],[min(abs(Xj(1,:))) max(abs(Xj(1,:)))],'r--')
xlabel('周波数 Hz')
ylabel('加速度/力 [(m/s^2)/N]')
message1=['共振周波数 ' num2str(fn_re) 'Hz   ,モード減衰比' num2str(zeta1)];
text(fn_re*0.5,max(abs(Xj(1,:)))*1.2,message1)

disp('【共振周波数の比較】')
disp(['真値: ', num2str(fn(1)), ' Hz'])
disp(['同定値: ', num2str(fn_re), ' Hz'])
disp(['誤差: ', num2str((fn(1)-fn_re)/fn(1)*100), ' %'])
disp(' ')
disp('【モード減衰比の比較】')
disp(['真値: ', num2str(modal_dampimg)])
disp(['同定値: ', num2str(zeta1)])
disp(['誤差: ', num2str((modal_dampimg-zeta1)/modal_dampimg*100), ' %'])

外部関数ファイル (eval_Mmatrix.m)

function [M]=eval_Mmatrix(m_vec)
    M=diag(m_vec);
end

外部関数ファイル (eval_Kmatrix.m)

function [K]=eval_Kmatrix(k_vec)
    K=zeros(length(k_vec));
    for ii1=1:length(k_vec)
        if ii1==1
            K(ii1,ii1)=k_vec(ii1);
        else
            K(ii1-1:ii1,ii1-1:ii1)=K(ii1-1:ii1,ii1-1:ii1)+[1 -1;-1 1]*k_vec(ii1);
        end
    end
end