Microsoft Word - Article 8 - 62-521-1-Galley finale.doc ACTA IMEKO December 2013, Volume 2, Number 2, 41 – 47 www.imeko.org ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 41 A study of Savitzky-Golay filters for derivatives in primary shock calibration Hideaki Nozato1, Thomas Bruns2, Henrik Volkers2, Akihiro Oota1 1 National Metrology Institute of Japan, Umezono 1-1-1, 3058563 Tsukuba Ibaraki, Japan 2 Physikalisch-Technische Bundesanstalt, Bundesallee 100, 38116 Braunschweig, Germany Section: RESEARCH PAPER Keywords: shock; acceleration; infinite impulse response filter; Savitzky-Golay filter; derivative Citation: Hideaki Nozato, Thomas Bruns, Henrik Volkers, Akihiro Oota , A study of Savitzky-Golay filters for derivatives in primary shock calibration, Acta IMEKO, vol. 2, no. 2, article 3, December 2013, identifier: IMEKO-ACTA-02 (2013)-02-03 Editor: Paolo Carbone, University of Perugia Received February 6th, 2013; In final form June 30th, 2013; Published December 2013 Copyright: © 2013 IMEKO. This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited Funding: This work was supported by JSPS (Japan Society for the Promotion of Science) International Program for Young Researcher Overseas Visits. Corresponding author: Hideaki Nozato, e-mail: hideaki.nozato@aist.go.jp 1. INTRODUCTION The international standardization document ISO 16063-13 [1] describes shock and complex sensitivities of accelerometers determined by shock calibration. For primary shock calibration, it is important to derive an almost undistorted acceleration waveform using laser interferometry and digital signal processing. A setup involving a homodyne laser interferometer requires a digital filtering process associated with a derivative filtering to calculate acceleration from the measured displacement. To obtain a smooth acceleration waveform, ISO 16063-13 recommends the use of a 4th-order Butterworth low- pass (4th BW) filter followed by a difference method applied twice in sequence. The 4th BW filter works to remove high- frequency noise from the photo detector output signal and may also be used to suppress the resonant frequency of mechanical parts such as the anvil [2]. Conversely, a Savitzky–Golay (S–G) filter [3] can be applied to directly obtain a 2nd derivative with smoothing to produce a less distorted waveform of acceleration than the 4th BW filter. To evaluate their performance, numerical simulations were conducted for a homodyne setup using known analytical excitation functions. This paper reports the characteristics of acceleration measurements using the S–G filter. For the comparison, the validation of acceleration measurements using laser interferometry is also presented. 2. BASIC INFORMATION 2.1. Low amplitude shock Figure 1 illustrates the assumed acceleration waveforms over time and their respective frequency domain spectra for low- shock. The low-shock with a duration of 0.5 ms has acceleration components with several metre per second square between 5 kHz and 15 kHz. In this case, an appropriate cutoff frequency for the 4th BW filter would be 10 kHz. Thus, consideration of the spectrum of acceleration is important in shock calibration. ABSTRACT This manuscript reports the result of an investigation into two types of digital filters for numerical differentiation of the displacement signal in the case of primary shock calibration of accelerometers. The first is a difference method with a 4th-order Butterworth low- pass (4th BW) filter; the other is a Savitzky–Golay (S–G) filter, which applies a moving polynomial approximation. The computational comparison was applied to low and high amplitude shocks using known excitation functions. Each sum of residuals was compared for the optimized conditions of the 4th BW and S–G filters. The results of computer simulation indicated that the S–G filters exhibited better performance for the derivatives than the 4th BW filters. In addition, an analytical comparison using experimental vibration data also indicated that the S–G filters exhibited better derivative characteristics than the 4th BW filters. ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 42 The low-shock is generated by a collision motion between two rigid bodies: a hammer and an anvil through a rubber pad. An accelerometer to be calibrated is attached on an edge surface of the anvil, and a low-shock motion with a shape approximating a half-sine squared function is induced to the accelerometer by the collision. Typical peak acceleration ranges from 100 m/s2 to 10000 m/s2 with a duration of several ms [2, 4]. In the computer simulations, equation (1) was used to define the mathematical function for low-shock.          tAta low 2sin)( , ( < t < 2) (1) The value of A is the peak shock acceleration. The low-shock waveform is given by the square of the sine function, and the duration is defined as . In the range 0 < t <  and 2 < t < 5, the low-shock waveform vanishes. 2.2. High amplitude shock The high-shock is based on elastic wave propagation inside a long thin bar and is induced to an accelerometer attached on an edge surface of the bar where a reflection of the elastic wave occurs [5]. The shape of the high-shock waveform is normally described by the 1st derivative of the Gaussian velocity function [6]. In the computer simulations, the mathematical function for high-shock was defined by equation (2).                       2 2 2 10 2 5.0 10)(      t high e t eAta , (0 < t < 20) (2) Figure 2 presents the assumed acceleration waveforms over time and their respective frequency domain spectra for high- shock. Typical peak accelerations range from 10000 m/s2 to 100000 m/s2 with a duration of several dozen microsecond. The shock duration is defined as 4. The calculation range of acceleration is set at values from −10 to 10, such that the initial displacement becomes sufficiently small compared to the half-wavelength of the Ne–He laser. 3. NUMERICAL SIMULATION 3.1. Simulation procedure Figure 3 illustrates the basic procedure followed in the computer simulations. Quadrature signals with 400 mVp-p, white noise with 10 mVp-p, and 16-bit quantization for ±16 V are typical specifications for the homodyne laser interferometer in NMIJ [2]. The displacement waveform is calculated by computer simulation from the original acceleration waveform. The typical quadrature signals of the p and s waves are generated in accordance with this displacement waveform [7], including the effects of white noise and 16-bit quantization. Figure 4 presents a typical waveform of acceleration, displacement, and interferometric quadrature signals obtained for a low-shock with peak acceleration of 5000 m/s2 and duration of 0.5 ms. The identification stage starts by applying the Heydemann correction [8] followed by an arctangent demodulation to obtain the reconstructed displacement waveform. Decimation is applied for reduction of data amount. Then, the acceleration waveform is reconstructed through the differentiation process using the 4th BW and S–G filters. Thus, performance of the differentiation is investigated by computational comparison between the original and reconstructed accelerations. 3.2. 4th order Butterworth filter Figure 5 presents the calculated differences between the 4th BW (red line) and S–G (blue line) filters in the case of low- shock with a peak acceleration of 5000 m/s2 and a duration of 0.5 ms. The black line represents an assumed acceleration waveform. From this computational result, the assumed Figure 1. Computer simulations of low-shock waveforms and corresponding frequency spectra. Figure 3. Procedure of numerical simulation. Figure 2. Computer simulations of high-shock waveforms and corresponding frequency spectra. ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 43 waveform is generated with a sampling frequency of 50 MHz and the demodulated displacement is decimated to a sampling frequency of 10 MHz. The second and third panels in figure 5 are residuals and their spectrum compared with the assumed waveform, respectively. The cutoff frequency of the 4th BW filter is optimized to obtain the smallest residuals and is 10 kHz (see figure 6). Each minimal value depends on the duration. Shock waveforms have broad frequency bandwidth, as illustrated in figures 1 and 2; the longer the duration, the narrower is the frequency bandwidth covered by the shock waveform. The S–G filter is also optimized with 7 orders and 1400 side points (see figure 9). The second panel in figure 5 implies that the distortion of the S–G filter is smaller than that of the 4th BW filter. In the third panel, the S–G filter exhibits smaller residuals in both frequency ranges (below 3 kHz and beyond 10 kHz in the spectrum) compared to those of the 4th BW filter. Consequently, such a small deviation of waveform would be effective in obtaining more accurate results in terms of shock and complex sensitivities [9, 10]. Figure 7 illustrates the dependence of the deviation of the assumed peak acceleration from the computed peak acceleration for the 4th BW filter on the selected cutoff frequency. This result indicates that the peak acceleration that Figure 4. Example of acceleration, displacement, and interferometric quadrature signals (p and s waves) in computer simulation of low-shock with peak acceleration of 5000 m/s2 and duration of 0.5 ms. Figure 5. Assumed and two computed acceleration waveforms using the 4th BW and S–G filters for low-shock with duration of 0.5 ms. Deviation and spectrum between assumed and two optimized acceleration waveforms. 10-4 10-3 10-2 10-1 0 5000 10000 15000 20000 25000 30000 0.1 ms 0.5 ms 2.0 ms N or m al iz ed s um o f r es id ua ls (a rb .) Cut-off frequency (Hz) Figure 6. Dependence of the normalized sum of residuals on cutoff frequency for 3 different durations of low-shock. ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 44 can be measured by a laser interferometer is more accurate for longer durations. For a duration of 0.5 ms, a deviation of less than 0.1% corresponds to a cutoff frequency below 10 kHz. According to figures 6 and 7, a cutoff frequency of 10 kHz is suitable to obtain both a less distorted waveform and a more accurate peak acceleration, for low-shock durations of 0.5 ms. In order to more accurately measure the acceleration waveform using a homodyne laser interferometer, the authors recommend selecting a duration greater than 0.5 ms. However, it should be noted that the appropriate duration and cutoff frequency depend strongly on the mechanical resonance of each shock calibration machine [2]. Figure 8 shows the computed differences between 4th BW (red line) and S–G (blue line) filters in cases of high-shock with peak acceleration of 100 km/s2 and duration of 160 s. Additionally, this graph visually implies that a S–G filter with 7 orders and 440 side points can achieve smaller deviations of the acceleration waveform than can be achieved with the optimal 4th BW filter. 3.3. Savitzky-Golay filter Figure 9 illustrates the dependence of the sum of residuals on the order and number of side points of the S–G filter for low-shock, peak acceleration of 5000 m/s2, and duration of 0.5 ms. The sum of residuals using the 4th BW filter is minimized by optimizing the cutoff frequency (10 kHz) as shown in figure 6. The sum of residuals indicates that the S–G filter with 7 orders and 1400 side points is optimal under the chosen conditions. Table 1 summarizes the optimal value of each cutoff frequency in 4th BW filters and side point in S–G filters. -1.0 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 0 5000 10000 15000 20000 25000 30000 0.1 ms 0.5 ms 2.0 ms Pe ak a cc el er at io n de vi at io n (% ) Cut-off frequency (Hz) Figure 7. Dependence of peak acceleration deviation on the cutoff frequency for low-shock and 3 different durations. Table 1. Optimal values of cutoff frequency in 4th BW filter and window width in S–G filter. Peak [ms2] Duration [ms] Decimation [MHz] Cut-off [Hz] Side points Window [ms] 5000 0.1 50 21000 NA 5000 0.2 50 16000 4500 0.18 5000 0.5 50 10000 NA 5000 1.0 50 7000 NA 5000 2.0 50 5000 NA 5000 0.1 10 23000 600 0.12 5000 0.2 10 16000 900 0.18 5000 0.5 10 10000 1400 0.28 5000 1.0 10 7000 2000 0.40 5000 2.0 10 5000 2800 0.56 5000 0.1 1 22000 60 0.12 5000 0.2 1 17000 80 0.16 5000 0.5 1 10000 140 0.28 5000 1.0 1 6000 220 0.44 5000 2.0 1 4000 320 0.64 5000 0.1 0.5 24000 30 0.12 5000 0.2 0.5 15000 40 0.16 5000 0.5 0.5 9000 70 0.28 5000 1.0 0.5 6000 100 0.40 5000 2.0 0.5 4000 160 0.64 1000 0.1 10 14000 800 0.16 1000 0.2 10 10000 1300 0.26 1000 0.5 10 7000 2100 0.42 1000 1.0 10 5000 3000 0.60 1000 2.0 10 3000 4300 0.86 Figure 8. Assumed and two computed acceleration waveforms using the 4th BW and S–G filters for high-shock with duration of 160 s. Deviation and spectrum between assumed and two optimized acceleration waveforms. ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 45 “window”— the numerical value given by equation (3) —refers to window widths for polynomial approximation of a S–G filter. Accordingly, peak accelerations of 1000 m/s2 require lower optimal cutoff frequencies compared to those for 5000 m/s2. Therefore, the effect of high-frequency noise becomes more pronounced in cases of low peak acceleration.   1000/)1points2(  decimationsidewindow (3) The relationship between “cut-off” and “window” for low- shock, on the basis of Table 1, is plotted in figure 10. The closed circles represent peak acceleration of 5000 m/s2 for 4 different decimations from 0.5 MHz to 50 MHz. The open circles represent peak acceleration of 1000 m/s2. These results indicate that the relationship between optimal cutoff frequency and window width is independent of three types of variables: peak acceleration, duration, and decimation. Correspondingly, figure 11 illustrates this relationship for low and high shocks. The open circles represent high-shock, which exhibits different waveforms from low-shock. The peak accelerations are 20 km/s2 and 100 km/s2, and the duration ranges from 60 s to 300 s. Thus, the relationship indicates that a S–G filter works on low- and high-shocks purely as a low-pass filter (with differentiation). 3.4. Calculation time The S–G filter is a kind of finite impulse response (FIR) filter in which the calculation time depends on the side points. On the other hand, since the 4th BW filter is an infinite impulse response (IIR) filter with less filter coefficients, its calculation time is generally shorter than that of the S–G filter. So, two conditions in table 2 were selected to evaluate the calculation time between the difference method with the 4th BW filter and the S–G filter. Here, the software and operating system were LabVIEW 2009 and Windows 7 with the central processing unit of an Intel core i7-2620M at 2.7 GHz. Then, the calculation time using the difference method with the 4th BW filter and S–G filter was 1.029 ms and 123.087 ms in the decimation of 10 MHz. Also, that was 0.171 ms and 3.429 ms in the decimation of 0.5 MHz. Information of the calculation time was obtained using the profile performance in tools. From the results, the S–G filter took longer time than the 4th BW filter, and each calculation time drastically decreased by the data reduction due to the decimation. 4. INVESTIGATION USING EXPERIMENTAL DATA Experimental data in shock and vibration calibrations were examined to indicate the applicability of S–G filters, even when the original waveform is unclear. Fundamentally, it is difficult to discern the correct amplitude and frequency from the shock 10-4 10-3 10-2 10-1 0 1000 2000 3000 4000 5000 5th order 7th order 9th order 4th BW N or m al iz ed s um o f r es id ua ls (a rb .) Side points Figure 9. Sum of residuals between optimized 4th BW filter and S–G filters of 3 different orders for low-shock and duration of 0.5 ms. Table 2. Conditions for comparison of calculation time between difference with 4th BW filter and S–G filter. Peak [ms2] Duration [ms] Decimation [MHz] Cut-off [Hz] Side points Window [ms] 5000 0.5 10 10000 1400 0.28 5000 0.5 0.5 9000 70 0.28 10-2 10-1 100 103 104 105 5000 m/s2 1000 m/s2 W in do w w id th ( m s) Cut-off frequency (Hz) Figure 10. Relationship between optimal cutoff frequencies and window widths in low-shock. 10-2 10-1 100 103 104 105 low shock high shockW in do w w id th (m s) Cut-off frequency (Hz) Figure 11. Relationship between optimal cutoff frequencies and window widths in low- and high-shocks. Figure 12. Typical experimental waveforms for shock and vibration calibrations. ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 46 waveforms, as shown in figure 12. Conversely, a fundamental frequency is known in case of vibration waveforms. Therefore, in order to obtain a reference to the fundamental frequency from the vibration waveforms, a fast Fourier transform was applied to the experimental data of the NMIJ high-frequency calibration system [11]. Figure 13(a) presents an example of the displacement (black line) and its acceleration waveforms that are derived experimentally using the difference method with 4th BW (red line) and S–G (blue line) filters at the fundamental frequency of 5 kHz. Here, the two acceleration waveforms were derived from the displacement waveform. The cutoff frequency of the 4th BW filter was 10 kHz and the side point of the S–G filter with a 7th-order polynomial was given according to the relationship illustrated in figure 11. Thus, the displacement waveform becomes the standard against which to evaluate whether a S–G filter can differentiate more accurately than a 4th BW filter. Figure 13(b) shows a deviation between nominal and differentiated acceleration waveforms using 4th BW and S–G filters. Firstly, the nominal acceleration waveform )(tD is derived by applying equation (4) to the displacement waveform. ftbftbbtD  2cos2sin)( 210  (4) where 0b , 1b and 2b are parameters for sine-approximation, and f is the fundamental frequency. Next, the nominal acceleration waveform is calculated by the second derivative of equation (4). As figure 13(b) shows, the S–G filter seems to have smaller deviation for the derivative than the 4th BW filter. Figure 14 represents the spectrum of the displacement and accelerations calculated on the basis of each waveform in figure 13(a). Figure 15 illustrates the results of comparison between the 2nd derivative of the 4th BW filter and that of the S–G filter on the basis of this spectrum. The reference value Dis (black circle) of the acceleration amplitude was calculated from the displacement waveform according to the following equation:   dfDis 22  (5) The parameter d is the displacement amplitude. The S–G filter is better than the 4th BW filter in two respects: the first is that the acceleration amplitude due to the S–G filter is closer to the reference value than that of the 4th BW filter over the entire range of cutoff frequencies (i.e., up to 50 kHz); the other is that the total harmonic distortion of the S–G filter is smaller than that of the 4th BW filter. On this basis, it can be concluded that the 2nd derivative of the S–G filter performed better than that of the 4th BW filter in obtaining the acceleration waveform for shock calibration. Figure 15. Comparison of acceleration amplitude and total harmonic distortion between the 4th BW and S–G filters. Figure 13. (a) Displacement and acceleration waveforms using 4th BW and S–G filters. (b) Deviation between nominal and differentiated acceleration waveforms using the 4th BW and S-G filters. 10-14 10-13 10-12 10-11 10-10 10-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 101 102 103 103 104 105 Displacement 4th BW S-G D is pl ac em en t ( m ) A cc el er at io n (m /s 2 ) Frequency (Hz) 5 kHz Figure 14. Spectrum of displacement and accelerations calculated using the 4th BW and S–G filters. ACTA IMEKO | www.imeko.org December 2013 | Volume 2 | Number 2 | 47 5. SUMMARY An approach for the calculation of derivatives using a S–G filter to obtain appropriate waveforms of acceleration in low and high shocks was proposed and discussed on the basis of computer simulations. The obtained waveform was compared with waveforms that were differentiated and smoothed twice using a 4th BW filter according to the procedure proposed in [1]. A relationship was derived from each optimal residual summation between 4th BW and S–G filters with 7th orders; this relationship did not depend on the decimation, peak acceleration, duration, or waveform in cases of low or high shocks. The results imply that the S–G filter produces better derivative and low-pass characteristics than the 4th BW filter for both low and high shocks. An analytical comparison between the derivatives of the 4th BW and S–G filters was investigated using experimental vibration data; fast Fourier transform was applied to the displacement and two acceleration waveforms of the 4th BW and S–G filters. This analytical comparison also indicated that the acceleration amplitudes of the S–G filter were closer to the reference values derived from the displacement than those of the 4th BW filter. Moreover, the S–G filter was considered superior in terms of total harmonic distortion. These results could lead to a more accurate measurement of shock and complex sensitivity for primary calibration of accelerometers implemented according to [1]. ACKNOWLEDGEMENTS The authors would like to thank the staff of the working groups “Realization of acceleration” and “Impact Dynamics” of PTB for their helpful partnerships and constructive suggestions. This paper was drawn up on the basis of an international collaboration between PTB and NMIJ. This work was supported by the JSPS (Japan Society for the Promotion of Science) International Program for Young Researcher Overseas Visits. REFERENCES [1] ISO 16063-13:2001 Methods for the calibration of vibration and shock transducers – Part 13: Primary shock calibration using laser interferometry. [2] H. Nozato, T. Usuda, A. Oota, T. Ishigami, Calibration of vibration pick-ups with laser interferometry: part IV. Development of a shock acceleration exciter and calibration system, Measurement Science and Technology, vol. 21, no. 6 (2010), 065107. [3] A. Savitzky, M.J.E. Golay, Smoothing and differentiation of data by simplified least squares procedures, Analytical Chemistry, vol. 36, no. 8 (1964) pp. 1627-1639. [4] Y. Huang, J. Chen, H. Ho, C. Tu, C. Chen, The set up of primary calibration system for shock acceleration in NML, Measurement, vol. 45, no. 10 (2012) pp. 2383-2387. [5] R.M. Devies, A critical study of the Hopkinson pressure bar, Philosophical Transactions of the Royal Society A, vol. 240, no. 821 (1948) pp. 375-457. [6] K. Ueda, A. Umeda, H. Imai, Uncertainty evaluation of a primary shock calibration method for accelerometers, Metrologia, vol. 37 (2000) pp. 187-197. [7] H. Nozato, T. Usuda, A. Oota, T. Ishigami, The methods for the calibration of vibration pick-ups by laser interferometry: part V. Uncertainty evaluation on the ratio of transducer’s peak output value to peak input acceleration in shock calibration, Measurement Science and Technology, vol. 22, no. 12 (2011), 125109. [8] P.L.M. Heydemann, Determination and correction of quadrature fringe measurement errors in interferometers, Applied Optics, vol. 20, no. 19 (1981), pp.3382-3384. [9] T. Bruns, A. Link, F. Schmähling, H. Nicklich, C. Elster, Calibration of acceleration using parameter identification – Targeting a versatile new standard, proc. of XIX IMEKO World Congress, Lisbon, Portugal (2009), pp. 2485-2489. [10] A. Link, A. Täubner, W. Wabinski, T. Bruns, C. Elster, Calibration of accelerometers: determination of amplitude and phase response upon shock excitation, Measurement Science and Technology, vol. 17, no. 7 (2006), pp. 1888-1894. [11] A. Oota, T. Usuda, H. Nozato T. Ishigami, Development of primary calibration system for high frequency range up to 10 kHz, proc. of IMEKO 20th TC3, 3rd TC16 and 1st TC22 International Conference, Merida, Mexico (2007), ID-103.