Digital Image Correlation Vibrometry with Low Speed Equipment

A low-cost method is presented which enables digital image correlation (DIC) with conventional cameras (i.e. not high-speed) to be used for determination of vibration deflected shapes via the use of a stroboscopic lamp and some simple ancillary circuits. For each natural frequency of the structure under consideration, a sequence of images is captured asynchronously with the vibrations using the DIC system and the resulting displacement fields are correlated with the excitation signal driving the vibration using a least-squares approach. Three approaches for performing this correlation are outlined, one of which is developed into the algorithm used for processing the present results to obtain the amplitude and phase of the vibration at each point on the specimen, allowing the deflected shape to be reconstructed. This process is illustrated using the example of a vibrating aluminium plate. The resulting shapes and frequencies agree well with finite element modal analyses of the plate.


Introduction
This work describes a method of obtaining displacements (and, potentially, strains) at high-speed and high resolution using a digital image correlation (DIC) system without the need for high-speed cameras. It is applied to the determination of vibration mode shapes (strictly, deflected shapes as approximations to mode shapes) as an alternative to laser vibrometry, and builds upon existing work which fits DIC data to assumed models of strain or displacement behaviour in order to minimise the effects of noise and improve the sensitivity of the DIC system.
The experimental determination of vibration mode shapes has involved a range of techniques over many years, ranging from the generation of Chaldni figures (where powder granules congregate on nodal lines) via the use of impact modal analysis (for example, multiple-input single-output tests) to the use of scanning laser Doppler vibrometry. Of these, only scanning laser Doppler vibrometry can achieve fine full-field measurement of vibration and unlike digital image correlation it cannot measure all points simultaneously. Conversely DIC has become a popular method for full-field experimental determination of displacement and strain, including the use of 3D DIC techniques which can capture out-of-plane displacements and eliminate the effects of such displacement from the strain field. Recent years have seen the use of DIC in conjunction with high-speed cameras to capture the mode shapes of vibrating objects via a continuously-captured sequence of images. For example, Helfrick et al. [1] compare and contrast modal analysis via high-speed DIC with impact tests and laser vibrometry, while Wang et al. [2] use high-speed DIC in conjunction with random excitation as the basis of establishing modal parameters. Beberniss et al. [3] describe the use of DIC with high-speed cameras to analyse the response to random excitation. An overview of non-contact modal analysis methods, including the use of high-speed DIC, is presented by Avitabile et al. [4]. The application of a Dantec Q-450 high-speed DIC system to operational modal analysis is described by Trebuňa et al. [5], focussing on the additional software to perform the necessary post-processing of the DIC data.
Both laser Doppler vibrometry for modal analysis, and DIC in conjunction with high-speed cameras, require capitalintensive equipment. Moreover, the displacements associated with vibration are often small, meaning that direct DIC measurements of vibration mode shapes will tend to be noisy, and high-speed cameras with low resolutions would probably be unable to resolve satisfactorily the small displacements encountered. The present approach aims to extend the capabilities of a high resolution DIC system to capture mode shape data with minimal outlay. The apparatus consists of a 3D DIC system with conventional cameras rather than high-speed cameras. Images are captured asynchronously with the vibrations and used to determine 3D displacements via image correlation in the usual manner. Least-squares fitting techniques are used to exploit the relationship between the measured displacements and the excitation signal driving the vibrations. These approaches build upon work carried out by Lu et al. [6,7] to determine elastic strain field distributions in composites, in which the elastic model was fitted to the strain results by performing a linear regression of strain against the load signal. The fitting of a linear elastic model has also been reported by other workers [8]. While Lu's work [6,7] involved static rather than dynamic strains or displacements, the underlying issue faced was the same as here: the measurements obtained from the DIC system are subject to significant noise as the values measured are approaching the strain or displacement resolution limit of the DIC system, but by fitting the measurements to an appropriate model such as linear elasticity or a sinusoidal response, the quality of the results can be significantly improved and useful results obtained where direct measurements would be obscured by noise.
There is significant overlap, at least in terms of broad principles, between the approach described here and the lock-in DIC approach developed by Fruehmann et al. [9], who describe the use of a low-speed DIC system to capture cyclic strain data using a lock-in amplifier. As with Fruehmann et al., the images for DIC processing are captured here at a lower frequency than that of the vibration or cyclic loading, resulting in an undersampled set of data. A reference signal related to the loading or excitation is captured to enable the images to be correctly sequenced, and a second reference signal in quadrature with the first is created in software to resolve the ambiguity in placing the images correctly in each vibration cycle. Like the present approach, Fruehmann's method has been applied to vibration up to a few hundred Hertz [10,11], in that case within a wider exploration of damage detection. Unlike that approach however, the generation of the quadrature signal in the present work is performed using hardware (and errors in the π/2 phase shift are corrected in software), a lock-in amplifier is not used, and a xenon stroboscopic lamp (in the form of a high-powered photographic flash unit), in conjunction with very short exposure, is used to Bfreeze^the images rather than using a steady light source. The present work also covers higher frequencies (exceeding 1 kHz) than those reported in references [9][10][11], and may be regarded as being based on asynchronous rather than synchronous lock-in principles [12].
Various other references also show similarities with the present approach but do not duplicate it. For example Serio et al. [13] use an LED-based stroboscopic method to illuminate MEMS at four equally-spaced phases while conducting 3D DIC in order to analyse the vibrations. Helfrick et al. [14] use a stroboscope in conjunction with a 3D DIC system to measure out-of-plane vibrations of rotating blades. Shih and Sung [15] describe the use of a 3D DIC system without highspeed cameras to perform modal analysis on a building subjected to earthquake loading; the frequencies are low (up to around 7.5 Hz) and only limited details are presented regarding the processing of the results.

