Frequency domain decomposition

Recovering a shape from vibration measurements

Two sensors on a vibrating structure record two time histories. We want to identify a mode shape: the relative amplitudes and phases with which the sensor locations move in one mode. The forces are not measured. What information in the recordings could reveal that shape?

First consider an ideal case in which only one mode contributes. Suppose that, whenever the first sensor moves 1 mm in its positive direction, the second moves 2 mm in its negative direction. Let $t$ denote time, $a(t)$ the first displacement, and $y(t)$ the column vector of both measurements. The fixed vector $v$ describes their relative motion:

The time dependence is entirely in $a(t)$. It can be irregular; the spatial ratio remains $1:-2$. In this ideal case, inspecting the two time histories would already reveal the shape. Real measurements contain several modes at once, so that ratio usually changes with time. FDD looks for a fixed pattern within a frequency band, where one mode may dominate.

To see what survives in the frequency domain, let $f$ denote frequency in hertz, $A(f)$ the Fourier coefficient of $a(t)$ over a record, and $Y_1(f),Y_2(f)$ the corresponding sensor coefficients. Linearity of the Fourier transform gives $Y_1=A$ and $Y_2=-2A$. An overbar denotes complex conjugation and vertical bars denote magnitude. At each frequency,

The separate power spectra preserve the amplitude ratio, but lose the sign: squaring gives the same result for $2A$ and $-2A$. The cross product preserves their opposite motion. This is why we retain the products between channels.

Averaging these products and applying the spectral normalization gives the power spectral density (PSD) matrix $G_y(f)$. Let $S_a(f)$ be the PSD of the first signal, its mean-square displacement per unit frequency, and let the superscript $T$ denote transpose. For zero-mean stationary records—whose mean is zero and whose expected products depend only on the time separation—the spectral definition below gives

Where $S_a(f)>0$, this matrix has only one independent column, so its rank is one. Multiplication shows that its eigenvector with nonzero eigenvalue is precisely the measured shape:

The eigenvalue depends on how strongly the motion is excited. The eigenvector gives the relative motion $1:-2$, regardless of that strength. For a rank-one PSD matrix, a singular value decomposition recovers this same direction.

Frequency domain decomposition (FDD) applies this calculation to the measured spectral matrix at each frequency. The remaining task is to establish when a mechanical mode produces an approximately rank-one spectrum, and when other modes or the excitation invalidate that interpretation. [3]

Modes of a mechanical system

Consider a linear mechanical system with $n$ displacement coordinates, $r$ force inputs, and $p$ measured channels. A coordinate is one degree of freedom; a channel is one recorded sensor signal. Let $q(t)$, $x(t)$, and $y(t)$ denote the displacement, force, and measured-displacement vectors.

The real matrices $M,C,K$ are the $n\times n$ mass, viscous damping, and stiffness matrices, with $M$ invertible. The matrices $B$ and $L$, of sizes $n\times r$ and $p\times n$, distribute the forces and select the measurements. A dot denotes a time derivative. The equations are

Linearity means that responses add and scale with the inputs; time invariance means that the matrices are constant. The distinction between $q$ and $y$ matters: a structure may have many coordinates and only a few sensors.

For free motion, set $x=0$. Physical displacement $q(t)$ is real. To find the modes, we temporarily allow a complex auxiliary solution, denoted $q_c(t)$. We seek a fixed spatial vector $\phi_m\ne0$, the shape of mode $m$, evolving at a complex rate $\lambda_m$ in inverse seconds:

Substituting into (4) and cancelling the nonzero exponential leaves a matrix polynomial. For a complex rate $s$, define the dynamic stiffness matrix $Z(s)=s^2M+sC+K$. Writing $\det$ for the determinant, the admissible rates and shapes satisfy the quadratic eigenvalue problem:

The last equality follows because $Z(\lambda_m)$ maps a nonzero vector to zero and is therefore singular. A numerical solution normally uses an equivalent first-order eigenproblem; the state-space derivation is given below. [1]

Solving (6) gives the rates $\lambda_m$. To interpret one of them, return to the time factor $e^{\lambda_m t}$ in (5). Let $\sigma_m$ and $\omega_{d,m}$ be its real and imaginary parts, and let $j$ be the imaginary unit, $j^2=-1$. Euler’s formula separates the exponential into an amplitude factor and an oscillation:

