## Abstract

Output-only modal analysis (OMA) is an indispensable alternative to experimental modal analysis for engineering structures while in operation. Conventional OMA often fails to identify the underlying modal structure with insufficient modal energy contribution. Such low modal participation is expected when the sampled response is subjected to sensor nonlinearity or when specific modes are not directly excited. A novel subband decomposition (SBD) method that resolves modal parameters even with biased modal energy distribution is proposed. It isolates the system response within a narrow frequency subband through a finite impulse response analysis filter bank. Whenever the filter subband captures a resonance, the filtered system response is close-to-singular and contains mainly the resonant mode contribution. A modal cluster metric is defined to identify the resonant normal modes automatically. The modal parameters are also identified and extracted within the subband possessing the locally maximal clustering measure. The proposed method assumes no a priori knowledge of the structure under operation other than the system should have any repeated natural frequencies. Therefore, the SBD algorithm is entirely data-driven and requires minimal user intervention. To illustrate the concept and the accuracy of the proposed SBD, numerical experiments of a linear cantilevered beam with various stationary and non-stationary loading are conducted and compared to other OMA methods. Furthermore, physical experiments on an aluminum cantilever beam examine the method’s applicability in field modal testing. Compared to traditional OMA methods, the numerical and physical experiments show orders of magnitude improvement in modal identification error using the proposed SBD.

## 1 Introduction

Modal analysis is a well-adopted engineering practice for structural design and the dynamic characterization of engineering structures. Knowing the modal parameters during the design prevents catastrophic failures due to resonant deformations of structures, reduces excessive wear of machinery components, and mitigates harm to whoever operates the machines. Moreover, the obtained modal model helps explore structural responses to various loading scenarios and guides the structural-model-order reduction.

Modal testing resolves modal parameters (i.e., natural frequencies, modal damping ratios, and mode shapes) using experimentally obtained system responses. Experimental modal analysis (EMA) is a well-established engineering practice to identify modal parameters using both the system’s input and output measurements. Typically, the artificial excitation is an impulse force via a modally tuned hammer or a known excitation force provided by shakers. Frequency response functions are estimated based on the obtained modal testing data sets, and frequency domain modal models can be fit to the estimated functions to identify the poles and the modal matrix of the structure. The most frequently used modal parameter identification methods include the peak picking method [1,2], the least square complex exponential method [3], and the least square orthogonal rational function (LSRF) method [4].

In general, EMA is almost only applicable during the design stage of the structure; otherwise, it will impose unnecessary downtime when the system is in operation. To tackle modal testing during operation, output-only modal analysis or operational modal analysis (OMA) was established by using only the response data. OMA has received much attention for its use in on-site field measurements and real-time applications such as structural health monitoring and vibration control [5,6]. OMA algorithms can be grouped into the time- and frequency-domain approaches. For the latter, the basic frequency domain or peak picking technique is used to estimate the natural modes and frequencies from the power spectral density matrix of the response measurements [7]. For closely-spaced modes, this method deteriorates, and the results become biased. The frequency domain decomposition (FDD) method compensates for the deficiencies in the basic frequency domain method without losing its simplicity [8]. In addition, the extension of this method improves the modal damping and frequency estimations and automates the peak picking procedure [9]. However, when the resonant peaks in the power spectral density function are contaminated by non-modal interference, the modal parameters cannot be resolved. The Ibrahim time domain and eigen system realization algorithm (ERA) [10] developed in the 1970s and 1980s are widely used time-domain algorithms. However, these methods assume that the system is excited by either white Gaussian noise or an impulse (including a free-response case). The general realization algorithm extends the ERA method to allow for modal parameter estimation for arbitrarily excited systems. However, it uses the input signal to construct a weighted Hankel matrix, rendering it an EMA method [11].

It is also well recognized that the OMA-resolved modal parameters are primarily affected by the mass distribution of the structure (i.e., we experience mass scaling of the mode shapes) and the extraneously low modal participation usually caused by excitation and observer biases. Such biases affect the use of the OMA methods that directly exploit the covariance/correlation matrices of the response data. A well-known example is the proper orthogonal decomposition (POD) which uses the singular value decomposition (SVD) of the response matrix [12] as a solution. Although widely used for modal identification, it is only suited for systems with uniform mass distribution. Otherwise, the corresponding non-uniform mass matrix must be known a priori to accurately estimate the linear normal modes (LNM) from the raw time series, which is impractical for complicated mechanical systems. This problem is resolved by the smooth orthogonal decomposition (SOD) using the generalized singular value decomposition (GSVD) to find a set of smooth orthogonal temporal coordinates of the sampled response data [13]. The SOD could accurately estimate an undamped system’s LNMs and natural frequencies without knowing the mass matrix. However, when the system under investigation is in operation, not all the LNMs are necessarily excited. The operation-biased system response impairs the modal identification quality when the decomposition is solely based on the observed energy (e.g., the covariance matrix). Therefore, energy-based OMA methods like the POD and SOD tend to fail when the observation is biased by the sensors themselves or by extraneous non-system events. This can be caused by biased sensing, localized electrical noise, and unexpected intermittent events in real-world scenarios. As a result, essential modal information cannot be directly extracted from the system response. Hence, there is a need for a new method to accommodate the deficiencies experienced in POD and SOD while preserving the simplicity of the algorithms.

Some early work tried to resolve the biased modal energy problems by isolating the modal identification process within a small frequency bandwidth. One method designs bandpass filters to isolate individual modes allowing the successful extraction of LNMs in structures with non-uniform mass distribution using POD [14]. The time domain decomposition uses similar ideas in which digital bandpass filters were designed to isolate a system’s responses, where each behaves like a single-degree-of-freedom system [15]. However, both proposed methods require user input to design the filters for each mode in the response data and require users’ judgment on where the modal frequency is located. Recently, filter banks have been applied in structural health monitoring and fault diagnosis for mechanical systems [16]. However, to the author’s best knowledge, no study has explored using filter bank structures to assist modal parameter identification.