Development of Analysis Approach
The detailed apparatus and procedure are outlined in BApparatus and procedure^section, but may be summarised as involving a free-running camera and flash system capturing images along with two voltage reference signals related to the excitation. The resulting data set consists of: & A set of (typically) 60 pairs of images captured over a very large number of vibration cycles. These are then processed using the Dantec ISTRA-4D software to determine the displacement at every facet on the model. & A set of samples of analogue voltage signal pairs, approximately in quadrature, captured simultaneously with each image pair.
Three approaches have been explored for undertaking a Bbest fitting^approach for relating the sampled displacement to the excitation voltage.
1. Linear regression of displacement against excitation voltage is a simplistic approach which was sufficient to demonstrate the basic feasibility of the present work in early trials but requires the displacements at every point to be in phase (or antiphase) with the excitation voltage, i.e. the various phase shifts in the transmission and response must add up to a multiple of π radians. Moreover this approach is unable to handle complex modes i.e. those where there is a phase difference (other than π) between motions at different points.

2.
Fitting of an ellipse to the data forming a special case of a Lissajous figure (the pattern formed when two sinusoidal signals are plotted against each other [16]) plotting the displacement at a given point against the excitation voltage, using an existing ellipse-fitting algorithm e.g. that of Fitzgibbon et al. [17,18]. This approach was used initially as it provides a simple method of comparing two signals of equal frequency and unknown phase, and it is possible (in principle) to extract the amplitudes of both the excitation and response signals and their phase relationship within a range of π radians. However, two problems arise from this approach: a. This method alone cannot extract phase angle within a range of 0 to 2π radians as there is an inherent phase ambiguity (positive and negative phase angles result in identical ellipses). One method of overcoming this problem would be to create a second signal which is approximately in quadrature with the original signal (for instance by using an integrating circuit), construct a second Lissajous figure relating displacement to the quadrature signal, then follow a decision tree to extract an unambiguous value of phase angle from parameters of the two fitted ellipse. b. It is not straightforward to take account of the differing amounts of scatter present on the excitation and response signals, noting that in practice the excitation signal is captured with very little noise while the response signal is close to the noise floor of the DIC system due to the very small displacement amplitudes involved. Published ellipse-fitting algorithms e.g. [17] generally aim to minimise the sum of squares of algebraic distances between the points and the fitted curve rather than (for example) minimising the errors in the y-coordinate. This leads to large errors (which vary from point to point) in fitting the data to the excitation signal. Rather than re-deriving the fitting algorithm with the assumption that error exists only in the y signal, an alternative solution was instead pursued as described in (3). However, the Lissajous figure approach was retained for finding the exact phase relationship between two reference signals nominally in quadrature. 3. Direct fitting of sinusoidal terms was the approach finally used and involved two stages: firstly to determine the value of circular time or angular position within a cycle ωt (in the range 0 to 2π radians) corresponding to each image in the sequence taking the excitation signal as the phase reference, then to perform a direct least squares fit of the displacement at each point to a sine wave. This approach inherently assumes that the error exists purely in the displacement terms. While in principle the value of ωt may be found from the arctangent of the directly-measured and integrated excitation signals, it is found in practice that the integrated signal does not lag behind the directly measured signal by precisely π/2 radians and further processing is needed to correct for this phase error before determining the value of ωt for each point.
As noted earlier, the above approach (3) has strong similarities with some implementations of lock-in thermography [12], with that concept providing one of the starting points for the present work, and (consequently) also has similarities with the lock-in DIC method described by Fruehmann [9] and applied to vibration by Waugh et al. [10,11] However, the present approach was derived independently of references [9][10][11] and some major differences were noted in the BIntroductionŝ ection. The full implementation is described in the next two sections.