If $\sigma_m<0$, the amplitude decays; if $\sigma_m>0$, it grows. A nonzero imaginary part produces oscillation. Since the system matrices are real, nonreal eigenvalues occur in conjugate pairs. Select the member with $\omega_{d,m}>0$. For a decaying mode, define $\alpha_m=-\sigma_m>0$, so $\lambda_m=-\alpha_m+j\omega_{d,m}$. The decay rate is therefore $\alpha_m$, and the damped angular frequency is $\omega_{d,m}$, in radians per second. The complex solution is a way to represent this motion; it is not a complex physical displacement.

Because the free equations are linear with real coefficients, the real and imaginary parts of any complex solution separately satisfy them. Equivalently, a conjugate pair combines into a real modal contribution. Let $c_m$ be a complex coefficient fixed by the initial conditions, $q_m(t)$ the displacement contribution of this pair, and $\operatorname{Re}$ denote real part:

The imaginary parts cancel. A complex mode shape records relative phases between coordinates; each coordinate still has a real displacement. Summing these real contributions gives the physical free response $q(t)$, with the coefficients chosen to match its real initial displacement and velocity.

Natural frequency and damping ratio

We already know the decay rate and oscillation frequency from the eigenvalue. The terms natural frequency and damping ratio come from a mass–spring–damper system. Consider a mass $\mu>0$, spring stiffness $k>0$, damping coefficient $b\ge0$, and displacement $u(t)$ from equilibrium. With no external force,

First remove the damper, setting $b=0$. Dividing by the mass gives $\ddot u+(k/\mu)u=0$. A cosine at angular frequency $\omega_n$ has second derivative $-\omega_n^2$ times itself, so this equation requires $\omega_n^2=k/\mu$. This is the undamped natural angular frequency: the rate of free oscillation set by mass and stiffness. Its value in cycles per second is $f_n=\omega_n/(2\pi)$.

Restore the damper. Before defining a damping ratio, find the value of $b$ at which oscillation disappears. Substituting $e^{st}$ into the equation of motion and cancelling the exponential gives the characteristic equation: its roots are the allowed exponential rates.

The expression under the square root determines the type of motion. When $b^2<4\mu k$, the roots are complex and produce oscillation; when $b^2>4\mu k$, both roots are real and negative, so the motion contains only decaying exponentials. At the boundary, the roots coincide. Denote this critical damping coefficient by $b_{\mathrm{crit}}$. Since $b\ge0$,

The symbol $\equiv$ here introduces a definition. We define the dimensionless damping ratio $\zeta$ as the actual damping divided by this critical value. Thus $\zeta=0.05$ means 5% of critical damping; $0<\zeta<1$ is underdamped, $\zeta=1$ is critically damped, and $\zeta>1$ is overdamped. The denominator and its factor of 2 come directly from the quadratic roots. [9]

Using $\omega_n=\sqrt{k/\mu}$, the damping coefficient per unit mass becomes $b/\mu=2\zeta\omega_n$. Dividing the motion equation by $\mu$ therefore gives the conventional oscillator equation and its characteristic polynomial:

Returning to the modal eigenvalue

The link between an eigenvalue and a polynomial root comes from the exponential solution. In the mechanical system, substituting $q_c(t)=\phi_m e^{\lambda_m t}$ gave $Z(\lambda_m)\phi_m=0$. Since $\phi_m\ne0$, the matrix is singular: $\det Z(\lambda_m)=0$. Thus $\lambda_m$ is both an eigenvalue of the mechanical problem and a root of its determinant polynomial.

For a scalar oscillator the same substitution gives a scalar equation. Denote its characteristic polynomial by $p(s)=s^2+2\zeta\omega_ns+\omega_n^2$. To test whether it admits the rate $\lambda_m$, substitute the auxiliary complex motion $u_c(t)=e^{\lambda_m t}$:

The last step follows because the exponential never vanishes. Therefore the oscillator admits $e^{\lambda_m t}$ precisely when $\lambda_m$ is a root of $p$. It is not automatically a root of an arbitrary oscillator’s polynomial. We are choosing an equivalent oscillator whose two exponential rates match the already computed mechanical pair $\lambda_m=-\alpha_m+j\omega_{d,m}$ and $\bar\lambda_m=-\alpha_m-j\omega_{d,m}$.

