# H1, H2, H3 and Hv estimation of FRF

## Introduction

There are several types of transfer functions (Frequency Response Function (FRF)): H1 estimation, H2 estimation, H3 estimation, and Hv estimation. The four methods mentioned above have been introduced in the past. While summarizing these articles, I would like to talk about coding the program.

H1 estimation of FRF
There are several types of transfer functions (Frequency Response Function (FRF)): H1 estimation, H2 estimation, H3 estimation, and Hv estimation. This article describes H1 estimation. MATLAB code is also included.
H2 estimation of FRF
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. It also explains a little about the coherence function (Coherebce Function). MATLAB code is also included.
H3 and Hv estimation of FRF

#### Concept of MATLAB code

This is the flow chart of my MATLAB function of FRF analysis.

## Theoretical expressions for H1, H2, H3, and Hv

I will not explain it here, but only introduce the formula.

$$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(ω)}$$

## Comparative verification in MATLAB

#### Case (1) where noise is included in the input (Execution file (1))

Probably the H2 estimation would be the best result – right? Now let’s see what happens.
As a practical matter, it is hard to imagine that there will be noise in the input, but even if there is, it will be less than 1%.
In the case of a 1% error (noise), the result is shown in the figure below.

Figure 1: Error of 1%.

There is almost no difference. Since there is no use of consideration if there is no difference, I will give an error of 300% for now.

Figure 2: Error of 300%.

The order of the smallest difference from the ideal FRF is as follows
H2 estimation > H3 estimation > Hv estimation > H1 estimation

#### Case (2) where noise is included in the response (Execution file (2))

It is very common to get noise in the response.
In this case, for now, we gave an error of 3% based on the maximum amplitude X(ω). The result is shown below.

Figure 3: Error of 3%.

The order of the smallest difference from the ideal FRF is as follows
H1 estimation > Hv estimation > H3 estimation > H2 estimation

#### (3) Cases where noise is included in the input and response (executable file (3))

Although not realistic, the following figure shows the result when 300% noise is introduced into the input and 3% noise into the response.

Figure 4: 300% error in input, 3% error in response

The order of the smallest difference from the ideal FRF is as follows
Hv estimate > H3 estimate > H2 estimate > H1 estimate
(We placed importance on the value at the resonant frequency and the absence of peaks at other frequencies)

Considering the normal experimental environment, I imagine that even if noise is mixed in, it is about 3% for the input and 5% for the response, and compare the results under these conditions.

Figure 5: 3% error in input, 5% error in response

The order of the smallest difference from the ideal FRF is as follows
H1 estimate > Hv estimate > H3 estimate > H2 estimate
(We placed importance on the value at the resonant frequency and the absence of peaks at other frequencies)

In most FFT analyzers that can measure FRF, the default setting for FRF estimation method is “H1”. As shown in Fig. 3 and Fig. 5, “H1” estimation is considered to have the highest estimation accuracy in a normal measurement environment.

Under conditions other than those shown in Figures 3 and 5 (noisy environment), it is necessary to select an estimation method for each case.

#### Executable file (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));
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)])



#### executable file (2)

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)])


#### Executable file (3)

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


MATLABパイセンをフォローする
MATLABパイセンが教える振動・騒音・音響・機械工学