MATLABで学ぶ振動工学・信号処理 【完全版】FRF(周波数応答関数)のプログラム

はじめに

伝達関数(周波数応答関数:Frequency Response Function(FRF))には、H1推定、H2推定、H3推定、Hv推定というようにいくつか種類があります。先に挙げた4手法については過去に紹介をしました。これらの記事をまとめつつ、プログラムのコーディングについて語りたいと思います。

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

 

スポンサーリンク

プログラムのコーディングについて(余談なので読み飛ばして良いです)

最近の若者のコードを見てて、人生のパイセンとして一言!
コードが汚ない!見にくい!

実行ファイルなら別に汚くても良いんですよ。
ただ、functionファイルは綺麗に書きましょうよ!

functionファイルは資産ですよ。

 

ということで、パイセンからfunctionファイルを作成するときの頭の中をイメージを言語化しながら語らせていただきます。

私の考えるfunctionファイル化するときの重要なことは下記です。
①functionファイル単体で機能する
②functionファイルの中に別のfunctionファイル
③数年後に見返してもわかる(教科書と同じ文字式にするなど、処理が遅くなっても理解しやすいように書く)

②に関しては、今回のFRFでは下記のようなイメージです。

どこの矢印でfunctionファイルのつながりが切れても、その位置からのfunctionファイルとして使用しても機能が損なわれないようにしています。
例えば、H1.mでも機能しますし、PowerSpectralDensity.mでも機能するように、機能を詰めすぎずに、整理することが重要です。

先に上図のような関係図を頭の中にイメージして、functionファイル化するのが大事だと思います。

③に関しては、処理速度を重視する人とは意見が割れると思います。
私の考えでは、処理速度よりも見やすさが重要です。
処理速度重視のときは、その時用にコードを少し書き換えれば良いだけ!との認識でいます。

それよりも、数か月後、数年後に自分で作ったコードを見返して、理解できないときにロスする時間が膨大です。

みなさんエンジニアなので、毎年のように論文や参考書を読んだり、製品設計などの業務をこなします。こんな日々の中で、過去に作ったコードの中身が汚いから解読できない!ってことを経験した人は多いはずです。

また、他人とコードを共有したいけど、自分のコードが汚いので、説明しないとダメor作り直さないとダメ!みたいな経験をした人もいるかと思います。

私はこのようなケースをできる限り回避できるようにコードを書くことが最も効率が良いと考えています。

なので、学生の方には今回のコードを参考にしていただけると嬉しいです。
(コードは最後に添付)

H1,H2,H3,Hvの理論式

ここでは説明はせずに、式だけ紹介します。

$$
H1(ω)= \frac{\sum_{n=1}^{N} W_{fx}(ω)}{\sum_{n=1}^{N} W_{ff}(ω)}
$$

$$
H2(ω)= \frac{\sum_{n=1}^{N} W_{xx}(ω)}{\sum_{n=1}^{N} W_{xf}(ω)}
$$

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

$$
Hv(ω)= \sqrt{H1(ω)} \sqrt{H2(ω)}
$$

スポンサーリンク

MATLABでの比較検証

①入力にノイズが含まれているケース(実行ファイル①)

おそらく、H2推定が最も良い結果になると思いますよねー。さてどうなるか確認してみましょう。
現実的な問題として考えると、入力にノイズが入るとは考えにくいけども、入ったとしても1%以下だと思います。
1%の誤差(ノイズ)の場合、下図の結果となりました。


図1 誤差1%の場合

ほぼ差がないですね。差がないと考察の使用がないので、とりあえず、誤差300%を与えてみます。

図2 誤差300%の場合

理想的なFRFとの差が小さい順番は下記になります。
H2推定>H3推定>Hv推定>H1推定

②応答にノイズが含まれているケース(実行ファイル②)

応答にノイズが入ることは良くありますよね。
今回はとりあえず、最大振幅X(ω)を基準として3%の誤差を与えてみました。結果は下図です。

図3 誤差3%の場合

理想的なFRFとの差が小さい順番は下記になります。
H1推定>Hv推定>H3推定>H2推定

③入力と応答にノイズが含まれているケース(実行ファイル③)

現実的ではないですが、入力に300%のノイズ、応答に3%のノイズが混入した場合の結果は下図です。

図4 入力に300%の誤差、応答に3%の誤差

理想的なFRFとの差が小さい順番は下記になります。
Hv推定>H3推定>H2推定>H1推定
(共振周波数の値と、それ以外の周波数でピークがないことを重要視しました)

 