This paper aims to isolate and identify modal structures located at or near resonance frequencies where the dynamics of the system are known to be singular (i.e., the system response is composed of or dominated by only one mode). In this new method, the total bandwidth of the response is partitioned into predetermined subbands via a uniform analysis filter bank. Afterward, the SOD is applied to each subband to resolve the underlying low-rank modal structure. Unlike the SVD applied in time domain decomposition, the basis found from SOD is generally not orthogonal. For bandlimited signals at a resonant frequency, the response should be dominated by only one mode, and one should expect the matrix of modes obtained by SOD to be close to rank deficient. By exploiting this property of SOD, resonance peaks and modal structures can be identified using a classification criterion to find the close-to-singular modal structures within the subbands. To justify the advantages of the proposed method, numerical and physical experiments will be conducted to emulate the excitation bias resulting from the limited bandwidth of the environmental loading (e.g., non-white Gaussian loading) or the operational loading correlated to the system response (e.g., harmonic excitation) and the observer bias due to sensor nonlinearity for specific frequency range(s).

The remainder of the paper is arranged as follows: Sec. 2 provides a brief overview of the preliminaries, particularly the POD and SOD algorithms, and directs the readers to the relevant literature for further details; Sec. 3 formulates the proposed methodology; Sec. 4 presents numerical and experimental examples are explored to compare the newly proposed methodology and some popular OMA methods; afterward, a thorough discussion and important remarks are followed by the cConclusions.

## 2 Popular Output-Only Modal Analysis Methods and Their Limitations

For a given vibratory mechanical system, its dynamics—when subjected to a small perturbation—can be assumed to be a linear combination (i.e., superposition) of a set of time-invariant spatial functions and their corresponding time coordinates. Therefore, one looks for such coherent spatiotemporal structures through the separation of variables. In practice, this separation of variables is conducted by decomposing a matrix that contains a sample of the underlying continuous field of the structural response in time. As a result, we consider a sampled scalar field in the form of a data matrix *X* ∈ ℝ^{m×p} for *m* time snapshots of position, velocity, or acceleration measurements at *p* distinct spatial points. Since no input information is directly given, for almost all output-only modal analysis methods, spatiotemporal dynamics separation is realized by investigating the correlation between the acquired system response and the transformed version of the correlation matrix (e.g., power spectral density estimates in FDD). The POD, the SOD, and the ERA are considered as examples to illustrate the limitations of some standard output-only modal analysis methods and show the necessity and advantage of the proposed modal analysis framework. These methods are widely recognized in OMA, are easy to implement, and do not require special processing of the data (e.g., transformation into frequency domain).

### 2.1 Proper Orthogonal Decomposition.

*X*, i.e.,

*U*are the proper orthogonal coordinates,

*Σ*= diag(

*σ*

_{i}) is a diagonal matrix where the square of the- diagonal elements, $\lambda i=\sigma i2$, are the proper orthogonal values, and the columns of orthonormal

*V*are the proper orthogonal modes (POMs). It is important to note that the

*λ*’s are ordered based on modal energy content, i.e.,

*λ*

_{1}≥

*λ*

_{2}≥, …, ≥

*λ*

_{p}. In the context of modal analysis, the most energetic mode in a given experiment or simulation will have the largest

*λ*. Therefore, the energy (or the variance) of data dictates the resolved spatial functions (modes) and the temporal functions (temporal coordinates); as a result, the modal participation affects the ability and accuracy of the resolved proper orthogonal modes from a given data set [12].

### 2.2 Smooth Orthogonal Decomposition.

*Φ*, where a projection of the original data matrix onto this basis arranges the time coordinates according to their maximal smoothness [13] while keeping their variances the same. The degree of smoothness of each time coordinate is described by the variance of the corresponding time derivative approximated by some finite difference operator, $D$. The SOD problem is solved by performing a generalized singular value decomposition (GSVD) on the matrix pair

*X*and $X\u02d9$, where $X\u02d9=DX$

*U*∈ ℝ

^{m×p}and

*V*∈ ℝ

^{m×p}are matrices with orthonormal columns,

*C*= diag(

*c*

_{i}) ∈ ℝ

^{p×p}and

*S*= diag(

*s*

_{i}) ∈ ℝ

^{p×p}are diagonal matrices, and the common spatial basis of the decomposition,

*Φ*∈ ℝ

^{p×p}are the

*smooth modes*(SMs). The SMs are invertible when

*X*is full-rank [19,20] (but not necessarily orthogonal), and their corresponding temporal coordinates in

*UC*are called the smooth orthogonal coordinates (SOCs). The SOCs are conventionally sorted according to their “smoothness,” characterized by $\mu i=ci2/si2$ smooth values. SOD resolved modes are both energetic and smooth. Therefore, the low-participation modes are more likely to be resolved as long as they are smoother than the more energetic modes or noise.

### 2.3 Eigensystem Realization Algorithm.

*A*and an initial state

*u*

_{0},

*O*is the observation matrix and $Du0=exp(A\Delta t)u0$ is the solution to the underlying linear time-invariant first-order system $u\u02d9=Au$ with initial condition

*u*

_{0}. The algorithm of ERA, in general, considers an ensemble of system responses

*Y*(

*n*) = [

*y*

_{1}(

*n*),

*y*

_{2}(

*n*), …] evaluated at the sample time

*n*and their corresponding initial conditions

*U*

_{0}= [

*u*

_{01},

*u*

_{02}, …]; where the following relationship holds,

*Y*, two Hankel matrices that are one sample time away from each other are created and decomposed as follows [10,22],

