MATLABで学ぶ振動工学 FRFのH1推定について(伝達関数のノイズ低減)

はじめに

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

伝達関数は(応答)/(入力)で求められるから、1種類じゃないの?

上記の理解は間違ってはいないです。
ただし、実測データにはノイズ成分が含まれています。この影響を低減するために平均化処理をします。FRFを算出時の平均化処理の仕方によって推定の名称が変わります。

また、ノイズの影響以外にも考えないといけないことがあります。
それは、実測では1回目の(応答)/(入力)と2回目の(応答)/(入力)に差異があることです。
(解析だと1回目も2回目も一緒になります)

なぜ、実測では差異が生じるのかについては、いろいろな原因があるので一概には言えませんが、代表的な例だと下記があげられます。
・ハンマリング試験でFRFを測定 → 加振力がばらつく → 加振力依存性の影響が顕在化
・加振試験でFRFを測定 → 時間とともになじむ(ねじ止めなどの影響) → 周波数特性が変化

原因はいろいろありますが、FRFは測定ごとに差異があるので、どれくらい信頼性があるのかを示す評価指標があり、それはコヒーレンス関数と呼ばれます。

本日はH1推定について説明します。
次回はH2推定について、次にH3推定、Hv推定、コヒーレンス関数について順次説明していこうと思います。

H1推定について

先にも述べましたが、FRFには、H1推定、H2推定、Hv推定というようにいくつか種類があります。その中で、H1推定は伝達関数を推定する方法で最も一般的な手法です。

 

まず、理想的なFRFを考えてみましょう。理想的なFRF G(ω)は下記のようになります。

ただ、現実問題として、このような理想通りの測定はできません。それはノイズの影響があるからです。

例えば、スピーカーを用いた音響測定であれば、入力信号がスピーカーの電圧で、応答がマイクロフォンで測定する音になります。しかし、音には暗騒音(ノイズ)が含まれてしまいますよね。
この例のように、入力信号はノイズの影響が小さくても、応答信号はノイズの影響があるということが多いです。(もしくは応答信号のノイズの影響が不明)

上記の例をブロック線図のように表すと下図になります。

上図のような応答信号にノイズ成分が含まれる状況で、FRF=(応答V)/(入力F)を推定する方法がH1推定です。

スポンサーリンク

H1推定の理論

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

$$
H1(ω)= \frac{W_{fx}(ω)}{W_{ff}(ω)}
$$

ここでWfxはFとXのクロススペクトル密度関数(Cross Spectral Density Function)、WffはFのパワースペクトル密度関数(Power Spectral Density Function)もしくはオートスペクトル密度関数(Auto Spectral Density Function)である。

ただ、教科書に書いてあるHz推定の式は少し間違っています。実際には下記になります。Nは平均化の回数です。

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

もしくは(下記で言いたいのは、クロススペクトル密度関数の平均とパワースペクトル密度関数の平均から求められるよってことです)

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

 

では、なぜ上式で応答に混入するノイズが低減できるのかを説明します。
なお、*は共役を指します。(*は×という意味ではないので注意)

$$
H1(ω)= \frac{\sum_{n=1}^{N} F(ω)^*X(ω)}{\sum_{n=1}^{N} F(ω)^*F(ω)}
$$
$$
H1(ω)= \frac{\sum_{n=1}^{N} F(ω)^*(V(ω)+N(ω))}{\sum_{n=1}^{N} F(ω)^*F(ω)}
$$
$$
H1(ω)= \frac{\sum_{n=1}^{N} (F(ω)^*V(ω)+F(ω)^*N(ω))}{\sum_{n=1}^{N} F(ω)^*F(ω)}
$$

ノイズ成分Nは無相関なノイズであることを前提としているため、平均化処理をすると0に近づきます。
イメージをつかみたい方は下記記事を参照ください。

MATLABで学ぶ信号処理 平均化によるノイズの低減(時刻歴編)
「平均化するとノイズの影響を低減できる」って良く聞きますよね? ただ、実測データを扱わない人(実験をしない人)はあまり理解できませんよね。 なぜかというと、そもそも扱うデータに誤差が混入していないからです。 本記事では、平均化によるノイズの低減(誤差の低減)について説明したいと思います。 今回は”平均化”についての1回目の記事ということで、はじめて勉強する方にイメージをつかんでもらうことを目的としています。次回以降でより専門的な内容について紹介していきます。

従って、Nを増加させて平均化すると、無相関ノイズNは0に収束するので下式になります。

$$
H1(ω)= \frac{\sum_{n=1}^{N} F(ω)^*V(ω)}{\sum_{n=1}^{N} F(ω)^*F(ω)}
$$

上記により、ノイズが低減できます

スポンサーリンク

MATLABでの検証

解析結果

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

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

図からわかるように、平均化処理をしていないFRFと比較して、理想的なFRFH1推定結果の誤差が小さくなっています。

ちなみに、応答Xの最大振幅の5%でノイズ成分を与えているので、もともと応答Xが小さかった周波数(例えば10Hz以上)では、相対的にノイズ成分が大きくなるので、平均化してもノイズが低減しにくいことがわかります。

コード

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

[h1] = H1(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(h1),'r--','linewidth',3)
hold off
grid on
legend('理想的なFRF','平均化処理をしないFRF','H1推定')
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

 

 

 

FRF MATLAB 信号処理
スポンサーリンク
MATLABパイセンをフォローする
MATLABパイセンが教える振動・騒音・音響・機械工学

コメント