Its characteristic polynomial must be zero at both rates. The factor $s-\lambda_m$ is zero at the first, and $s-\bar\lambda_m$ is zero at the second. Their product is therefore a quadratic with exactly the required roots. We choose its leading coefficient to be 1, as in the oscillator equation after division by mass. For example, roots 2 and 3 give $(s-2)(s-3)=s^2-5s+6$; the same factorization rule applies to complex roots.

Substitute the two eigenvalues into these factors and expand. Since $j^2=-1$, the imaginary terms cancel:

Matching this with $s^2+2\zeta_m\omega_{n,m}s+\omega_{n,m}^2$ requires $\zeta_m\omega_{n,m}=\alpha_m$ and $\omega_{n,m}^2=\alpha_m^2+\omega_{d,m}^2$. These define the oscillator parameters associated with mode $m$. With $|\lambda_m|$ denoting magnitude,

The eigenvalue now gives three quantities: $\alpha_m$ describes how quickly the amplitude decays, $\omega_{d,m}$ describes how fast it oscillates, and $\zeta_m$ expresses the damping relative to the critical value. The plot below shows how these quantities change together.

How a mode enters the forced response

Free motion identifies the possible rates and shapes. To connect them to measurements under random forcing, we need the transfer matrix.

The Laplace transform converts a differential equation in time into an algebraic equation. Write $s=\sigma+j\omega$, where $\sigma$ is a real exponential weighting rate and $\omega$ is angular frequency. Denote the Laplace transform by $\mathcal L$, and let $Q(s)$ be the transform of displacement, taken entry by entry:

This is a Fourier transform of the exponentially weighted signal on $t\ge0$. Choosing $\sigma$ sufficiently large makes the integral converge for signals that grow no faster than an exponential. The variable $s$ ranges over complex values; the eigenvalues $\lambda_m$ found earlier are particular rates belonging to the system.

The useful property follows by integration by parts. The endpoint term at infinity vanishes in the region of convergence, giving

Thus, with zero initial displacement and velocity, each time derivative becomes multiplication by $s$. [10] We impose these initial conditions to isolate the response caused by the force input; a response due to initial motion can be added separately.

Let $X(s)$ and $Y(s)$ be the transforms of the force $x(t)$ and measured displacement $y(t)$. Applying the transform to (4) gives

The same matrix $Z(s)$ appeared in the free-motion eigenproblem. Wherever it is invertible, solve for $Q(s)$ and substitute into $Y(s)=LQ(s)$. The resulting matrix mapping transformed forces to transformed measurements is the transfer matrix, denoted $H(s)$:

Separating the modal contributions

The transfer matrix contains all the modes at once. To see their individual contributions, we expand it into fractions whose denominators contain one eigenvalue each.

Return to the mass–spring–damper oscillator, now driven by a force $F(t)$. Its displacement is $u(t)$, and its mass, damping coefficient, and stiffness are $\mu,b,k$. Divide by the mass and define the input $g(t)=F(t)/\mu$, the force per unit mass. With the frequency and damping parameters defined above,

Let $U(s)$ and $G(s)$ be the Laplace transforms of $u(t)$ and $g(t)$. With zero initial displacement and velocity, transforming gives

Here $h(s)$ maps force per unit mass to displacement. Its numerator is 1 because the normalized input appears with coefficient 1. For a transfer function from the actual force $F$, the numerator would instead be $1/\mu$.

The denominator is the same characteristic polynomial that determines the free oscillation. For an underdamped oscillator, its two roots are $\lambda$ and $\bar\lambda$. A quadratic with leading coefficient 1 factors as the product of $s$ minus each root, so $s^2+2\zeta\omega_ns+\omega_n^2=(s-\lambda)(s-\bar\lambda)$. Partial fractions then give

As $s$ approaches $\lambda$, the first term becomes unbounded while the second stays finite. We call $\lambda$ a pole and the coefficient $1/(\lambda-\bar\lambda)$ its residue. The pole locates the singularity; the residue specifies its contribution. This decomposition also connects to time: the Laplace transform of $e^{\lambda t}$ is $1/(s-\lambda)$, wherever the integral converges. Each fraction therefore corresponds to one of the exponential rates found in the free-motion problem.

Now consider one entry of the transfer matrix. Let $H_{ik}(s)$ map force input $k$ to sensor $i$. Each entry is a ratio of polynomials: the inverse in $H=LZ^{-1}B$ has denominator $\det Z(s)$, although factors can cancel. Its partial fractions therefore have the mechanical eigenvalues as their possible poles.