*P*= [

*O*,

*OD*,

*OD*

^{2}, …]

^{T}is the observability matrix and

*Q*= [

*U*

_{0},

*DU*

_{0},

*D*

^{2}

*U*

_{0}, …] is the controllability matrix. The unbiased estimates of

*P*and

*Q*are usually obtained through the SVD of $H(0)=UH0\Sigma H0VH0T$ and by letting $P^=UH0r\Sigma H0r1/2$ and $Q^=\Sigma H0r1/2VH0rT$, respectively [24,25]. Then, the estimated system dynamics can be obtained as

^{†}being the Moore-Penrose pseudo-inverse, and

*r*is the order of realizations of the system matrix (where

*r/*2 is the number of modes considered in the above full-rank decomposition). Then, the eigen decomposition to the estimated

*r*-realization system dynamic matrix can be accomplished via [24,25]

Although widely adopted in the contemporary OMA, the aforementioned methods may fail to extract the modal parameters when there is low modal participation, observer-biased modal energy distribution, or energy leakage from non-impact excitation. For example, if one of the modal structures is not properly excited during the recording and its energy contribution to the overall system response is far less than other modes, the POD may completely fail to extract this mode. Even though smooth orthogonal decomposition was proposed to resolve non-orthogonal modes by rendering the resolved temporal basis functions (temporal coordinates) to be maximally smooth, it may still suffer from the above-mentioned issues. Numerical examples will be illustrated in regard to these problems in Sec. 4.

## 3 The Proposed Subband Decomposition Algorithm

The idea behind the subband decomposition (SBD) algorithm is to isolate the modal structure within a certain frequency range that is close to a resonant frequency. When at the damped natural frequency, the system’s frequency response function is expected to be close to singular (and it is singular when there is no damping) and the system response is dominated by the mode associated with that resonant frequency [24,26]. Past research has used filter theory to isolate these modal structures but relies on the user to manually locate the resonances in the frequency domain beforehand [14,15]. This section proposes an algorithm to partition the vibrational bandwidth into many subbands and identify the subbands containing singular modal structures based on a new and robust criterion.

*R*and the poles

_{l}*p*=

_{l}*jω*; and (·)

_{dl}^{*}denotes the complex conjugation. This assumption implies that whenever the system response is close to a resonant frequency [13], the resolved system dynamics is composed of close to only one modal structure and its corresponding damped natural frequency by assuming no closely-spaced modes are present.

If a set of collected spatiotemporal response data *X* is filtered to present only in a limited bandwidth, such filtered response *Y _{i}* in the frequency domain can be treated as a linear combination of a limited number of frequency components (ideally, only one frequency component), each of which shares modal contribution from all the modes present in this limited bandwidth. Whenever the subband approaches a resonance, the contribution from that resonance dominates the current filtered data set. By applying the SOD to the subband system response matrix

*Y*which is centered at a resonant frequency, the SOD will identify mode shapes that are similar to each other; in other words, the resolved SMs will span a closely-spaced modal subspace (if not lying on the same modal subspace when there is no noise and only one mode is active).

_{i}### 3.1 Finite Impulse Response Bandpass Filter Design.

Ideally, the proposed method decomposes the obtained system response data into a sum of bandpass-filtered response matrices without energy leakage. Therefore, a proper digital bandpass filter must approximate the ideal filter. The filter design has to isolate the frequency component of the response data so that when the center frequency of the bandpass filter approaches the resonance, close-to-singular system dynamics can be identified by the SOD. Although a direct application of bandpass filtering in the frequency domain (e.g., using the fast Fourier transform) may be used, in the proposed algorithm, a filter bank is designed before the subband decomposition without explicitly transforming the response data between time and frequency domain to reduce the computational overhead. The narrow-band decomposition of the trajectory matrix is conducted in the time domain via convolution by first designing a finite impulse response (FIR) analysis filter bank. An FIR filter bank is preferred due to its balance of stability, linearity, and ease of design.

*ω*

_{1}and

*ω*

_{2}are the lower and upper cutoff frequencies of the bandpass filter, respectively.

*h*[

_{D}*n*] is the discrete-time ideal impulse response of the filter. To obtain the actual filter coefficients, we must multiply

*h*[

_{D}*n*] by some window function,

*w*[

*n*].

*A*, is a function of the passband and stopband ripples,

*δ*and

_{p}*δ*, respectively. This stopband attenuation is found through

_{s}*A*= −20log

_{10}(

*δ*) where

*δ*= min(

*δ*

_{s},

*δ*

_{p}). The passband ripple,

*δ*, does not have any strict specifications for our application since we are only concerned with rejecting out-of-band modal contributions. Hence, the stopband ripple plays an important role in isolating vibration modes and must be selected to satisfy the stopband attenuation,

_{p}*A*. After which, the filter-order,

*N*(i.e., number of coefficients), can be found through

*N*≥ (

*A*− 7.95)/(14.36Δ

*f*) where Δ

*f*is the normalized transition width. The transition width is specified by the user and can be selected to be as “sharp” as needed. Note that a sharp transition width prevents modal contributions at the edge of the bandpass filter from leaking into the subband. The piecewise relationships that govern the Kaiser window are given by [27] and stated below

*I*

_{0}is the zero-order modified Bessel function of the first kind and

*N*is the number of filter coefficients needed to satisfy the design parameters. The ripple control parameter is then computed as,

*A*is the stopband attenuation. The FIR filter coefficients can then be obtained by multiplying the ideal impulse filter response and the window function,

The original trajectory matrix is in the form of a tall-and-skinny matrix with *m* ≫ *p*. Each column in *X* is convolved with the filter coefficients in Eq. (13) to filter the entire response,

*Y*∈ ℝ

^{(m+p−1)×p}is the output of the filtered data matrix

*X*and

*N*is the number of coefficients needed for the proposed filter to satisfy the design specifications.

### 3.2 Filter Bank Construction and Analysis.

We aim to partition the bandwidth of a system response signal into *M* subbands, where each subband contains a bandlimited spectrum governed by the lower and upper cutoff frequencies (i.e., *ω*_{1} and *ω*_{2}) of the *M* bandpass filters in the array, see Eq. (9). This step determines the quality of the resolved subband SMs, which affects the accuracy of determining the modal structures from the obtained data. To help better conceptualize the proposed subband decomposition algorithm, Fig. 1 visually summarizes the proposed decomposition method.

*M*impulse response functions separately. The window function and design criteria remain constant throughout the banks; however, we must slide the passband region over the usable bandwidth,

*i*corresponds to the

*i*th passband region for that bandpass filter and

*h*[

_{i}*n*] corresponds to the FIR filter coefficients of each bandpass filter in the filter bank. An overlap ratio can be specified between adjacent subbands if desired. Then the discrete-time convolution is performed between the data matrix,

*X*∈ ℝ

^{m×p}, and the FIR coefficients

*h*[

_{i}*n*] for each subband,

*Y*,

_{i}The common basis for each bandlimited response, *Φ _{i}*, contains the SMs, which are expected to be rank deficient if the subband is centered on a resonant frequency. Only the SM matrices for each subband are stored in the aforementioned procedure. Thus, the algorithm does not have demanding memory requirements. The authors considered a maximum 3 × 10

^{4}subbands without any memory storage issues. However, the maximum number of subbands allowed will depend on the user's computer system. The algorithm can be parallelized for accelerated computation for clusters with enough random access memory.

### 3.3 Subband Cluster Measure for Modal Identification.

Since the original trajectory matrix is partitioned into many subbands, it is important to have an automated selection process to identify the resonant subbands. When decomposing a subband that overlaps with a resonant frequency using SOD, the SMs are expected to exhibit similar displacement amplitudes. This obviously assumes minimal leakage of other modes into the subband. The SVD (or POD) does not share this property since the right singular vectors in *V* have to be orthogonal. Therefore, if one were to apply the POD to the same subband, regardless of whether the subband is centered at a resonant frequency or in between two resonances, the resolved mode shapes would always be orthogonal, and the proposed identification metric will fail. The clustering measure exploits this characteristic of SOD to identify the resonant subbands across the entire partition discussed in Sec. 3.2. The resolved SMs in a resonant subband are expected to be closely-spaced or clustered. Therefore, this clustering of resolved SMs is exploited to identify the underlying modal structure within each subband.

If the identified SMs in *Φ _{i}* are highly correlated, the underlying dynamics is close to singular (or rank-1). In contrast, there are no resonant modes whenever the resolved subband SMs are less clustered. Since SOD is less sensitive to the variance of the data, the subband SOD helps isolate the mode in the current subband by negating the leakage effects due to the highly energetic out-of-band modes.

To characterize the clustering behavior of the subband SMs, a weighted angular similarity between identified SMs in *Φ* ∈ ℝ^{p×p} (i.e., column-wise) is considered. Note that this cluster measure is to be computed for $\Phi i\u2200i=1,2,..,M$. First, column-wise cosine similarity is computed

*l*

^{2}norm, see Fig. 2(a). The absolute value is used to constrain the output between 0 and

*π*/2 degrees. It is desirable to normalize

*S*

_{jk}by

*π*/2

*Θ*

_{jk}∈ ℝ

^{p×p}is a normalized angular similarity matrix whose elements are bounded between zero and one. An element with a similarity of zero corresponds to a pair of orthogonal vectors being compared, while a value of one corresponds to a set of linearly dependent vectors. It is important to note that

*Θ*

_{jj}= 1 and

*Θ*

_{jk}=

*Θ*_{kj}, which is the result of finding the inner product of the same vector. When the SMs are a complete set of orthogonal vectors in the column space, the resultant

*Θ*will be the

_{jk}*p*×

*p*identity matrix,

*I*. Furthermore, when the SMs are all linearly dependent, then the resultant

_{p}*Θ*will be the

_{jk}*p*×

*p*matrix filled with ones,

*J*.

_{p}To differentiate the subbands encompassing a resonant frequency even further, an exponential function of the form $h(\Theta jk)=aexp\u2212b\Theta jk$ is used to weight the elements in *Θ _{jk}*. To use the same scale, the following condition is imposed:

*h*(0) = 1, which results in

*a*= 1. The parameter

*b*is found from the condition

*h*(1) =

*ɛ*for

*ɛ*≪ 1 and is left for the user to tune in to the post-processing stage if needed.

Transforming the similarity measure *S _{ij}* to the angle

*Θ*provides a linear metric that characterizes the distance between the two SMs within a subband. The imposition of the weight

_{ij}*h*(·) provides a family of similarity metrics that emphasizes the location of the modal structure by penalizing the similarity exponentially. Figure 2(b) shows similarity functions with and without the weighting as the angle between two arbitrary SMs approaches

*π*/2. For the studies in Sec. 4, we choose

*b*= 10 for the weighting function, and for the experimental study in Sec. 5, we choose

*b*= 5. This parameter

*b*is what further distinguishes subbands containing a resonance from non-resonance. It is important to note that the weighting will have little effect when all SMs are linearly dependent on a subband. This weighting only affects SM matrices that do not contain a close to singular structure, thus distinguishing the non-resonant subbands from the resonant subbands.

*Θ*are then summed together to reduce the classification criterion to a scalar. The maximum of this sum will depend on the number of measurement points used in the study, which is

_{jk}*p*

^{2}for

*p*distinct measurement locations. Thus, we propose to divide the sum by

*p*

^{2}, which normalizes the criteria such that the maximum is 1, i.e.,

*w*is bounded between

_{i}*p*is the column dimension of the SM matrix,

*Φ*. The criterion is applied to the SM matrix for every subband. The result is an array of elements, $wi\u2200i=1,2,\u2026,M$ that quantify the linear dependence (singularity) of the SMs in the effective bandwidth of the response. An ideal

*w*as a function of the frequency (or subband index) for a system with three distinct modes is shown in Fig. 2(c). A graphical illustration of obtaining the subband clustering measure is summarized in Fig. 2.

### 3.4 Linear Normal Mode and Damped Natural Frequency Identification.

*w*is directly proportional to the damped natural frequency of the resonance. Those with the largest magnitude correspond to a close-to-singular structure and are identified as resonance subbands. Then, the natural frequency (Hz) information is given as

_{i}*BW*is the bandwidth of the bandpass filters with units in Hz, and

*i*is the subband index when

*w*is a local maximum. In ideal cases, the subbands corresponding to resonance will have

*w*=

_{i}*p*

^{2}. However, this is seldom the case due to damping and noise in the system. Despite this, the information in the

*w*vector indicates where the resonances occur. The corresponding modes are contained in the close-to-singular SM matrix,

_{i}*Φ*. Since there are

*p*candidate vectors in

*Φ*, the SVD of

*Φ*=

*UΣV*

^{T}is performed, and the first column in the left singular vectors,

*U*(i.e., dominant mode shape), is selected as the LNM. Alternatively, we could also average all columns in close-to-singular

*Φ*to estimate the corresponding LNM. For free decay responses, the damping ratios can be obtained from the smooth orthogonal coordinates within each identified resonant subband using standard methods such as the logarithmic decrement. For the forced vibration method, the half-power method could be applied to the FRF of the smooth orthogonal coordinates.

## 4 Numerical Experiments

In this section, the authors explore the application of the proposed SBD algorithm for numerous situations in which the conventional OMA methodologies usually fail. We consider a uniform cantilever beam that is clamped at *x* = 0 and free at *x* = *L*. The parameters of the beam are kept identical to the ones used in Ref. [28], where the stiffness is *EI* = 1, the mass per unit length *m* = 1, and a unit length of *L* = 1. A total of ten Euler Bernoulli beam elements are used to construct the system’s mass and stiffness matrices, *M* and *K*, respectively. A Rayleigh damping model is used in numerical experiments. The damping matrix, *C*, can be approximated by *C* = *α*_{1}*M* + *α*_{2}*K*. The first and tenth modes are given the modal damping ratio of *ζ*_{1} = 5.0 × 10^{−3} and *ζ*_{10} = 1 × 10^{−3}, respectively. After which, the Rayleigh coefficients are computed to be *α*_{1} = 5.6 × 10^{−3} and *α*_{2} = 1.36 × 10^{−5}. The natural frequencies and damping ratios of the finite element cantilever beam model are given in Table 1. The prescribed modal damping is aimed to emulate an intrinsic low modal participation scenario. This is a consequence of the exceedingly large modal damping of the first mode, whose damping ratio is at least five times greater than any other mode. As a result, the first mode is expected to have very low energy associated with its modal motion.

Modal parameters for finite element beam | ||
---|---|---|

Mode # | f_{n} (Hz) | ζ |

1 | 0.560 | 0.005 |

2 | 3.507 | 0.000821 |

3 | 9.822 | 0.00058 |

4 | 19.261 | 0.000276 |

5 | 31.888 | 0.000304 |

6 | 47.773 | 0.000383 |

7 | 67.028 | 0.000497 |

8 | 89.799 | 0.000641 |

9 | 116.109 | 0.000813 |

10 | 144.340 | 0.001 |

Modal parameters for finite element beam | ||
---|---|---|

Mode # | f_{n} (Hz) | ζ |

1 | 0.560 | 0.005 |

2 | 3.507 | 0.000821 |

3 | 9.822 | 0.00058 |

4 | 19.261 | 0.000276 |

5 | 31.888 | 0.000304 |

6 | 47.773 | 0.000383 |

7 | 67.028 | 0.000497 |

8 | 89.799 | 0.000641 |

9 | 116.109 | 0.000813 |

10 | 144.340 | 0.001 |

*i*, and

*x*is the magnitude of the normalized FEM mode shape for the same spatial point; and

_{i}*q*= 10 stands for the total number of nodes considered in this study.

### 4.1 Cantilever Beam: Impulse Excitation.

A single impulse is applied to the first transverse displacement node of the beam. The system’s impulse response is recorded for 60 s with a sampling frequency of 1000 Hz. The SOD, POD, ERA, and the proposed SBD methodologies are applied to the response matrix to estimate the LNMs of the structure. For the case of SBD, the filter bank parameters in Table 2 are used to construct the coefficients for each filter in the filter bank. The system response is then convolved with the filter bank using Eq. (16). In this example, the first four modes are identified. The error between the underlying FEM modal structure and natural frequencies is compared to those identified by the algorithms. Table 3 summarizes the error between the estimates of the LNMs and damped natural frequencies with that of the FEM model. The proposed identification criterion and identified first LNM are shown in Fig. 3.

RMSE lightly damped cantilever beam | ||||||||
---|---|---|---|---|---|---|---|---|

Mode # | $E\varphi $ | $E\omega $ | ||||||

SOD | SBD | POD | ERA | SOD | SBD | POD | ERA | |

1 | 7.7 × 10^{–2} | 7.2 × 10^{–4} | 1.24 × 10^{–1} | 4.35 × 10^{–1} | N/A | 9.6 × 10^{–3} | N/A | N/A |

2 | 6.5 × 10^{–3} | 7.0 × 10^{–4} | 2.2 × 10^{–1} | 8.67 × 10^{–2} | 1.4 × 10^{–1} | 4.3 × 10^{–2} | N/A | 2.9 × 10^{–1} |

3 | 2.9 × 10^{–3} | 1.4 × 10^{–3} | 2.43 × 10^{–1} | 6.4 × 10^{–3} | 6.8 × 10^{–2} | 2.8 × 10^{–2} | N/A | 7.2 × 10^{–3} |

4 | 5.4 × 10^{–3} | 1.8 × 10^{–3} | 2.16 × 10^{–1} | 3.5 × 10^{–3} | 1.8 × 10^{–2} | 1.1 × 10^{–2} | N/A | 1.5 × 10^{–4} |

RMSE lightly damped cantilever beam | ||||||||
---|---|---|---|---|---|---|---|---|

Mode # | $E\varphi $ | $E\omega $ | ||||||

SOD | SBD | POD | ERA | SOD | SBD | POD | ERA | |

1 | 7.7 × 10^{–2} | 7.2 × 10^{–4} | 1.24 × 10^{–1} | 4.35 × 10^{–1} | N/A | 9.6 × 10^{–3} | N/A | N/A |

2 | 6.5 × 10^{–3} | 7.0 × 10^{–4} | 2.2 × 10^{–1} | 8.67 × 10^{–2} | 1.4 × 10^{–1} | 4.3 × 10^{–2} | N/A | 2.9 × 10^{–1} |

3 | 2.9 × 10^{–3} | 1.4 × 10^{–3} | 2.43 × 10^{–1} | 6.4 × 10^{–3} | 6.8 × 10^{–2} | 2.8 × 10^{–2} | N/A | 7.2 × 10^{–3} |

4 | 5.4 × 10^{–3} | 1.8 × 10^{–3} | 2.16 × 10^{–1} | 3.5 × 10^{–3} | 1.8 × 10^{–2} | 1.1 × 10^{–2} | N/A | 1.5 × 10^{–4} |

### 4.2 Cantilever Beam: Single Harmonic Excitation.

In certain situations, modal tests cannot be conducted using broadband excitation signals (e.g., a sinusoidal sweep or random signals), especially when the excitation system has limited bandwidth or is generated by an open-loop control (the excitation may not be what is prescribed from the input source). Moreover, in operating reciprocating machinery, harmonic forcing, and environment loading appear simultaneously [25]. In such cases, it is desirable if a narrow-band excitation or a single harmonic excitation can be used for modal identification purposes. Here, a single harmonic excitation is used to excite the base of the cantilever beam model. The excitation is of the form *f*_{e}(*t*) = *A*sin(2*π ft*), where different excitation frequencies, *f*, are chosen to understand the effect on the proposed SBD’s identification process (Sec. 3.3). For all cases, the excitation amplitude is set to *A* = 1. The system is simulated for 100 s using a sampling frequency of 1000 Hz. The whole time series for each forcing frequency is used in the analysis.

The natural frequencies for the modes corresponding to the ten transverse degrees-of-freedom are between 0 Hz and 150 Hz for the unit parameters. For modal tests in the field, the structure’s resonance frequencies may not be known ahead of time, and therefore the selection of an appropriate forcing frequency may not be possible. We consider four different forcing frequencies when *f* = [25, 75, 125, and 400] Hz. Figure 4 shows the proposed resonant subband identification criterion for each case of forcing frequency. These results demonstrate the forcing frequency has a significant effect on the performance of the criterion. In all cases, no LNMs above the forcing frequency cannot be identified. Moreover, for the 75 Hz and 125 Hz cases, the surrounding LNMs are affected by the forcing frequency. For the 400 Hz harmonic input, the criterion indicates that the first nine modes are clearly identified and are successfully extracted by the algorithm. The estimates for modal vectors and natural frequencies in this section are based on the 400 Hz excitation. Figure 5(a) shows the error between the estimated modal vectors for each mode. The combined errors for both the estimated modal vectors and natural frequencies for each methodology are summarized in Table 4.

Average of RMSE under single harmonic excitation | ||||||
---|---|---|---|---|---|---|

Mode # | $E\varphi $ | $E\omega $ | ||||

SOD | ERA | SBD | SOD | ERA | SBD | |

All | 8.43 × 10^{–2} | 6.26 × 10^{–2} | 4.3 × 10^{–3} | 5.45 × 10^{–1} | 1.43 × 10^{–1} | 2.65 × 10^{–2} |

Average of RMSE under single harmonic excitation | ||||||
---|---|---|---|---|---|---|

Mode # | $E\varphi $ | $E\omega $ | ||||

SOD | ERA | SBD | SOD | ERA | SBD | |

All | 8.43 × 10^{–2} | 6.26 × 10^{–2} | 4.3 × 10^{–3} | 5.45 × 10^{–1} | 1.43 × 10^{–1} | 2.65 × 10^{–2} |

### 4.3 Cantilever Beam: Non-Stationary Excitation.

*A*is the forcing amplitude,

*c*= (

*f*

_{1}−

*f*

_{0})/

*T*is the chirp rate defined by difference between the final frequency,

*f*

_{1}, and starting frequency,

*f*

_{0}, divided by the total sweep time. The phase of the sinusoid is

*ϕ*

_{0}. A simple case of

*f*

_{0}= 0 and

*ϕ*

_{0}= 0 is considered here. The simulation is sampled at 500 Hz and the input sweeps through the entire effective bandwidth (0 → 250 Hz) in

*T*= 100 s. The sweep rate for this signal is

*c*= 2.5. The beam is excited by a full sweep through the vibrational signal’s entire effective bandwidth (i.e., 250 Hz). Figure 5(b) shows the errors computed by Eq. (23) for each mode compared to the FEM LNMs.

## 5 Experimental Study

An experimental study was performed on a uniform cantilever beam to investigate the effectiveness of the SBD algorithm in an experimental setting. A uniform aluminum beam that is clamped at *x* = 0 and free at *x* = *L* (assuming cantilever boundary conditions) was considered. The beam dimensions are 610 × 25.4 × 3.175 mm. Six accelerometers (PCB U33B42) that are equally spaced across the beam span are used to record the acceleration response of the beam under impulses. The system is excited by an impulse train via a modally tuned hammer (PCB 086C01) across the mid-span of the beam. The roving hammer technique [29] introduces redundancy and more sample responses for averaging when performing modal parameter identification. This guarantees more accurate frequency response function estimates and leads to a more reliable experimental benchmark for comparison.

The results obtained by the LSRF method serve as the benchmark for comparison for the SBD and SOD algorithms, therefore, the identified modes and frequencies using this method are taken as the ground truth. The same procedure is applied to the experimental data from Sec. 4 using the filter bank parameters shown in Table 5. Since only six accelerometers are used during the experiments, a maximum of six modes can be resolved from the obtained data. The sampling frequency is set to 1000 Hz, and an anti-aliasing filter was set to 500 Hz to ensure no higher-order modes or high-frequency noise interferes with the modal parameter identification. The first two modes of the resolved modal structure and identification criterion are shown in Fig. 8.

Filter bank parameters for the experimental study | ||||
---|---|---|---|---|

Bandwidth (Hz) | Stopband attenuation (dB) | Stopband ripple (dB) | Transition width (Hz) | Overlap |

0.1 | −100 | 1 × 10^{–5} | 0.05 | 0 |

Filter bank parameters for the experimental study | ||||
---|---|---|---|---|

Bandwidth (Hz) | Stopband attenuation (dB) | Stopband ripple (dB) | Transition width (Hz) | Overlap |

0.1 | −100 | 1 × 10^{–5} | 0.05 | 0 |

## 6 Discussion

The numerical and physical experiments demonstrate the applicability of the SBD algorithm for several practical situations that arise in OMA. The objective of the above studies is to illustrate situations in which the considered conventional OMA algorithms discussed in Sec. 2 fall short of identifying all of the LNMs. The identified shortcomings of the considered OMA algorithms are due to: (1) the excitation-induced low modal participation or (2) the biased response variance introduced by sensor nonlinearity during operational modal testing. The following arguments relating to the considered experiments justify the benefits of the new proposition:

The first numerical example explored the performance of the SOD, POD, ERA, and proposed SBD algorithms by identifying the lower-order LNMs of a cantilever beam for a single impulse applied close to the base. Although the system is lightly damped, both the SOD and ERA still cannot resolve the first mode, which has extraneously low modal participation due to the excitation location. Instead, the identification is biased by the second mode which is observed in Fig. 3. Although the POD appears to identify the first mode better than the SOD and ERA algorithms, the identified mode shape shown in Fig. 3 still does not correspond with the FEM mode shape. Furthermore, the errors for the other three modes are orders of magnitude higher than those identified by the proposed SBD. A train of impulses across the span of the beam may help better excite the first mode for SOD to identify; however, excessive excitation may not always be feasible, especially for large or complicated mechanical systems. Furthermore, the proposed SBD is still expected to work for this situation and therefore has an advantage over the other considered conventional OMA methods.

Two forced vibration examples were considered to understand how the proposed method outperforms others under operational conditions. For the case of single harmonic excitation, the bottom-right subplot in Fig. 4 shows that SBD can accurately identify the first ten modes of the system despite being excited with a high-frequency input signal. SOD and ERA do not share this ability for the first two modes, as shown in Fig. 5. When the frequency is within the bandwidth of the modes of interest, the proposed criteria cannot identify modes higher than the fifth order. This is demonstrated in Fig. 4. In practice, a harmonic excitation with a frequency above the bandwidth of interest (i.e., half of the sampling frequency) can maximize the ability to identify the modes. In this way, all of the modes that are below the effective bandwidth will be appropriately identified. It should be noted that despite the contamination near the excitation band in the criterion shown in Fig. 4 for the 25, 75, and 125 Hz cases, the lower modes are still correctly identified by the criterion. The SBD algorithm generally outperforms both SOD and ERA aside from a few of the identified modes, where only minor differences across the methods can be observed. The ability for SBD to accurately identify the majority of the LNMs using only a single harmonic excitation is valuable for engineering applications where complicated forcing techniques are unavailable or prohibited. Additionally, the proposed SBD algorithm provides reliable damped natural frequency estimates, as shown in Table 4.

Being another popular branch of modal testing method, the non-stationary excitation poses challenges to the considered OMA candidates. For ERA, the non-stationarity of the forcing violates its basic assumptions. Regarding SOD, such non-stationarity makes the lower-order modes less energetic; as a result, SOD cannot resolve these modes with extraneously low modal energy. Logarithmic, exponential, or slow linear sweep signals are often employed to overcome this problem; however, the resulting large data set slows down identification significantly. In contrast, the SBD is insensitive to this constraint. As long as an approximately singular modal structure exists in a particular subband, the SBD algorithm can identify it regardless of the type of excitation signal, as shown by the clustering measure in Fig. 6. Figure 5 shows that the proposed SBD has the smallest absolute error attributed to it aside from a few cases from the third to the sixth modes with RMSE shown in Table 6.

The accelerometers used in the experimental hammer tests are biased against lower modes due to their nonlinearity. Fig. 7(b) shows the PSD of the accelerometer response subject to a bandlimited white Gaussian noise (0–500 Hz). The 0–50 Hz region response experiences nonlinear sensor distortion. As a result, LNMs in this region will have lower energy attributed to their modal motion despite being physically energetic (Table 7). Thus, this represents an example of sensor-induced low modal participation for the low-frequency modes. Despite this unfavorable condition, the SBD algorithm can still identify all four modes within range of the results obtained from LSRF, which was used as the ground truth in this example. The same is not valid for SOD, whose first identified mode appears to be a modulation of the first and second actual LNMs of the cantilever beam, as shown in Fig. 8.

Average of RMSE under swept-sine excitation | ||||||
---|---|---|---|---|---|---|

Mode # | $E\varphi $ | $E\omega $ | ||||

SOD | ERA | SBD | SOD | ERA | SBD | |

All | 5.39 × 10^{–2} | 6.42 × 10^{–2} | 4.0 × 10^{–3} | 7.7 × 10^{–1} | 2.3 × 10^{–1} | 3.0 × 10^{–2} |

Average of RMSE under swept-sine excitation | ||||||
---|---|---|---|---|---|---|

Mode # | $E\varphi $ | $E\omega $ | ||||

SOD | ERA | SBD | SOD | ERA | SBD | |

All | 5.39 × 10^{–2} | 6.42 × 10^{–2} | 4.0 × 10^{–3} | 7.7 × 10^{–1} | 2.3 × 10^{–1} | 3.0 × 10^{–2} |

## 7 Conclusion

A novel subband decomposition (SBD) method is proposed to identify modal parameters from response data saddled with low modal participation during operational modal analysis (OMA). This method first decomposes the raw response trajectory matrix into a series of subband response matrices via an FIR uniform analysis filter bank, where each subband matrix is bandlimited. With the in-band frequency-isolated trajectory matrices, the non-orthogonality of the smooth modes (SMs) is exploited to identify the subbands that are close to resonances where the system response is close to singular. This close-to-singular nature of the subband response motivates a new modal clustering measure, developed to robustly identify the resonant subbands without any prior knowledge of the underlying modes. Once the modal frequency or subband is identified, the most-energetic basis of the subspace spanned by the subband SMs is extracted as the corresponding LNM. Compared to various popular OMA methods, SBD successfully identifies the lower energy LNMs even if the subband response is contaminated by energy leakage from the higher-energy out-of-band modes.

Numerical experiments of a cantilever beam are considered to examine the modal identification accuracy of SBD when the collected response is biased by low modal participation under various stationary and non-stationary loading conditions. For identification accuracy, the results show as much as a two-orders-of-magnitude reduction in SBD estimation error compared to the conventional OMA methods. These case studies also validate the applicability and versatility of the proposed SBD in applications that are limited to narrow-band and non-stationary excitation techniques. Furthermore, to validate the applicability of the SBD in modal testing, physical experiments on a cantilever beam are conducted. The experimental results show a significant advantage of the SBD-based identification when low-frequency modes are subject to sensor nonlinearity-induced low modal participation. The experimental results confirm the advantage of the proposed SBD over various popular OMA methods in identifying LNMs with low modal energy.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Nomenclature

*c*=linear sinusoidal sweep rate

*j*=imaginary unit or matrix row index

*m*=number of rows of the trajectory matrix

*n*=discrete-time variable

*p*=number of columns in the trajectory matrix

*A*=stopband attenuation

*C*=diagonal matrix whose squared values,

*c*^{2}, are numerators of SOVs_{i}*D*=system dynamic matrix in ERA

*E*=root-mean-square error between the estimated modes and the FEM-resolved modes

*H*=Hankel matrix in ERA

*S*=diagonal matrix whose squared,

*s*^{2}, values are denominators of SOVs_{i}*T*=total chirp time during the chirp excitation

*U*=proper orthogonal coordinate matrix

*V*=proper orthogonal mode matrix

*X*=system trajectory matrix

- $A$ =
discrete-time state transition matrix

- $D$ =
temporal derivative operator

*u*_{0}=initial state vector in ERA

*w*=_{i}scalar modal clustering measure

*Y*=_{i}*i*th subband trajectory matrix- $f^i$ =
SBD-estimated damped natural frequency of the

*i*th mode*I*_{0}=zero-order modified Bessel function of the first kind

*S*=_{jk}cosine similarity matrix whose (

*j,k*) nomenclature characterizes the cosine similarity between the*j*th and the*k*th smooth modes- (·)
_{H}_{0}, (·)_{H}_{0r}= the resultant SVD-decomposed matrices and their reduced forms in ERA

- (·)
^{*}= complex conjugation

- (·)
^{†}= Moore–Penrose pseudo inverse

- $D^,P^,Q^,p^i$ =
ERA estimated system dynamics matrix, observability matrix, controllability matrix, and poles

*h*[_{D}*n*] =ideal discrete-time impulse response function

*h*[_{i}*n*] =discrete-time impulse response function of the

*i*th subband*w*[*n*] =discrete-time windowing function

*x*,$x^i$ =_{i}FEM mode amplitude and estimated mode amplitude of the

*i*th mode*H*(*jω*)*, R*=_{l}, p_{l}system transfer function, system residues, and system poles

*BW*=bandwidth of the bandpass filters

*O, P, Q*=observation matrix, observability matrix, and controllability matrix in ERA

### Greek Symbols

*α*_{1},*α*_{2}=Rayleigh damping coefficients

*β*=ripple control parameter

*ζ*=_{i}damping ratio of the

*i*th mode*Θ*=_{jk}angular similarity matrix derived from the cosine similarity matrix

*S*_{ϕjk}*λ*=_{i}square of the

*i*th proper orthogonal value*µ*=_{i}the

*i*th smooth orthogonal value*σ*=_{i}the

*i*th proper orthogonal value*Σ*=proper orthogonal values matrix

*Φ*=_{i}*i*th subband smooth mode matrix*Ψ*,*Ψ*′ =the eigenvectors of the estimated system dynamics matrix in physical and observation space, respectively

*ω*=angular frequency