Apparatus and Procedure
The method was used to measure the deflected vibration s h a p e s o f a n a l u m i n i u m p l a t e o f d i m e n s i o n s 350 mm × 225 mm and thickness 2 mm. The plate was supported on four polyurethane foam pads located at approximately a quarter of each side length in from the edges. Excitation was provided via a loudspeaker driver unit located immediately underneath the plate, driven via a GW-Instek analogue function generator with 75Ω output. The plate displacements were measured using a Dantec Q-400 3D digital image correlation system equipped with two 5 megapixel cameras with a maximum frame rate in the order of 6 fps. The excitation signal was fed into analogue input 1 of the Q-400 system, and was integrated using the circuit in Fig. 1 (based on that in ref. [19]) and fed into analogue input 2. Frequency measurement and monitoring of the signal amplitudes and phases was provided using a GW-Instek GDS-1022 oscilloscope. Lighting was provided using a single unit of a Bowens model ProLite 60 flash lamp system, triggered from the Q400's Bsync^output signal using a simple interface circuit (Fig. 2). The Q400's existing high intensity LED array light source was used only during focussing and calibration, not within the measurement procedure. The apparatus is shown schematically in Fig. 3 and a photograph of the experimental set-up is shown in Fig. 4. Strictly the method presented here measures operational deflected shapes rather than mode shapes. However, the structure is lightly damped and only loosely coupled to the excitation driver, so the vibration shapes are hypothesised to be close approximations of the shapes of the modes excited. The validity of this assumption will be explored via use of the Modal Assurance Criterion [20]. The procedure adopted was as follows: 1. The plate was prepared with a stochastic speckle pattern using black and white aerosol paint sprays following standard practice. The speckle pattern is shown in Fig. 5. 2. The plate was placed on the foam pads as described earlier with the loudspeaker below the centre of the plate but just out of contact. 3. The cameras were focussed and calibrated using the DIC system's internal procedures, with illumination provided by the DIC system's LED light source. The field of view (measured along the image centrelines) was 387 mm × 344 mm and the camera resolution 2452 × 2056 pixels, giving an average spacing between pixels of around 0.16 mm. 4. The loudspeaker was then driven at one of the plate's resonant frequencies with a sinusoidal signal from the frequency generator; maximum amplitude of the generator output was generally used to maximise the vibration amplitude obtained. The frequency was tuned manually to correspond to the natural frequencies of the plate, which were identified by ear and with the aid of a  small plastic bead which Bdanced^on the plate at resonance. 5. With the frequency set, a sequence of 60 images was then captured at a rate of 0.5 Hz to allow the flash unit to recharge between images. An electronic shutter speed (camera integration time) of 1/10,000 s was used. 6. Steps 3 and 4 were repeated until the data for the first six detected modes had been collected.
Although the approach shares with that of Fruehmann [9] the need to avoid the sampling rate being an exact fraction of the loading frequency, it was found in practice that the ratio of the two frequencies is so large (ranging from around 160 to around 2000) that an exact ratio is highly unlikely to occur and no special care was needed.

