# H2 estimation of FRF

## Introduction

There are several types of transfer functions (Frequency Response Function (FRF)): H1 estimation, H2 estimation, H3 estimation, and Hv estimation. This article describes H2 estimation.

If you want to see all my lessons in English, please see below.

TopPage 振動騒音 研究所
MATLABとNVH(Noise, Vibration and Harshness)をこよなく愛する人生の先輩、それが私（MATLABパイセン）です。NVHを極めるにあたり、周辺の知識も網羅的に勉強することになり、その知識を共有すべく、本HPを運営しています。日本の製造業を応援すべく、機械エンジニアの「車輪の再開発」の防止と業務効率化の手助けをモットーに活動。専門である振動騒音工学や音響工学について、プログラムを示しながら解説。大学の授業よりもわかりやすく説明することを目指す。質問頂ければ、回答したいと思います。

## H2 estimation of FRF

As mentioned in the previous article, H1 estimation is the most common method of estimating the transfer function.

Frankly, it is not likely to be used at present: ……

First, let us consider an ideal FRF. The ideal FRF G(ω) is as follows. (入力＝input、応答＝response)

However, as a practical matter, such an ideal measurement is not possible. In the case of H1 estimation, we considered the case where noise enters the response as shown in the figure below. (入力＝input、応答＝response、ノイズ＝noise)

In the case of H2 estimation, consider the case of noise at the input as shown in the figure below. H2 estimation is a method to estimate FRF = (response V)/(input F) in a situation where the input signal contains noise components, as shown in the figure above.

At this point, you may be asking yourself the following question.

Background
If the input is given by generating vibration in a flow such as (signal generator) → (amplifier) → (shaker) as shown in the figure below, the input F can be measured from anywhere. Specifically, there are three possible ways to do this
(1) Branch the wire connecting the signal generator and amplifier, and measure the input F.
(2) Branch the wire connecting the amplifier to the shaker and measure the input F.
(3) Place an acceleration sensor on the shaker, or install a load cell in the shaker, or measure the input F by using a terminal for reference signal on the shaker itself, etc.

In this case, with ①, the signal is small and may be affected by noise due to the fact that it is branched for measurement. So, it is “not without”.

In the past, there were no power-supplied sensors such as ICPs, etc., so if input F was used as a sensor signal before passing it through an amplifier, etc., would you have used H2 estimation? (I’m guessing quite a bit here.)

## Theory of H2 Estimation

The formula for H2 estimation in the textbook itself is very simple. The theoretical formula is based on the following.

$$H2(ω)= \frac{W_{xx}(ω)}{W_{xf}(ω)}$$

where Wxf is the Cross Spectral Density Function of X and F and Wxx is the Power Spectral Density Function or Auto Spectral Density Function) of X.

However, the formula for Hz estimation in the textbook is a bit wrong. It is actually the following, where N is the number of times of averaging.

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

Or (what I mean below is that it can be obtained from the average of the cross spectral density function and the average of the power spectral density function)

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

Now we will explain why the above equation can reduce noise in the response.
Note that * refers to conjugation. (Note that * does not mean x.)

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

Since the noise component N is assumed to be uncorrelated noise, it approaches 0 after the averaging process. Therefore, when N is increased and averaged, the uncorrelated noise N converges to 0, resulting in the equation below.

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

Noise can be reduced by the above.

## Validation in MATLAB

#### Result

The results of the analysis are shown first. (100 times average) The ideal FRF with no noise component is shown by the black line.
The FRF with noise components and without averaging is represented by the blue line, and
（The H2 estimation result (with averaging) is shown by the red dotted line.

As can be seen from the figure, the error between the ideal FRF and the H2 estimation result is smaller than that of the FRF without averaging.

Note that the noise component is given by ±85% of F as a reference.

As we will explain in the next section, the coherence function γ^2 can be easily obtained once the H1 and H2 estimation can be programmed.

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

The coherence function in this study is shown in the figure below. #### MATLAB code

###### Executable file
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*0.85; % Fを基準として±85%のノイズを加える
end

[h1] = H1(ref,res,freq);
[h2] = H2(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),'g--','linewidth',3)
loglog(freq,abs(h2),'r--','linewidth',3)
hold off
grid on
legend('理想的なFRF','平均化処理をしないFRF','H1推定','H2推定')
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('FRF m/N');

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

###### function file
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 MATLABパイセンが教える振動・騒音・音響・機械工学