Half-power method


I am re-reviewing the bible-like book “Introduction to Modal Analysis,” which every Japanese vibration engineer has read at least once.

I would like to explain what I have learned in the re-review, showing the MATLAB code. This time, we will review from section 6.2.1, “Half-power method,” p. 344.

The following link is a summary to a description of vibration theory in English, if you are interested.

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

6.2.1 Method using the magnitude of the frequency response function (Half-power method)

The compliance (frequency response function of displacement/force) of a multi-degree-of-freedom system is Eq. (3.116). near the resonant frequency of the rth order, this rth order vibration mode is the dominant vibration component. Therefore, this resonance peak (around the peak resonance frequency) can be regarded as an approximate one-degree-of-freedom system. For example, in the case of compliance as shown in the figure below, each resonance peak can be regarded as a one-degree-of-freedom system.

G(ω)= x_j / F_i e^{jωt} = (\sum_{r=1}^{N} \frac{ φ_{ri} φ_{rj}  }{ -m_r ω^2 + j c_r ω + k_r }) (r=1~N)  \tag{3.116}

Figure 1: Compliance when resonance peaks are far apart

Note that when the two resonance frequencies are close and the resonance peaks are close to each other, as shown in the figure below, the system cannot be regarded as a one-degree-of-freedom system because of the influence of each other’s vibration components.

Figure 2 Compliance when resonant peaks are close together

The fact that equation (3.116) can be regarded as one degree of freedom means that ∑ and the subscript r can be omitted, so it can be expressed in equation (6.1).

G(ω) = \frac{ 1/K  }{ 1 – β^2 + 2jζβ }  , β=\frac{ ω }{Ω} \tag{6.1}

where the magnitude of compliance in equation (6.1) is equation (6.2).

|G(ω)| = \frac{ 1/K  }{ \sqrt{ (1 – β^2)^2 + (2ζβ)^2 } }   \tag{6.2}




“Mode identification” and “Curve-fitting” are applied to measurement data (experimental data). Therefore, mode identification for frequency response functions (compliance, etc.) obtained from an n-degree-of-freedom spring-mass model, as in this article, is not the original usage. (In the case of the spring mass model, it is nonsense to obtain mode characteristics such as resonance frequencies and mode vectors from the mass and stiffness matrices, so it is not necessary to go to the trouble of mode identification of frequency response functions.)
However, when understanding the principle of “mode identification,” it is easier to verify the results by using a spring mass model, so the spring mass model will be used in this explanation.

From here, the half-power method is explained. The half-power method is well known as a method for estimating the modal damping ratio from the frequency response function, but once the modal damping ratio is known, all of the modal characteristics can be obtained.
Figure 3 shows an enlarged view of the first-order resonance peak in Figure 1.

Figure 3 Enlarged view of the first-order resonance peak (m1 self-compliance)

Read the frequencies f1 and f2 at the two points that are 1/√2 of the maximum value|G|max in equation (6.2) (the two points that are -3 dB from the maximum value) and find the mode attenuation ratio ζ from the resonance frequency f0.

ζ= \frac{ f2-f1  }{ 2f_0 }   \tag{6.7.1}

Note that the resonance frequency f0 can be obtained without using the resonance frequency f0 using equation (6.7.2).

ζ= \frac{ f2-f1  }{ f2+f1 }   \tag{6.7.1}

The maximum value of compliance can also be expressed by Equation (6.5).

|G(ω)|_{max} ≒ \frac{ 1  }{ 2Kζ }   \tag{6.5}

Since the amplitude of the mode vector is 1 based on the case where this compliance G(ω) is self-compliant, K in Eq. (6.5) can be treated as equivalent to the modal stiffness k. Self-compliance is, for example, the transfer function obtained from the response of m1 when m1 is vibrated.

Once the modal stiffness k is obtained, the modal mass can be easily obtained from the following equation.

Ω^2 m = k

So far, we know the resonant frequency (the frequency of the maximum peak value, so we know it immediately), the modal damping ratio ζ, the modal stiffness k, and the modal mass m. All that remains is to find the mode vector, and then we will have obtained all the modal characteristics.

In the case of self-compliance, the component of the mode vector is now 1.
For example, in the case of Figure 3, the component of the mode vector corresponding to m1 is 1 at the resonance frequency of 9.432 Hz.
Then how do m2, m3, and …. are obtained?

It is obtained by measuring the mutual compliance Gis (the excitation point is i and the response point is s). φi=1, so the component φs of the mode vector becomes equation (6.11).

K_{is}= \frac{ k  }{ φ_s }   \tag{6.11}

The relationship to the maximum value of mutual compliance is equation (6.12).

|G_{is}(ω)|_{max} ≒ \frac{ 1  }{ 2K_{is}ζ }   \tag{6.12}

Therefore, the component φs of the mode vector can be obtained using equation (6.13).
φ_s ≒ 2kζ|G_{is}(ω)|_{max}   \tag{6.13}