Processing of Results
The overall flow diagram for the process including the collection of experimental data is illustrated in Fig. 6. The algorithm for processing the results is described as follows. Firstly, standard DIC procedures are followed to convert the camera images into displacement fields, with the displacements being referenced (for convenience) against the first image in the sequence. A facet size of 17 × 17 pixels and a step size of 17 pixels were used, as shown in Fig. 5. The result of this is a set of files in HDF5 format, one relating to each image pair captured, each containing displacement maps for the x, y and z directions, along with the analogue data for the directly-measured and integrated quadrature signals.
The data from the directly-measured and integrated excitation signals are retrieved from the HDF5 files for each time point. For convenience, the integrated signal is re-inverted in software to correct the inversion inherent in the circuit in Fig. 1. The signal values are stored as a set of coordinate pairs, through which an ellipse is fitted (using the algorithm of Fitzgibbon et al. [17]) corresponding to the Lissajous figure relating these two sinusoidal quantities of equal frequency. This ellipse has its axes nearly, but not exactly, orthogonal to the Cartesian axes, and is characterised in terms of its offsets x 0 and y 0 , half major axis a, half minor axis b and its orientation θ to the x (directly-measured) axis ( Fig. 7(a)). An actual example of an ellipse from the experiments is shown in Fig. 7(b). The amplitudesv x andv y of the direct and integrated signals are: Various approaches exist for calculating the phase angle from a Lissajous figure (e.g. [16]); here, for generality, the phase angle ϕ of the integrated signal with respect For clarity this omits the procedure for rejection of outliers, which takes place within the sine fitting of facet displacements to the directly-measured one was calculated from the following: where: Additional decision-making is needed to ensure that the resulting phase angle is placed in the correct half-plane, with ϕ being replaced with π−ϕ if ϕ lies outside the range −π/2 ≤ ϕ ≤ π/2.
The value of ωt at each time point j is then calculated from the instantaneous values of v x and v x : This equation is solved using the two-argument inverse tangent function atan2 available in MATLAB and most other programming languages. The magnitude and phase of the i direction component of displacement vector u at each (k,l) cell position (i.e. pixel position on the displacement map) relative to the excitation signal is now calculated by firstly expressing the displacement vector at each time point j as a linear combination of sine, cosine and constant terms, with the (k,l) subscripts being implicit in the u, A, B and C terms: The least-squares estimates of A i , B i and C i are found from the set of n measurement readings for a given pixel position over the image sequence by solving the following overdetermined set of simultaneous equations using standard methods [21]: or for which the least-squares solution x i ≡[A i B i C i ] T is given by: The magnitude |u i | and phase angle ϕ i of each vector component i are then straightforwardly found: It is also straightforward within this analysis to calculate the coefficient of multiple regression in order to give a measure of the quality of fit of the estimated vibration data to the captured displacements: where J is an n×n matrix of ones. An alternative measure of the quality of results is the signal-to-noise ratio (SNR), which in the present publication is expressed as a simple ratio rather than in dB: Once the magnitude and phase of the vibrations are known for each pixel of the image, it is straightforward to create an animation of the deformed shape or a snapshot at maximum deflection. In the present case a sequence of 36 frames representing a complete cycle of displacements is stored in a new HDF5 file, along with maps of amplitude, phase, R 2 and SNR. In practice, it is found that the high frequency measurements in particular contain a significant level of noise. Grossly outlying points arising from failure of the DIC system to correlate a particular image correctly are straightforwardly eliminated if their displacement values greater than 200 mm, less than −200 mm or (with the exception of the first image in the sequence, which is taken as the datum) is equal to zero, noting that zero displacement values are used in the displacement maps to denote a point outside the region of interest. A separate, iterative procedure was used to in an attempt to eliminate less extreme outliers. This involved: 1. Initially calculating the least squares fit sine curve 2. Calculating the residual e j for each time point j and hence estimating the standard deviation σ from the mean squared error (noting that the denominator takes account of the constant term and two coefficients fitted): 3. Eliminating points lying further than a specified multiple of σ from the fitted curve 4. Re-fitting the sine curve to the remaining points 5. Repeating steps 2-4 a specified number of times.
In practice a value of five iterations each rejecting points lying further than 2 σ from the fitted curve was found to give a good compromise between excessive rejection of consistent data and retention of outliers.
For comparison purposes a finite element model has been constructed to allow the accuracy of the approximated mode shapes to be determined. The geometry was modelled in the CAD software CREO Parametric and the model analysed in CREO Simulate. The parameters were determined by measuring the plate's dimensions (350 mm × 225 mm × 2 mm) and by calculating the structure's density (2794 kg/m 3 ).