Assume the $2n$ roots of $\det Z(s)$ are distinct and form decaying, oscillatory conjugate pairs. Choose one root, $\lambda_m$. In the expansion of $H_{ik}$, let $r_{ik,m}$ be the coefficient of $1/(s-\lambda_m)$, and let $g_{ik,m}(s)$ collect the fractions associated with all the other poles, including the conjugate pole $\bar\lambda_m$. Only the term for $\lambda_m$ is separated out. Then

This is an exact separation of the sum. The remainder $g_{ik,m}$ stays finite as $s$ approaches $\lambda_m$, because its denominators vanish at other, distinct roots. To recover the coefficient $r_{ik,m}$, multiply by $s-\lambda_m$ and take the limit:

The last term tends to zero: it is a vanishing factor times a bounded function. This is why the limit isolates the residue.

Repeat this for every sensor–input pair. Arrange the constants $r_{ik,m}$ into a matrix $R_m$, and the remainder functions $g_{ik,m}$ into a matrix $H_{\mathrm{rest},m}$, keeping their row and column positions. The entrywise equations become

Thus $R_m$ is a constant $p\times r$ matrix, with one coefficient for each of the $p$ sensors and $r$ inputs. It describes only the chosen pole’s contribution; $H_{\mathrm{rest},m}$ contains the others, including the conjugate pole. A zero coefficient means that this particular entry has no contribution from the chosen pole.

Why the residue contains the mode shape

First distinguish the physical shape from the measured shape. The vector $\phi_m$ describes mode $m$ at all $n$ displacement coordinates. The matrix $L$ maps those coordinates to the $p$ sensor readings, so define $v_m=L\phi_m$. For example, if $\phi_m=[1,3,-2]^T$ and the sensors measure only coordinates 1 and 3, then $v_m=[1,-2]^T$. The sensors see only that part of the shape.

Why must different force inputs give multiples of this same vector? Let $T(s)=Z(s)^{-1}B$ be the transfer matrix from forces to all physical displacements, so $Q=TX$. The transfer matrix $H$ from (8) maps forces to the sensor readings: since $Y=LQ=LTX$, we have $H=LT$. Denote the residue of $T$ at $\lambda_m$ by $W_m$. From $Z(s)T(s)=B$, multiplying by $s-\lambda_m$ and taking the limit gives

Thus every column of $W_m$ solves the same equation as the mode shape: multiplication by $Z(\lambda_m)$ gives zero. At the assumed simple root, this nullspace has dimension one, so every such column must be a scalar multiple of $\phi_m$. Let $\gamma_{m,k}$ be that scalar for force input $k$. The corresponding sensor column is therefore $L(\phi_m\gamma_{m,k})=v_m\gamma_{m,k}$.

Now use $H=LT$ to return to the measured response. Because $L$ is constant, it can be taken outside the residue limit:

Thus the result for the physical residue $W_m$ gives the sensor residue $R_m$, which is the quantity needed in the expansion of $H$.

For example, if two inputs have coefficients 2 and $-3$, and the measured shape is $v_m=[1,-2]^T$, their residue columns are

The columns differ in strength and sign, but both retain the sensor ratio $1:-2$. This statement concerns the residue of one pole; the columns of the full transfer matrix need not be proportional because they include other modes.

Keep the two-input example and let $X_1(s),X_2(s)$ be the transforms of the applied forces, so $X(s)=[X_1(s),X_2(s)]^T$. Multiplying the residue by this input vector combines its two columns:

Both forces affect the response through the single quantity $2X_1-3X_2$. The first sensor receives that quantity and the second receives $-2$ times it. This is the fixed spatial pattern.

The general notation expresses exactly this calculation. For $r$ inputs, collect their coefficients into the row $\gamma_m^T=[\gamma_{m,1},\ldots,\gamma_{m,r}]$. Factoring the common vector out of the residue columns gives

Here $X_k(s)$ is the transform of input $k$. The product $\gamma_m^TX$ is one scalar because it is a row times a column. In the example, $\gamma_m^T=[2,-3]$. These coefficients belong to the system; $X$ describes the forces actually applied. The appendix derives the coefficients explicitly.

Recall that the residue is only the numerator of the pole term. Its contribution to $H(s)$ is $R_m/(s-\lambda_m)$. Since $Y=HX$, the contribution of this one term to the transformed sensor response is

