MATLABで学ぶ振動工学コヒーレンス関数について - MATLABパイセンが教える振動・騒音・音響・機械工学

はじめに
コヒーレンス関数
MATLABでの検証
1. 関連

はじめに

伝達関数（周波数応答関数：Frequency Response Function(FRF)）には、H1推定、H2推定、H3推定、Hv推定というように数種類あることを過去に紹介しました。

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

伝達関数（周波数応答関数：Frequency Response Function(FRF)）には、H1推定、H2推定、H3推定、Hv推定があります。これらの理論式を紹介し、それぞれの推定手法の比較検証を行いました。MATLABのコードも載せています。

本記事ではFRFのコヒーレンス関数について説明します。

コヒーレンス関数は、”コヒーレンス”や”関連度関数”とも呼ばれます。
ざっくり説明すると、FRFがどの程度信用できるのかを示す指標です。

コヒーレンス関数

振動工学を学習している・していた方の中にはコヒーレンス関数をプログラミングしようとして、うまくいかなかったという人がいるかと思います。

コヒーレンス関数のややこしい点は2点あります。

ややこしい点 1つ目

ややこしい点の1つ目が、教科書や参考書ではコヒーレンス関数の数式が多数あることです。（ただし、展開していくと同じになる？はず？）

例えば、”モード解析入門”の中でコヒーレンス関数は5つの数式で表現されています。（下図は参考）

コヒーレンス関数の詳細情報は下記を参照ください。

<br />

$γ^2_{fx}= \frac{W_{vv}}{W_{xx}}=\frac{\bar{V}V}{W_{xx}}=\frac{\bar{F H1}F H1}{W_{xx}}=\frac{\bar{H1}FF H1}{W_{xx}}=\frac{\bar{H1}W_{ff} H1}{W_{xx}}$

$γ^2_{fx}= \frac{(\bar{W_{fx}} / \bar{W_{ff}})W_{ff}(W_{fx} / W_{ff})}{W_{xx}}=\frac{\bar{W_{fx}} W_{ff}W_{fx}}{W_{xx}\bar{W_{ff}} W_{ff}}$

$γ^2_{fx}= \frac{\bar{W_{fx}}W_{fx}}{W_{ff}W_{xx}}= \frac{<W_{fx}W_{fx}>}{W_{ff}W_{xx}}$

$γ^2_{fx}= \frac{\|W_{fx}\|^2}{\|F\|^2\|X\|^2}$

$γ^2_{fx}(ω)= \frac{H1}{H2}$

ややこしい点 2つ目

ややこしい点の2つ目が、どこで平均化するかが不明確であることです。

例えば、下の式ではWfxなどを単体で平均化するのか、分母と分子ごとに平均化するのか、下記の等式を平均化回数分求めてその後平均化処理するのか、不明確ですよね。

$γ^2_{fx}=\frac{\bar{W_{fx}} W_{ff}W_{fx}}{W_{xx}\bar{W_{ff}} W_{ff}}$

オススメの方法

下式ではH1とH2さえしっかり求めておけばコヒーレンス関数は簡単に求まります。

$γ^2_{fx}(ω)= \frac{H1}{H2}$

H1推定やH2推定がわかっていない方は下記を参考にしてください。

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での検証

解析結果

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

誤差がないとコヒーレンス関数は1となり、誤差が増加するにしたがってコヒーレンス関数は0に近づきます。

また、共振周波数では応答が高いので、SN比（信号とノイズの比）が改善されるので、コヒーレンス関数は1に近くなることがわかります。

コード

実行ファイル

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));
res2=zeros(N,length(freq));
res3=zeros(N,length(freq));
res4=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.01; % 最大変位Xを基準として1%のノイズを加える
    res2(ii1,:)=X(3,:)+( rand(1,length(freq)) - 0.5 )*max(abs(X(3,:)))*0.05;
    res3(ii1,:)=X(3,:)+( rand(1,length(freq)) - 0.5 )*max(abs(X(3,:)))*0.20;
    res4(ii1,:)=X(3,:)+( rand(1,length(freq)) - 0.5 )*max(abs(X(3,:)))*0.50;
    ref(ii1,:)=F(1,:);
end
 
[h1_1] = FRF(ref,res,freq,'H1');
[h1_2] = FRF(ref,res2,freq,'H1');
[h1_3] = FRF(ref,res3,freq,'H1');
[h1_4] = FRF(ref,res4,freq,'H1');
[r2_1] = Coherence(ref,res,freq);
[r2_2] = Coherence(ref,res2,freq);
[r2_3] = Coherence(ref,res3,freq);
[r2_4] = Coherence(ref,res4,freq);
[r2_0] = Coherence(ref(1,:),X(3,:),freq);
 
figure(1)
subplot(211)
loglog(freq,abs(FRF_ideal_x3_f1),'k','linewidth',7)
hold on
loglog(freq,abs(h1_1),'r--','linewidth',3)
loglog(freq,abs(h1_2),'g--','linewidth',3)
loglog(freq,abs(h1_3),'b--','linewidth',3)
loglog(freq,abs(h1_4),'c--','linewidth',3)
hold off
grid on
legend('誤差0%','誤差1%','誤差5%','誤差20%','誤差50%')
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('FRF m/N');
 
subplot(212)
semilogx(freq,r2_0,'k','linewidth',3)
hold on
semilogx(freq,r2_1,'r','linewidth',3)
semilogx(freq,r2_2,'g','linewidth',3)
semilogx(freq,r2_3,'b','linewidth',3)
semilogx(freq,r2_4,'c','linewidth',3)
 
xlabel('周波数 Hz');ylabel('Coherence');
grid on
ylim([0 1.05])
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