Results and Discussion
Raw DIC z (out-of-plane) displacement output for the five modes investigated are given in Fig. 8(a-f); frames representing the largest displacements have been chosen in order to show the kind of results which are obtainable with the use of the flash but without the use of the fitting software. At low frequencies the mode shapes are clearly visible but at 460 Hz and above the data are becoming swamped by noise and are of little direct use. Surprisingly, the DIC output at 1165 Hz appears to be of slightly better quality than those at 460 Hz and 781 Hz. The results of the analysis are presented in Fig. 8(g-l) and are compared with the corresponding mode shapes from the finite element analysis (scaled to be the same amplitude as the experimental displacements) in Fig. 8(m-r). Typical graphs of displacement vs ωt are shown in Fig. 9(a-d) both for positions at typical antinodes (where displacements are large and the SNR is also large) and as close as possible to nodal lines where displacement and SNR are both small. Maps showing the signal-to-noise ratio are shown in Fig. 9(e-f). It is noteworthy that values of SNR range from several hundred where displacements are large to values much less than 1 close to nodal positions. As a further comparison, Fig. 10 (a-d) shows the experimental deflected shapes for four frequencies plotted along the long centreline of the plate (corrected for perspective in the DIC results using an algorithm described by Chan [22]), and compared against the FE mode shapes scaled (using the same scale factors as in Fig. 8(a, c, d and f)) to match the experimental results.
It can be clearly seen that, for the first four modes at least, the system calculates deflected shapes which agree very well with finite element predictions of mode shape, with scatter typically being within around 5 % of predicted displacement at 83 Hz and 10 % at 358 Hz. A further comparison of the deflected shapes and the mode shapes is provided by applying the Modal Assurance Criterion (MAC) [20] to perspective-corrected maps of experimental deflected shape and interpolated maps of predicted modal shape as shown in Table 1. Taking the widely-used value of 0.9 as the threshold for good correlation between shapes, it is seen that the first four frequencies show good to excellent correlation between measured deflections and predicted mode shapes (MAC values between 0.92 and 0.99), the fifth is borderline (0.88) and the sixth is relatively poor (0.63). Visual examination both of the maps in Fig. 8(g-r) and of the graphs in Fig. 10(a-d) shows qualitatively that the deterioration in correlation is primarily due to local noise in the results rather than serious overall differences between the experimental deflected shapes and the predicted mode shapes, and it can be concluded that the operational deflected shapes in this experiment are generally very good approximations to the mode shapes. Some moirélike artefacts are visible in the raw DIC results and hence appear in the mode shape plots; these do not arise from the processing method described here but neither are  Fig. 8 (a-f) Raw displacement results from DIC system after correlation (g-l) Deflected shapes at maximum amplitude from present method (m-r) Finite element predictions of mode shapes. FE amplitudes are scaled to fit experiment they eliminated by it. This phenomenon is understood from the manufacturers 1 to be due to the presence of speckles smaller than 1-2 pixels in size, and can be reduced or eliminated by applying Gaussian smoothing to the raw images. This facility to achieve this is available in later releases of the ISTRA-4D software but not in the system used for these experiments. It is encouraging that the least squares fitting procedure enables the vibration displacement data to be recovered for frequencies up to around 1 kHz despite the raw displacement results being dominated by noise. The frequencies covered were at least as high as for previous reported work using a DIC system with high-speed cameras e.g. 70 Hz for reference [2] and 1060 Hz for reference [3]. It is noted that in general the low speed cameras have higher resolutions, and are thus able to measure smaller displacements resulting in less noise in the returned data. However, it is visible that the quality of measurement deteriorates as the highest attempted frequency of 1165 Hz is approached, and the signal-to-noise ratio for the 1165 Hz case is very poor, only exceeding 1 at local regions. Two reasons potentially account for the deteriorating quality.
1. The displacements become very small, typically less than 5 μm at antinodes going down to zero at nodes. Helfrick et al. [1] state that the sensitivity of a 3D DIC system for out of plane displacements is in the order of 0.03 pixels which here equates to around 5 μm, so the measurements are consequently at the limit of the system and close to the noise floor of the DIC. This limitation could be ameliorated by the use of larger excitation amplitude (at the risk of encountering nonlinearity in the system's behaviour) or higher resolution cameras. 2. The time period of the vibrations reduces to being approximately equal to the duration of the flash period (characterised in the specification as a half-life or t0.5 time of 1/750 s), potentially resulting in the displacement measurement being smeared over the flash period. In fact this issue was overcome in the present experiments by using a short exposure (integration) time in the cameras namely 1/10,000 s, but this was only possible because of the intense brightness of the flash. Higher specification flash units typically have shorter flash periods (e.g. 1/2900 s for the Bowens Gemini Pro 500 unit) so there is scope for investigating the use of such units for measurement of higher frequency vibrations, either to Bfreeze^the motion stroboscopically or to provide a more intense light pulse permitting even shorter exposure times.
When the present approach was proposed, it was not clear whether or not the use of a short-duration flash with a conventional scientific video camera would be effective, noting that this is not a normal use case for the camera. In fact there appears to be no problem in achieving this for an exposure time of 1/10,000 s and direct synchronisation of the flash and the camera trigger signal is all that is required. Although the present paper describes the use of flash units in conjunction with DIC as a means of measuring mode shapes, the technique is intended to be applicable more generally to measurements of stresses and strains in cyclic situations. Specifically, it could be used to extend the applicability of the dynamic DIC technique presented by Fruehmann [9].  A limitation of the present approach is that the structure's vibration modes must be found manually by searching for resonance. Although the details have not been implemented, there is scope for extending the method to enable it to find modal behaviour automatically. One approach would be to excite the structure with a recurring pattern of vibration containing a range of frequencies e.g. a swept sine, chirp or pink noise, and capture a large number of samples at time-stepped points within the repeating cycle. Effectively this would involve capturing the response sample by sample; conventional modal analysis methods would then be used to analyse the data. The shortcoming of this approach would be the very large number of images needed and the consequent storage issue. Further investigation is needed to determine whether such an approach, or something similar, could be implemented in practice.
One practical limitation of the present experiment was that the placement of the loudspeaker (close to the middle of the plate) resulted in only modes with two planes of symmetry being excited, so that resonances with a centrally-placed node, or with a nodal line on a plane of symmetry, remained undetected either audibly or using the Bdancing bead^and consequently no measurements were taken for those modes. This limitation is unrelated to the vibration measurement approach developed here, but it would be useful to avoid this issue in future experiments.

Conclusions
A technique for the measurement of small amplitude vibrations has been demonstrated for frequencies up to around 1 kHz. A method loosely based on the asynchronous method of lock-in thermography has been adopted. The novel features of the technique are the use of a flashgun synchronised with a free-running DIC system, and the use of an integrating amplifier in conjunction with a quadrature correction algorithm to provide the phase reference signal. Measured operational deflected shapes agree very closely (within a few percent) with predicted mode shapes at low frequencies though for higher frequencies there is considerable noise with values approaching 50 % of the vibration amplitude. A Modal Assurance Criterion (MAC) analysis shows that the measured deflected shapes are generally very good approximations of the mode shapes. Reasons for the frequency limit of the present system, and possible approaches for extending it, have been identified. The method is also intended to be applicable to measurement of strains in cyclic situations at higher frequencies than currentlyavailable techniques.