MATLABで学ぶ振動工学　FRFのH3推定およびHv推定について

2021.06.02

はじめに

伝達関数（周波数応答関数：Frequency Response Function(FRF)）には、H1推定、H2推定、H3推定、Hv推定というようにいくつか種類があります。
これまではH1推定とH2推定を紹介しました。本記事ではH3推定とHv推定について説明します。

MATLABで学ぶ振動工学　FRFのH1推定について（伝達関数のノイズ低減）
伝達関数（周波数応答関数：Frequency Response Function(FRF)）には、H1推定、H2推定、H3推定、Hv推定というようにいくつか種類があります。 本記事ではH1推定について説明します。 MATLABのコードも載せています。
mech-eng-uni.com
2021.05.28
MATLABで学ぶ振動工学　FRFのH2推定について
伝達関数（周波数応答関数：Frequency Response Function(FRF)）には、H1推定、H2推定、H3推定、Hv推定というようにいくつか種類があります。 本記事ではH2推定について説明します。また、コヒーレンス関数（Coherebce Function）についても少しだけ説明しています。 MATLABのコードも載せています。
mech-eng-uni.com
2021.05.31

次回はコヒーレンス関数について順次説明していこうと思います。

H3推定とHv推定について

これまでに、入力にノイズ成分が含まれるときはH2推定を使用するのが望ましく、
応答にノイズ成分が含まれるときはH1推定を使用することが望ましいと説明しました。

H3推定やHv推定は、入力＆応答にノイズが含まれるときに使用する推定方法です。

ちなみに、私の愛読している長松先生著の”モード解析入門”には下記のような説明がありました。

入力と出力の両方に誤差が混入するときの周波数応答関数は、近似的に次式で定義したH3、厳密に定義したHvなどがある。

H3は近似的で、Hvは厳密なんですね。
ちなみに”モード解析入門”にはHv推定の理論式は紹介されてません。
（長松先生ーーなんでー！加振機とか打撃試験の内容じゃなくて、理論書いてよーって思ってます）

スポンサーリンク

理論式

H3推定

教科書に書いてあるH2推定の式自体はめちゃめちゃシンプルです。

$$
H3(ω)= \frac{H1(ω)+H2(ω)}{2}
$$

めっちゃシンプルですね。説明するところないですね。

Hv推定

Hv推定の理論式は下記です。
$$
Hv(ω)= \sqrt{H1(ω)} \sqrt{H2(ω)}
$$

こちらもめっちゃシンプルですね。

スポンサーリンク

MATLABでの検証

解析結果

まずは解析結果を先に示します。（100回平均）

ノイズの成分のない理想的なFRFを黒線、
ノイズ成分を含み平均化処理をしていないFRFを青線、
（平均化処理をした）H3推定結果を緑点線、
（平均化処理をした）Hv推定結果を赤点線で表しています。

どうですかね？
理想的なFRFと比較すると、応答の小さい10Hz以上でH3推定もHv推定も誤差が大きいですね。
長松先生が仰っていたように、Hv推定の方がH3推定よりも精度が高そうですね。

 

コード

実行ファイル
cleaclear all
close all
% % % % 初期設定
freq=0.5:0.1:20;
m_vec=[1 1 1 1];
[M]=eval_Mmatrix(m_vec);
k_vec=[1 1 1 1]*10^3;
[K]=eval_Kmatrix(k_vec);
C=zeros(length(k_vec));
F=zeros(length(m_vec),length(freq));
F(1,:)=ones(1,length(freq));

% % % % 変位X（ノイズなし）
X=eval_direct_x_2ndedition(M,K,C,F,freq);
% % % % 理想的なFRF
FRF_ideal_x3_f1=X(3,:)./F(1,:);

N=100; %平均化回数 100回
res=zeros(N,length(freq));
ref=zeros(N,length(freq));
for ii1=1:N
    res(ii1,:)=X(3,:)+( rand(1,length(freq)) - 0.5 )*max(abs(X(3,:)))*0.03; % 最大変位Xを基準として3%のノイズを加える
    ref(ii1,:)=F(1,:)+( rand(1,length(freq)) - 0.5 )*2*0.10; % Fを基準として±10%のノイズを加える
end

[h3] = H3(ref,res,freq);
[hv] = Hv(ref,res,freq);

figure(1)
loglog(freq,abs(FRF_ideal_x3_f1),'k','linewidth',7)
hold on
loglog(freq,abs(res(1,:)./ref(1,:)),'b','linewidth',3)
loglog(freq,abs(h3),'g--','linewidth',3)
loglog(freq,abs(hv),'r--','linewidth',3)
hold off
grid on
legend('理想的なFRF','平均化処理をしないFRF','H3推定','Hv推定')
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('FRF m/N');

スポンサーリンク

functionファイル
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 x=eval_direct_x_2ndedition(M,K,C,F,freq)
% % 2021.05.28 Fを修正 
x=zeros(size(F,1),length(freq));
j=sqrt(-1);
for ii1=1:1:length(freq)
    w=2*pi*freq(ii1);
    x(:,ii1)=inv(K+j*w*C-w^2*M)*F(:,ii1);
end

end

function [Wff]=PowerSpectralDensity(F,freq)

if size(F,1)==length(freq)
    F=F.';
end

Wff=conj(F).*F;


end

function [Wfx]=CrossSpectralDensity(F,X,freq)

if size(F,1)==length(freq)
    F=F.';
end
if size(X,1)==length(freq)
    X=X.';
end



Wfx=conj(F).*X;


end

function [h1] = H1(ref,res,freq)

% h1=zeros(1,freq);

[Wfx]=CrossSpectralDensity(ref,res,freq);
[Wff]=PowerSpectralDensity(ref,freq);

h1=mean(Wfx,1)./mean(Wff,1);

end


function [h2] = H2(ref,res,freq)

% h1=zeros(1,freq);

[Wxf]=CrossSpectralDensity(res,ref,freq);
[Wxx]=PowerSpectralDensity(res,freq);

h2=mean(Wxx,1)./mean(Wxf,1);

end

function [h3] = H3(ref,res,freq)

[h1] = H1(ref,res,freq);
[h2] = H2(ref,res,freq);

h3 = (h1 + h2)/2;
end

function [hv] = Hv(ref,res,freq)

[h1] = H1(ref,res,freq);
[h2] = H2(ref,res,freq);

hv = sqrt(h1).*sqrt(h2);
end

コメント