- MATLABのSignal Processing Toolboxにある
modalfit関数の実用的な精度 - 3つのカーブフィット手法（’lsce’, ‘lsrf’, ‘pp’）の比較結果
- モード解析（同定）において最も精度が高くオススメの手法
はじめに：MATLABの「modalfit」関数とは？
過去の記事では、モード同定法（カーブフィット）について理論やMATLABプログラムを解説してきました。
例えば、半値幅法やモード円適合、プロニー法（ポリレファレンス法）、偏分反復法、RFP法などについて、バネマスモデルを用いて同定精度を検証しました。
これらは理論を理解する上では非常に有効ですが、読者の中には「自作のコードではなく、MATLAB公式のライブラリ関数を利用してサクッと解析したい」という方も多いでしょう。
そんな方にオススメなのが、Signal Processing Toolboxに収録されている関数 modalfit です。
この関数では、以下の3つのカーブフィット手法を指定してモーダルパラメータを推定できます。
- ‘lsce’：最小二乗複素指数法（時間領域）
- ‘lsrf’：最小二乗有理分数法（周波数領域）
- ‘pp’：ピークピッキング法（簡易手法）
しかし、カーブフィットマニアの間では「modalfitは精度があまり良くない」という噂も耳にします。私自身、これまで modalfit を本格的に使ったことがなかったため、今回バネマスモデルを用いてその推定精度を徹底検証してみました。
モード解析における検証内容とモデル設定
下図のような4自由度のバネマスモデルを想定します。今回はカーブフィットアルゴリズムの純粋な推定精度を検証したいため、後から「真値の共振周波数」を設定しやすい手法をとります。
具体的には、とりあえず適当な質量 $m$ と剛性 $k$ を決めてモードベクトルを算出しておきます。モードベクトルは共振周波数によらず形状が不変であり、モード質量は1、モード剛性は固有値（角共振周波数の二乗）となります。そのため、モードベクトルさえ先に計算しておけば、後から任意の共振周波数を設定しても矛盾なく物理モデルを構成できます。
今回は結果をわかりやすくするため、真値の共振周波数を 10Hz, 100Hz, 1000Hz, 1500Hz とキリの良い数値に設定しました。
MATLABによる初期設定コードは以下の通りです。
N = 4; % 共振周波数の数
m_vec=ones(1,N);
k_vec=ones(1,N)*10^7;
[M]=eval_Mmatrix(m_vec); %質量行列
[K]=eval_Kmatrix(k_vec); %剛性行列
[V,D]=eig(K,M); %固有値D、固有ベクトルV
fn=[10 100 1000 1500]; %後から共振周波数を設定する
wn=2*pi*fn;
kr=wn.^2; % モード剛性
mr=ones(1,N);% モード質量
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
function [M]=eval_Mmatrix(m_vec)
M=diag(m_vec);
end
まずはモード減衰比を0.5%に設定し、共振周波数とモード減衰比に微小なランダムノイズ（0.5%以下）を加えた際のFRF（1加振1応答）を用いて、各手法の推定精度がどうなるかを検証します。
続いて、モード減衰比を段階的に大きくし（それに比例してノイズも増加）、精度がどのように推移するかを確かめます。
検証結果1：モード減衰比 約0.5%の場合
さっそく結果を見てみましょう。下図は横軸にモード次数、縦軸に真値と modalfit で推定した値との誤差（%）を示しています。
比較対象として、過去の記事で私が自作したRFP法（直交多項式法）の結果もプロットしています。
一目瞭然ですが、手法 ‘lsce’（最小二乗複素指数法）の誤差が異常に大きくなっています。 3次以上のモードではNaN（計算不能）となってしまい、グラフ化すらできませんでした。’lsce’ の誤差が大きすぎるため、他の手法の差が全く見えません。
そこで、’lsce’ を除外してグラフを再描画したものが下図です。
手法 ‘lsrf’（最小二乗有理分数法）と ‘pp’（ピークピッキング法）は、実用上十分な精度を保っています。
- 手法 ‘lsrf’：共振周波数の誤差 最大0.0003%、減衰比の誤差 0.2%
- 手法 ‘pp’：共振周波数の誤差 最大0.06%、減衰比の誤差 5.7%
- 自作RFP法：最も誤差が小さく、最高精度を記録
なお、MathWorksの公式ドキュメントによると、手法 ‘pp’ は入力するFRFの数（複数応答点など）を増やすことで精度が改善する可能性があります。今回は検証工数の都合上、1加振1応答のデータのみで検証しています。
検証結果2：モード減衰比 約0.5～10%に変化させた場合
次に、構造物の減衰がより大きい（ピークがなだらかになる）ケースを想定し、モード減衰比を変化させた場合の共振周波数の推定誤差を下図に示します。色がモード減衰比の違いを表しています。
やはりここでも ‘lsce’ は誤差が大きく不安定です。一方、’lsrf’ と ‘pp’ はモード減衰比が大きくなるにつれて推定誤差が増加する傾向にありますが、減衰比10%の過酷な条件下でも ‘lsrf’ の共振周波数誤差は0.8%程度に収まっており、実用上は全く問題ないレベルです。
（そしてここでも、自作のRFP法が最も安定した精度を叩き出しています）
続いて、モード減衰比自体の推定精度を下図に示します。
共振周波数以上に、減衰比の推定において手法間の差が顕著に出ました。
‘lsce’ に加え、簡易手法である ‘pp’ でも減衰比の推定精度が大きく悪化しています。近接モードや減衰が大きい系ではピークピッキングの限界が露呈しています。
一方で 手法 ‘lsrf’ は、モード減衰比の推定誤差が最大でも9.5%程度に収まりました。
真値が減衰比10%のとき、誤差が10%含まれたとしても推定値は「9%〜11%」の範囲に収まります。工学的な実用範囲を考えれば、十分信頼できる値と言えるでしょう。
結論と考察：modalfitでオススメの手法は？
今回の検証結果と、各手法のアルゴリズム的背景を踏まえると、以下の結論が得られます。
MATLABの
modalfit 関数を使用する場合は、迷わず 'FitMethod', 'lsrf'（最小二乗有理分数法） を指定するのが無難かつ高精度です。考察（なぜ差が出たのか？）：
- lsce（時間領域手法）：FRF（周波数領域）データからインパルス応答（時間領域）への逆変換を伴うため、計算過程でノイズや打ち切り誤差の影響を強く受け、不安定になったと考えられます。
- pp（ピークピッキング）：シンプルで速いですが、減衰が大きい（ピークが鈍い）場合やモードが近接していると、他モードの影響（裾野の被り）を分離できず精度が低下します。
- lsrf（周波数領域手法）：FRFデータを直接、有理分数多項式でカーブフィットするため、今回のような周波数応答データの解析と最も相性が良く、ノイズに対してもロバスト（頑健）に機能しました。
高精度なモード同定を行いたい方へ（無料ツールのご案内）
今回の検証で、MATLAB標準の modalfit (‘lsrf’) は実用に耐えうる精度を持つことがわかりました。
しかしグラフからもわかる通り、当サイトで解説してきたRFP法（直交多項式を用いた手法）が、ノイズ環境下でも最も推定誤差が小さく、非常に優秀なアルゴリズムであることも改めて証明されました。
「MATLABのライセンスを持っていない」
「より高精度かつ直感的にモード同定（カーブフィット）を行いたい」
「複雑なコードを書かずにGUIで解析したい」
そんな方に向けて、私自身が開発した無料のモード解析ツール「PolyResi」を公開しています。今回最高精度を叩き出したRFP法やPolyMAX法をベースにしており、クリック操作だけで高精度なモーダルパラメータの抽出とアニメーション描画が可能です。
ぜひ一度ダウンロードして、ご自身のデータで精度を体感してみてください！
MATLAB 実行コード（全文）
ご自身で検証を再現したい方向けに、今回使用したMATLABコードの全文を掲載しておきます。
clear all;clc;close all
N = 4; % 共振周波数の数
freq=1:1:2000; %計算周波数の定義
df=freq(2) - freq(1);
w=2*pi*freq; %角振動数
m_vec=ones(1,N);
k_vec=ones(1,N)*10^7;
[M]=eval_Mmatrix(m_vec); %質量行列
[K]=eval_Kmatrix(k_vec); %剛性行列
[V,D]=eig(K,M); %固有値D、固有ベクトルV
for ii1=1:length(m_vec)
V(:,ii1)=V(:,ii1)/V(1,ii1); % 加振点m1のモードベクトルの成分が1になるように正規化
end
% モードベクトルVは上の結果を使って、検証しやすいようにモード質量とモード剛性は再設定
noise_coef=[0.005 0.01:0.01:0.1];
mr=ones(1,N);
F=zeros(length(m_vec),1); F(1)=1; %外力ベクトル(全周波数帯を1で加振する)
N2=100;
modalfit_check.lsce.fn=zeros(N2,length(noise_coef),N);
modalfit_check.lsce.md=zeros(N2,length(noise_coef),N);
modalfit_check.lsrf.fn=zeros(N2,length(noise_coef),N);
modalfit_check.lsrf.md=zeros(N2,length(noise_coef),N);
modalfit_check.pp.fn=zeros(N2,length(noise_coef),N);
modalfit_check.pp.md=zeros(N2,length(noise_coef),N);
RFPM.fn=zeros(N2,length(noise_coef),N);
RFPM.md=zeros(N2,length(noise_coef),N);
for iin=1:length(noise_coef)
for iin2=1:N2
noise = rand(1,N)*noise_coef(iin) + ones(1,N) ;
fn=[10 100 1000 1500].*noise;
wn=2*pi*fn;
kr=wn.^2;
c_cr=2*sqrt(mr.*kr);
modal_dampimg=noise_coef(iin)*noise(1); %モード減衰比
cr=c_cr.*modal_dampimg; %モード減衰係数
Xj=zeros(1,length(freq)); %変位ベクトル（計算前にベクトルを定義してメモリを確保）
% モード法
ii=1; %加振点
for ii1=1 %応答点
for ii2=1:length(wn)
Xj(ii1,:)=Xj(ii1,:)+V(ii1,ii2)*V(ii,ii2)./( -mr(ii2).*w.^2 + 1i*cr(ii2)*w +kr(ii2) )*F(ii);
end
end
FRF_E=Xj(1,:);
ideal_md=modal_dampimg; ideal_fn=fn;
[fn2,dr2] = modalfit(FRF_E.',freq.',freq(end)*2,4,'FitMethod','lsce');
modalfit_check.lsce.md(iin2,iin,:)=abs(ideal_md - dr2)/ideal_md;
modalfit_check.lsce.fn(iin2,iin,:)=abs(ideal_fn - fn2.')./ideal_fn;
[fn2,dr2] = modalfit(FRF_E.',freq.',freq(end)*2,4,'FitMethod','lsrf');
modalfit_check.lsrf.md(iin2,iin,:)=abs(ideal_md - dr2)/ideal_md;
modalfit_check.lsrf.fn(iin2,iin,:)=abs(ideal_fn - fn2.')./ideal_fn;
[fn2,dr2] = modalfit(FRF_E.',freq.',freq(end)*2,4,'FitMethod','pp');
modalfit_check.pp.md(iin2,iin,:)=abs(ideal_md - dr2)/ideal_md;
modalfit_check.pp.fn(iin2,iin,:)=abs(ideal_fn - fn2.')./ideal_fn;
[FRF_recalcu,modal_parameter]=Rational_Fraction_Polynomial_Method(FRF_E,freq,4);
RFPM.md(iin2,iin,:)=abs(ideal_md - modal_parameter(:,2))/ideal_md;
RFPM.fn(iin2,iin,:)=abs(ideal_fn - modal_parameter(:,1).'/(2*pi))./ideal_fn;
end
end
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
function [M]=eval_Mmatrix(m_vec)
M=diag(m_vec);
end
function [FRF_recalcu,modal_parameter]=Rational_Fraction_Polynomial_Method(FRFdata,freq,N)
% modal_parameter = Modal Parameters [freq,damp,Ci,Oi]:
w=2*pi*freq;
[r,c]=size(w);
if r<c, w=w.'; end
[r,c]=size(FRFdata);
if r<c, FRFdata=FRFdata.'; end
nom_w=max(w);
w=w./nom_w;
m=2*N-1;
n=2*N;
[phi,coeff_A]=Complex_Orthogonal_Polynomials(FRFdata,w,1,m);
[theta,coeff_B]=Complex_Orthogonal_Polynomials(FRFdata,w,2,n);
[r,c]=size(phi);
P=phi(:,1:c);
[r,c]=size(theta);
Theta=theta(:,1:c);
T=sparse(diag(FRFdata))*theta(:,1:c-1);
W=FRFdata.*theta(:,c);
X=-2*real(P'*T);
G=2*real(P'*W);
D=-inv(eye(size(X))-X.'*X)*X.'*G;
C=G-X*D;
D=[D;1];
FRF_recalcu=zeros(length(w),1);
for n=1:length(w)
numer=sum(C.'.*P(n,:));
denom=sum(D.'.*Theta(n,:));
FRF_recalcu(n)=numer/denom;
end
A=coeff_A*C;
[r,c]=size(A);
A=A(r:-1:1).';
B=coeff_B*D;
[r,c]=size(B);
B=B(r:-1:1).';
[R,P,K]=residue(A,B);
[r,c]=size(R);
for n=1:(r/2)
Residuos(n,1)=R(2*n-1);
Polos(n,1)=P(2*n-1);
end
[r,c]=size(Residuos);
Residuos=Residuos(r:-1:1)*nom_w;
Polos=Polos(r:-1:1)*nom_w;
fn=abs(Polos);
damp=-real(Polos)./abs(Polos);
Ai=-2*(real(Residuos).*real(Polos)+imag(Residuos).*imag(Polos));
Bi=2*real(Residuos);
const_modal=complex(Ai,abs(Polos).*Bi);
Ci=abs(const_modal);
Oi=angle(const_modal).*(180/pi);
modal_parameter=[fn, damp, Ci, Oi];
end
function [P,coeff]=Complex_Orthogonal_Polynomials(FRFdata,w,phitheta,kmax)
if phitheta==1
q=ones(size(w));
elseif phitheta==2
q=(abs(FRFdata)).^2;
end
R_m1=zeros(size(w));
R_0=1/sqrt(2*sum(q)).*ones(size(w));
Rminus=[R_m1,R_0];
coeff=zeros(kmax+1,kmax+2);
coeff(1,2)=1/sqrt(2*sum(q));
Rplus=[];
R=[Rminus Rplus];
for k=1:kmax
Vkm1=2*sum(w.*R(:,k+1).*R(:,k).*q);
Sk=w.*R(:,k+1)-Vkm1*R(:,k);
Dk=sqrt(2*sum((Sk.^2).*q));
Rplus=[Rplus Sk/Dk];
coeff(:,k+2)=-Vkm1*coeff(:,k);
coeff(2:k+1,k+2)=coeff(2:k+1,k+2)+coeff(1:k,k+1);
coeff(:,k+2)=coeff(:,k+2)/Dk;
R=[Rminus Rplus];
end
R=R(:,2:kmax+2);
coeff=coeff(:,2:kmax+2);
P=R.*1i.^(0:kmax);
jk=1i.^(0:kmax);
coeff=(jk'*jk).*coeff;
end
コメント