はじめに
伝達関数（周波数応答関数：Frequency Response Function(FRF)）には、H1推定、H2推定、H3推定、Hv推定というようにいくつか種類があります。先に挙げた4手法については過去に紹介をしました。これらの記事をまとめつつ、プログラムのコーディングについて語りたいと思います。
スポンサーリンク
プログラムのコーディングについて（余談なので読み飛ばして良いです）
最近の若者のコードを見てて、人生のパイセンとして一言！
コードが汚ない！見にくい！
実行ファイルなら別に汚くても良いんですよ。
ただ、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
コメント