The denominator changes the frequency dependence, but it divides every sensor entry by the same number. The relative sensor amplitudes and phases remain those of $v_m$. A nonzero residue has rank one because its columns are all multiples of this vector. [1] [3]

Adding the other poles

So far we have kept only the term for $\lambda_m$. The full transfer matrix also includes its conjugate pole and every other mode. For real system matrices, the residue at $\bar\lambda_m$ is $\bar R_m$. Under our assumption of $2n$ distinct roots in conjugate pairs, summing the $n$ pairs reconstructs $H$:

The index $m$ runs over the pairs, choosing the eigenvalue with positive imaginary part in each pair. This is the full partial-fraction expansion, not an assumption that one mode dominates.

A mode can be absent from the transfer matrix: if $v_m=0$, the sensors do not observe it; if $\gamma_m^T=0$, none of the modeled inputs excites it. In either case $R_m=0$, so that root contributes no pole to $H$. [2] Even a visible mode may receive no contribution from a particular force combination for which $\gamma_m^TX=0$.

For example, suppose a mode has shape $\phi_m=[1,0,-1]^T$: the first and third coordinates move in opposite directions, while the middle coordinate stays still. If the only sensor measures that middle coordinate, then $L=[0,1,0]$, giving $v_m=L\phi_m=0$ and hence $R_m=0$. The mode can be vibrating elsewhere, but this sensor records none of its motion. Moving the sensor to the first coordinate makes this mode visible.

Equation (9) separates the transfer matrix into modal terms. The measured response still contains their sum. We next examine when one term dominates its spectral matrix strongly enough for FDD to recover the shape.

The response spectral matrix

We have expressed the forced response as a sum of modal contributions. In an operational measurement, however, the forces are unknown and fluctuate from record to record. We therefore need a description computed from the sensor signals themselves. The spectral matrix introduced in the opening example does this by averaging products of their Fourier coefficients.

From sensor records to a spectral matrix

Let $T$ be the duration of a record and $Y_T(f)$ the vector obtained by Fourier-transforming each sensor signal at frequency $f$. At one frequency, abbreviate this vector by $Y$. For two sensors its entries are $Y_1,Y_2$. The superscript $H$ means conjugate transpose: transpose the vector and conjugate its entries. It is an operation, distinct from the transfer matrix $H(s)$. The outer product is

The diagonal contains squared amplitudes. The off-diagonal entries compare the channels. For channels $i$ and $k$, let $\theta_i,\theta_k$ be their Fourier phases. Then $Y_i\bar Y_k=|Y_i||Y_k|e^{j(\theta_i-\theta_k)}$: conjugation leaves the phase difference. This is the information that distinguishes in-phase motion from opposite motion in our two-sensor example.

For random excitation, the Fourier coefficients vary between records. Write $\mathbb E$ for expectation, an average over possible records generated under the same statistical conditions. Averaging $YY^H$ retains the cross-channel products that persist across those records; unrelated phase fluctuations can cancel.

We assume zero-mean, wide-sense stationary measurements: each signal has zero mean and finite mean square, and its expected product with another channel depends only on the separation between the two times. Assume also that a spectral density exists in the band of interest. The response power spectral density (PSD) matrix $G_y(f)$ is

The factor $1/T$ converts the squared-transform energy density of a record into a mean-square density. “Two-sided” means that both positive and negative frequencies are included. Integrating a diagonal entry over all frequencies gives that channel’s mean-square displacement, which equals its variance because the mean is zero. A displacement PSD therefore has units of displacement squared per hertz. [5]

In practice, we have finite data, so we estimate this average using windowed segments of a recording. The estimation appendix gives the normalization and averaging procedure.

Why this matrix admits the FDD decomposition

The outer product has two useful properties. Swapping its row and column indices conjugates the entry, so $YY^H$ equals its conjugate transpose. Such a matrix is called Hermitian. Averaging preserves this property.

Next, choose any complex column vector $a$ of sensor weights. The scalar $a^HY_T$ is a weighted combination of the Fourier coefficients. Its squared magnitude cannot be negative, and neither can its average:

This second property is called positive semidefiniteness: every weighted combination has a nonnegative PSD. Together, the two properties ensure that $G_y$ has real, nonnegative eigenvalues and an orthonormal eigenvector basis. This is the matrix structure used by FDD.

Connecting the spectrum to the mechanical system