Now let us solve Figure 3 as an example.

Answer(theoretical value)

The example problem in Figure 3 is a 4-DOF spring-mass model with the following modal characteristics
First-order resonance frequency: 9.43205759496419 Hz
Modal damping ratio: 0.1%
Modal mass: 1
Modal stiffness: 3512.14651375420
Mode vector:

Note that the mode vector, mode mass and mode stiffness values will change from the “answer” set here because the half-width method normalizes the excitation point of the mode vector as 1.

So, let’s vibrate m1 and make sure that the component of the mode vector of m1 is 1.
Modal damping ratio: 0.1%
Modal mass: 8.38384691873869
Modal stiffness: 29445.2987274969
Mode vector:


Value obtained by the half-power method

Resonance frequency f0: 9.4320Hz
Frequency at which the maximum value becomes 1/√2 at f1: 9.4226Hz
Frequency at which the maximum value becomes 1/√2 at f2: 9.4414Hz

From equation (6.7.1), the modal damping ratio ζ is about 0.1%.

zeta_hanchi=(f2 - f1)/(f0*2)
zeta_hanchi =


From equation (6.5), we obtain the modal stiffness k and the modal mass m.

Modal stiffness: 29545.6919999358
Modal mass: 8.41253424855546

From equation (6.13), let us find the mode vector.

phys =


We can see that we have correctly identified the mode characteristics.


Cautions for using the half-width method with experimental data

In this article, the frequency response function is obtained in MATLAB, so the frequency resolution can be fine. However, please note that it may not be possible to achieve fine resolution when using experimental data. If the frequency resolution is not fine enough, the frequency at half value is not known and the half-width method cannot be applied.

If you want to obtain a finer frequency resolution with experimental data, measure longer time history data and increase the number of FFT points (number of lines). If impulse excitation was performed and time history data was measured for only about 1 s, the frequency resolution can be increased by adding pseudo-silence (data with zero amplitude).

MATLAB code(Example)

 Executable code

Try other cases by changing the frequency resolution of the calculated frequency (freq) and the modal damping ratio (modal_dampimg).

clear all;clc;close all

freq=9.43-0.1:0.0001:9.43+0.1; %計算周波数の定義
w=2*pi*freq; %角振動数

[M]=eval_Mmatrix(m_vec); %質量行列
[K]=eval_Kmatrix(k_vec); %剛性行列

[V,D]=eig(K,M); %固有値D、固有ベクトルV
for ii1=1:length(m_vec)
    V(:,ii1)=V(:,ii1)/V(1,ii1); % 加振点m1のモードベクトルの成分が1になるように正規化
wn=sqrt(diag(D)); %固有角振動数
fn=wn/(2*pi); %固有振動数

mr=diag(V.'*M*V); %モード質量行列
kr=diag(V.'*K*V); %モード剛性行列
modal_dampimg=0.1*0.01; %モード減衰比(とりあえず0.1%と設定)
cr=c_cr*modal_dampimg; %モード減衰係数

F=zeros(length(m_vec),1); F(1)=1; %外力ベクトル(全周波数帯を1で加振する)
Xj=zeros(length(m_vec),length(freq)); %変位ベクトル(計算前にベクトルを定義してメモリを確保)

% % % % モード法
ii=1; %加振点 (加振点が複数ある場合はfor ii=1:kength(F)のループを追加すればよい )
for ii1=1:length(m_vec) %応答点
    for ii2=1:length(wn)
        Xj(ii1,:)=Xj(ii1,:)+V(ii1,ii2)*V(ii,ii2)./( -mr(ii2).*w.^2 + 1i*cr(ii2)*w +kr(ii2) )*F(ii);
%       横ベクトル =横ベクトル  +  スカラー  ×スカラー   ./( スカラー×横ベクトル + スカラー×横ベクトル +スカラー)×スカラー

hold on
xlabel('周波数 Hz')
ylabel('変位 [m]')

hold on
% legend('モード法')
xlabel('周波数 Hz')
ylabel('変位 [m]')
xlim([freq(1) freq(end)])
semilogy([8.5 10],[1 1]*max(abs(Xj(1,:)))/sqrt(2),'--k')

% % % % 半値幅法
zeta_hanchi=(f2 - f1)/(f0*2);

% Gis=max(abs(Xj(:,:)),[],2);
% phys=2*k_hanchi*zeta_hanchi*Gis
[Gis,id]=max( abs( imag(Xj(:,:)) ) ,[],2 );
Gis=Gis .* imag(-Xj(:,id(1)))./abs(Xj(:,id(1))) ;

function ffile

function [M]=eval_Mmatrix(m_vec)


function [K]=eval_Kmatrix(k_vec)

for ii1=1:length(k_vec)
    if ii1==1
        K(ii1-1:ii1,ii1-1:ii1)=K(ii1-1:ii1,ii1-1:ii1)+[1 -1;-1 1]*k_vec(ii1);