通常の実験環境を考えると、ノイズが混入したとしても、入力に3%程度、応答に5%程度なのかなーと想像し、この条件での結果を比較してみます。

図5 入力に3%の誤差、応答に5%の誤差

理想的なFRFとの差が小さい順番は下記になります。
H1推定>Hv推定>H3推定>H2推定
(共振周波数の値と、それ以外の周波数でピークがないことを重要視しました)

FRFを測定できるFFTアナライザでは、デフォルトの設定でFRFの推定手法が”H1″となっていることが多いです。図3や図5の結果のように、通常の測定環境では”H1″推定が最も推定精度が高いと考えられます。

図3や図5以外の条件(ノイズ環境)では、各々で推定手法を選択する必要があります。

スポンサーリンク

コード

実行ファイル①
clear 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,:);
    ref(ii1,:)=F(1,:)+( rand(1,length(freq)) - 0.5 )*2*3; % Fを基準として±300%のノイズを加える
end

[h1] = FRF(ref,res,freq,'H1');
[h2] = FRF(ref,res,freq,'H2');
[h3] = FRF(ref,res,freq,'H3');
[hv] = FRF(ref,res,freq,'Hv');
[r2] = Coherence(ref,res,freq);

figure(1)
loglog(freq,abs(FRF_ideal_x3_f1),'k','linewidth',7)
hold on
loglog(freq,abs(h1),'r--','linewidth',3)
loglog(freq,abs(h2),'b--','linewidth',3)
loglog(freq,abs(h3),'g--','linewidth',3)
loglog(freq,abs(hv),'m.','linewidth',3)
hold off
grid on
legend('理想的なFRF','H1推定','H2推定','H3推定','Hv推定')
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('FRF m/N');

figure
semilogx(freq,r2,'k','linewidth',3)
xlabel('周波数 Hz');ylabel('Coherence');
grid on
ylim([0 1])
xlim([freq(1) freq(end)])

実行ファイル②
clear 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,:);
end

[h1] = FRF(ref,res,freq,'H1');
[h2] = FRF(ref,res,freq,'H2');
[h3] = FRF(ref,res,freq,'H3');
[hv] = FRF(ref,res,freq,'Hv');
[r2] = Coherence(ref,res,freq);

figure(1)
loglog(freq,abs(FRF_ideal_x3_f1),'k','linewidth',7)
hold on
loglog(freq,abs(h1),'r--','linewidth',3)
loglog(freq,abs(h2),'b--','linewidth',3)
loglog(freq,abs(h3),'g--','linewidth',3)
loglog(freq,abs(hv),'m.','linewidth',3)
hold off
grid on
legend('理想的なFRF','H1推定','H2推定','H3推定','Hv推定')
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('FRF m/N');

figure
semilogx(freq,r2,'k','linewidth',3)
xlabel('周波数 Hz');ylabel('Coherence');
grid on
ylim([0 1])
xlim([freq(1) freq(end)])
実行ファイル③
clear 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.05; % 最大変位Xを基準として5%のノイズを加える
    ref(ii1,:)=F(1,:)+( rand(1,length(freq)) - 0.5 )*2*0.03; % Fを基準として±3%のノイズを加える
end

[h1] = FRF(ref,res,freq,'H1');
[h2] = FRF(ref,res,freq,'H2');
[h3] = FRF(ref,res,freq,'H3');
[hv] = FRF(ref,res,freq,'Hv');
[r2] = Coherence(ref,res,freq);

figure(1)
loglog(freq,abs(FRF_ideal_x3_f1),'k','linewidth',7)
hold on
loglog(freq,abs(h1),'r--','linewidth',3)
loglog(freq,abs(h2),'b--','linewidth',3)
loglog(freq,abs(h3),'g--','linewidth',3)
loglog(freq,abs(hv),'m.','linewidth',3)
hold off
grid on
legend('理想的なFRF','H1推定','H2推定','H3推定','Hv推定')
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('FRF m/N');

figure
semilogx(freq,r2,'k','linewidth',3)
xlabel('周波数 Hz');ylabel('Coherence');
grid on
ylim([0 1])
xlim([freq(1) freq(end)])

 

functionファイル
function [frf] = FRF(ref,res,freq,FRFtype)


switch FRFtype
    case 'H1'
        [frf] = H1(ref,res,freq);
    case 'H2'
        [frf] = H2(ref,res,freq);
    case 'H3'
        [frf] = H3(ref,res,freq);
    case 'Hv'
        [frf] = Hv(ref,res,freq);
end

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
function [r2] = Coherence(ref,res,freq)

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

r2=h1./h2;

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 [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

コメント