The measured matrix $G_y$ describes the response. To interpret its directions as mode shapes, we must connect it to the transfer matrix derived earlier. For a stable system in stationary operation, its frequency response is $H(j2\pi f)$. Define $G_x(f)$ as the force PSD matrix, using the same Fourier-product convention as (10).

The input–output relation multiplies the force Fourier components by $H$. Consequently, their outer products transform by multiplication on both sides. Suppressing the frequency arguments and writing $X,Y$ for input and output Fourier components, the algebra is

At a fixed frequency, $H$ is determined by the system and does not vary between realizations, so it can be taken outside the expectation. Applying the spectral normalization gives the exact relation

The force spectrum determines how the system is excited; the transfer matrix contains its modal dynamics. [3] We now substitute the modal expansion of $H$ to determine when $G_y$ is dominated by one measured shape.

When the spectrum is approximately rank one

Suppose the force PSD varies little across a resonance. Denote its approximately constant value by $Q_x$. A constant PSD is called white; the matrix need not be diagonal, so different force channels may still be correlated.

Write $\omega=2\pi f$ for angular frequency and let $E_m(j\omega)$ contain every transfer contribution except the selected pole, including its conjugate partner. Then

The selected term becomes large when $\omega$ is close to $\omega_{d,m}$ and $\alpha_m$ is small. If it dominates the response to the actual excitation, we can omit $E_m$ in (12). Let $d_m(f)$ denote the resulting scalar spectral weight and $\approx$ denote approximate equality:

The numerator is nonnegative: with $b=\bar\gamma_m$, it is $b^HQ_xb$, a quadratic form of the positive-semidefinite force PSD. The frequency dependence is therefore confined to $d_m$; the direction $v_m$ is fixed. This is the general version of the two-sensor example.

The assumption concerns the excited response, not damping alone. A nearby mode can remain important, and a lightly damped mode may be weakly forced or poorly observed. The residual bound below makes this distinction precise. [3]

Extracting the mode shape

In practice, $G_y$ is replaced by an estimate $\widehat G_y$, with the hat indicating a quantity computed from data. FDD applies a singular value decomposition (SVD) separately at each frequency. For a Hermitian positive-semidefinite matrix, the singular vectors can be chosen as unit eigenvectors with zero pairwise inner products (orthonormal eigenvectors). The singular values are the corresponding eigenvalues.

At frequency $f_i$, let $u_k(f_i)$ and $s_k(f_i)$ be these vectors and values, ordered largest first. Assemble the vectors into the columns of $U_i$ and the values into the diagonal matrix $S_i$. Then

For a vector $v$, write $\|v\|=\sqrt{v^Hv}$ for its Euclidean norm. The rank-one matrix in (14) has a single nonzero eigenvalue, since

Let $\theta$ denote an arbitrary common phase angle. Normalizing the eigenvector gives

The first singular value describes the spectral strength of the mode; the first singular vector estimates its relative amplitudes and phases at the sensors. For the opening example, these are $5S_a(f)$ and $[1,-2]^T/\sqrt5$. The same calculation holds with any number of measured channels. [3]

Neither the absolute shape scale nor the force amplitude is identified. Multiplying the shape by a constant and compensating in the spectral weight gives the same matrix. Consequently, output-only FDD does not determine modal mass, the effective mass associated with a chosen shape normalization.

Peak frequency is not pole frequency

Let $f_{\mathrm{peak}}$ denote the frequency of a spectral maximum and $\omega_{\mathrm{peak}}=2\pi f_{\mathrm{peak}}$. Even a single-degree-of-freedom (SDOF) oscillator distinguishes this frequency from $\omega_n$ and $\omega_d$. With $h(j\omega)$ denoting its force-to-displacement transfer function, and modal subscripts omitted, its squared magnitude is given below. The symbol $\propto$ suppresses a frequency-independent factor; under white forcing, the displacement PSD has this same frequency dependence:

Differentiating the denominator gives $4\omega[\omega^2-\omega_n^2+2\zeta^2\omega_n^2]$; its nonzero root yields the stated maximum when $\zeta<1/\sqrt2$. At light damping the three frequencies are close, which explains the usefulness of peak picking. Colored forcing, acceleration measurements, and finite spectral resolution can shift a measured maximum further.

Two modes in the same band

Consider two independent modal responses: knowing one leaves the probability distribution of the other unchanged. Let their measured shapes be $v_1,v_2$ and their spectral weights $d_1,d_2$. Their contributions add:

