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.
So far, we have introduced H1 and H2 estimation. This article describes H3 and Hv estimation.

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

MATLABパイセンが教える振動・騒音・音響・機械工学 | Mechanical Engineering University

H3 and Hv estimation of FRF

We have explained that it is desirable to use H2 estimation when the input contains a noise component, and
H1 estimation should be used when the response contains a noise component.

H3 estimation or Hv estimation is the estimation method to use when the input & response contain noise.

Incidentally, I found the following explanation in “Introduction to Modal Analysis” written by Dr. Nagamatsu, which I have been reading with great interest.

The frequency response function when errors are mixed in both input and output is H3, which is defined approximately as follows, and Hv, which is strictly defined.

So H3 is approximate and Hv is exact.
Incidentally, “Introduction to Modal Analysis” does not present any theoretical formula for Hv estimation.

Theoretical expression

H3 estimation

The formula for H2 estimation in the textbook itself is very simple.

$$
H3(ω)= \frac{H1(ω)+H2(ω)}{2}
$$

Hv estimation

The formula for Hv estimation is below.
$$
Hv(ω)= \sqrt{H1(ω)} \sqrt{H2(ω)}
$$

Validation in MATLAB

Results

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 blue line represents the FRF without averaging and with noise components.
(The H3 estimation result (with averaging) is represented by the green dotted line, the Hv estimation result (with averaging) is represented by the red dotted line.
(The Hv estimation result (with averaging) is shown by the red dotted line.

What do you think?
Compared to the ideal FRF, both H3 and Hv estimation errors are large above 10 Hz, where the response is small.
As Dr. Nagamatsu said, Hv estimation seems to be more accurate than H3 estimation.

 

MATLAB code

Executable file
cleaclear 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,:)+( rand(1,length(freq)) - 0.5 )*2*0.10; % Fを基準として±10%のノイズを加える
end

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

 

 

コメント