Multiplying by the first physical shape gives

The second term generally points away from $v_1$. Therefore a singular vector need not equal either physical shape. If both weights are positive and the shapes are linearly independent, the leading two vectors recover their span—the set of their linear combinations. Recovering each mode separately requires more than a rank-two decomposition. [4]

Orthogonal measured shapes, satisfying $v_1^Hv_2=0$, are a special case. They should not be confused with mass-orthogonal physical modes, satisfying $\phi_1^HM\phi_2=0$. The latter does not imply the former. Equal singular values also leave the basis within the corresponding subspace undetermined.

Estimating damping

Enhanced frequency domain decomposition (EFDD) uses the spatial estimate to select the spectrum associated with one mode, then transforms that spectrum into a correlation decay. The selected part is called a spectral bell. Its width contains damping information that a peak location alone cannot supply. [4]

The modal assurance criterion (MAC) measures agreement between two nonzero vectors. Let $\hat v_m$ be a reference shape taken near a resonance and $u_k$ a candidate singular vector at a neighboring frequency. Their MAC is

The Cauchy–Schwarz inequality bounds an inner-product magnitude by the product of the norms, so MAC lies between 0 and 1. It is 1 for proportional vectors and 0 for orthogonal vectors; a common complex scale has no effect.

Starting at the reference peak, follow the matching spatial direction outward on both sides. This direction can move between ordered singular vectors when singular values cross. A MAC threshold of 0.8 is a practical starting value, not a universal constant. Stop when the match fails, set unselected bins to zero, and inspect the resulting bell for another mode’s contribution. [7]

A superscript $+$ denotes a one-sided PSD, containing only nonnegative frequencies with the positive interior contributions doubled. Let $b^+(f)$ be the selected bell and $\mathcal R_m(\tau)$ its autocorrelation at lag $\tau$, the time separation between two response values. For a discrete calculation, let $f_s$ be the sampling rate, $N$ an even transform length, and $\Delta f=f_s/N$ the frequency spacing. Write $b^+[k]$ for the sample at frequency $k\Delta f$, and $\ell$ for the lag index. Then

The cosine transform combines equal positive- and negative-frequency contributions. Its discrete implementation requires the one-sided scaling described under spectral estimation.

For an isolated mode under the EFDD approximation, the correlation oscillates with a decaying envelope. Let $a_m(t)$ be that positive envelope, $A_0>0$ its extrapolated initial amplitude, and $T_d$ the full damped period. Here $t$ measures elapsed lag. The logarithmic decrement $\delta$ is the natural logarithm, denoted $\ln$, of the amplitude ratio per period. For peaks separated by $k$ full periods,

The exponential makes the logarithmic ratio independent of the starting time. Define $f_{d,m}=\omega_{d,m}/(2\pi)$. Substituting (7) gives $\delta=2\pi\zeta_m/\sqrt{1-\zeta_m^2}$, and hence

A fit to several log peak amplitudes is more useful than one ratio. Same-sign peaks are a full period apart; alternating extrema are half a period apart. The fit should exclude initial distortion, late noise, and beating, the slow amplitude modulation from nearby modes. Bell truncation and spectral leakage—the spreading caused by a finite window—can bias the decay, so the damping estimate should be checked across reasonable spectral and fitting settings. [4] [7] [8]

Using the method on measurements

The calculation is short: estimate the full cross-spectral matrix, decompose it at every frequency, inspect the dominant vectors around candidate peaks, and, if damping is needed, select and fit the corresponding correlation decay. The interpretation depends on a few checks.

• Measurement consistency. Channels must be synchronized, with known units and sensor directions. Arbitrary channel rescaling changes the spatial geometry seen by the SVD.
• Spectral estimation. Use the same segments and window in every channel pair. Too few averages can produce an artificially small matrix rank. Zero-padding gives a denser frequency grid, not finer resolution.
• Modal interpretation. A peak with a stable spatial direction is a candidate mode. A sinusoidal force can also produce a rank-one peak; the decomposition alone cannot distinguish the two.
• Reported quantities. State the sensor coordinates, shape normalization, cross-spectrum convention, and spectral settings. For EFDD, include the selected bell and the fitted lag interval.

An absent peak does not establish an absent mode. The relevant shape may be weakly excited or invisible at the available sensors. These limitations follow from the same input and measurement factors that appear in the transfer matrix. [3]