- Email: [email protected]
Phase spectrogram and frequency spectrogram as new diagnostic tools
ARTICLE IN PRESS Mechanical Systems and Signal Processing Mechanical Systems and Signal Processing 21 (2007) 125–137 www.elsevier.com/locate/jnlabr/ymssp
Phase spectrogram and frequency spectrogram as new diagnostic tools Franc- ois Le´onard Institut de recherche d’Hydro-Que´bec, 1800 boul. Lionel-Boulet, Varennes, Que´bec, Canada J3X 1S1 Received 11 April 2005; received in revised form 2 August 2005; accepted 5 August 2005 Available online 23 September 2005
Abstract The power spectrogram has found wide use in diagnostics based on mechanical vibrations. It makes use of the amplitude values of the Short Time Fourier Transform (STFT) but does not consider the information contained in the phase values. In the phase spectrogram, on the other hand, it is precisely the phase values of the STFT that are processed so that they can be presented graphically and are useful. The time–frequency distribution of the phase of a specific component supplies information about the phase modulations around a reference point, which determines both the reference phase and the reference frequency for this component. The frequency spectrogram is calculated from the phase difference between each time slice of the STFT. The frequency spectrogram shows the drift on the instantaneous frequency of each spectral component. These tools have proven to be useful complements to the power spectrogram. In crack detection, for example, they accentuate any frequency drift that occurs in the damped free-decay response of a cantilever beam, which appears to be a function of the crack depth. When identifying a sporadically pulsed vibration source, the ability to show the phase coherence from one impact to the other in time on the spectrogram makes it possible to draw a conclusion about the presence of a stable-vibration source hidden in the background. A dual-channel phase spectrogram is also presented with two application examples. Published by Elsevier Ltd. Keywords: Spectrogram; Phase spectrogram; Frequency spectrogram; STFT; Phase unwrapping; Modulation; Instantaneous frequency; Dual-channel time-frequency distribution; Crack detection
1. Introduction The time–frequency representation was first mentioned in the work of Gabor [1] and Ville [2] in the mid1940s. Gabor represented the signal in the time–frequency plane as a distribution of ‘‘quanta of information’’ whereas Ville was in search of a means of obtaining an energy density. Since then, many efforts have been invested to combine and consolidate the different approaches under Cohen’s class [3]. However, although certain algorithms make use of the phase in the time–frequency plane [4–6], no author has come up with a way to graphically represent the phase in this plane.
E-mail address: [email protected]. 0888-3270/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.ymssp.2005.08.011
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
126
Nomenclature i index of the spectral line m,k index of the time slice mr index of the reference time slice, fixed value fs sampling frequency, fixed value f m;i ðPÞ instantaneous frequency of the wrapped frequency spectrogram fU m;i ðP; mr Þ instantaneous frequency of the unwrapped frequency spectrogram l number of time samples between two consecutive time slices, fixed value P, Pk integer gain, constant value Fm;i ðP; mr Þ instantaneous unwrapped phase related to the reference time slice mr W; Wm;i integer number of phase rotations in the trigonometric circle Sm;i STFT result, amplitude complex value y; ym;i unwrapped phase y; ym;i wrapped phase yG m;i ðP; mr Þ instantaneous phase of the generalised unwrapped phase spectrogram
To start with, we need a time–frequency distribution of the amplitude and phase of the signal. A fair number of algorithms are able to provide this distribution but, for reasons of popularity and simplicity, we take as our starting point the Short Time Fourier Transform (STFT). Usually the term ‘‘spectrogram’’ is reserved for presentation based on a power distribution of the STFT which, on a time signal s(t), is written: Z 1 F s ðt; nÞ ¼ sðtÞ h ðt tÞej2pnt dt, (1) 1
with F s ðt; nÞ denoting the complex distribution function and h ðt tÞ, the conjugate of the spectral window used as a time–frequency kernel. The resulting product F s ðt; nÞ F s ðt; nÞ gives the spectrogram. In this text, we will refer to the ‘‘power spectrogram’’ when we refer to this product, as opposed to the phase spectrogram and frequency spectrogram, which will be introduced later. The STFT is a reversible transform since we need only a surface integration in the time–frequency plane [6] such as Z Z 1 sðtÞ ¼ F s ðt; nÞ hðt tÞej2pnt dt dn (2) 1
in order to find the time signal s(t). On the other hand, it is much more difficult to reconstitute the signal exactly using just the power spectrogram [7] because the latter does not contain all the information needed to reconstruct the original time signal: what is missing is the phase information contained in the ‘‘phase’’ argument of the STFT. It should be noted that interference between the components present can be observed on the power spectrogram. Actually, it is the phase information that this interference contains that allows us to estimate the time signal from the power spectrogram. Although the latter does contain a certain amount of phase information, it is nevertheless incapable of providing an adequate illustration of this phase. This paper explains how the phase spectrogram, using the argument of the STFT, provides access to a source of information which completes that offered by the power spectrogram. The information in a phase spectrogram can take different forms. The easiest to understand intuitively is the phase variation of a sine wave compared to a reference sine wave, the latter being the same sine wave taken at a given time and not having changed its frequency since. On the other hand, in deriving this phase, we generate a frequency spectrogram that gives us a different perception of the instantaneous frequency.
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
127
2. Phase unwrapping Before we tackle the algorithms of the phase and frequency spectrogram, we should clarify the mathematics involved in phase unwrapping. When we extract a phase value from a complex amplitude in a Fourier domain, we obtain its projection y in the trigonometric circle. The unwrapped phase value y ¼ y þW 2p
with W 2 Z
(3)
differs by W rotations in the trigonometric circle; we know the position of the phase on the circle but not the number of rotations. Phase unwrapping amounts to finding this number, which is done progressively by tracing round the trigonometric circle repetitively, usually starting from a chosen reference point. The unwrapped value y only has any meaning if we can associate it with a reference that has physical significance. For example, in cepstral and modal analysis, phase unwrapping is done on the frequency axis with y ¼ 0 at 0 Hz as the reference point. This is unavoidable because the d.c. component has a null imaginary value and, also, because a time delay applied to the 0 Hz frequency does not result in a phase shift. On the other hand, the choice of reference is more arbitrary for the phase spectrogram where phase unwrapping follows the time axis in the time-frequency plane. The phase here never seems smooth without appropriate unwrapping [8]. Several publications deal with phase unwrapping on the frequency axis [7,9–12] but no author presents the phase unwrapping on the time axis in the time–frequency plane. However, whether it be the frequency or the time axis, phase unwrapping consists in selecting the shortest path on the trigonometric circle. Usually, the unwrapping can only be achieved if the phase does not vary by more than p radians between two successive time slices. The problem is that two successive time slices of a STFT can be several time samples apart. Let l be the step size in number of samples separating two adjacent time slices of the STFT. From Nyquist criterion, the highest frequency that sampling without spectral aliasing can withstand is equal to a half sampling frequency, which produces an increase of p radians between each time sample for this component. The latter will therefore have a phase increase of lp radians between two successive time slices of the STFT. It seems, therefore, that the phase increase exceeds p for frequencies higher than f s =ð2lÞ, where f s is the sampling frequency. At first sight, phase unwrapping appears difficult if not impossible on the time axis of a STFT but, as we will see later, the difficulty can be overcome. In this text, we use underscored notations to mean that the y phase y is contained within the trigonometric circle. The underscored operator n ¼ mod2p fn þ pg p
(4)
therefore removes the 2p multiple value in the operand. We should also make it clear that the phase ym;i comes from the trigonometric function atan2 (*,*) applied to the result of the STFT such that ym;i ¼ atan 2ðReðSm;i Þ; ImðS m;i ÞÞ
(5)
with ym;i 2 p; p, where m is the index of the time slice and i that of the spectral line. In calculating the STFT, application of a spectral window comprising N samples will give N/2 spectral lines such that i 2 f0; . . . ; N=2 1g is the spectral-line index. The phase-unwrapping algorithm in time is recursive and usually starts with the first time slice. A simple form of this algorithm, starting with the first sample, consists in writing ynm;i ¼ ym;i þWm;i 2p
with Wm;i 2 Z,
(6)
where Wm corresponds to the number of rotations or turns completed in the trigonometric circle. The integer Wm;i is established recursively such that 8 W for po ym;i ym1;i pp > < m1;i 1 þ W for ym;i ym1;i p p m1;i and W1;i ¼ 0, Wm;i ¼ (7) > : Wm1;i -1 forpo ym;i ym1;i so that any phase shift with an amplitude exceeding p can be attributed to the fact that the boundary p=p of the trigonometric circle has been crossed. If the phase sampling condition is respected, then ym;i ¼ ym;i . For
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
128
simplification, we write ym;i ¼ unwrap ym;i
(8)
to express the unwrapping of the phase distribution ym;i . The operator unwraps the phase without worrying about there being possibly n extra 2p radians between two consecutive samples. Lastly, in order to refer the unwrapped phase Wm;i to the reference time slice mr , all we need to do is to subtract the latter phase such that ym;i ymr ;i gives the unwrapped phase in relation to the reference slice. 3. Frequency spectrogram The wrapped frequency spectrogram f f m;i ðPÞ s P ym;i ym1;i 2plP
(9)
is obtained from the phase difference between successive time slices of the STFT where P 2 Z. A mathematical demonstration has been developed and can be found in [13]. The gain P multiplies the sensitivity expressed as the fringes per Hz in such a way that fs fs ; f m;i ðPÞ 2 (10) 2lP 2lP defines the frequency axis scaling. All we need to do is display the frequency on a scale of colours where the two extremes f s =2lP and f s =2lP have the same colour, in order to watch the frequency develop without any break in the progression. The gain P is adjusted so that the progression can be monitored and we can easily count the fringes on the screen. The sensitivity f s =lP is expressed in Hz/fringe. that this definition of the frequency spectrogram depends largely on the equality We should emphasise
P ym;i ym1;i ¼ P ym;i ym1;i þW2p
for integers W and P. With an integer value gain, we amplify the
phase difference between two time slices without any disturbance from an error in the estimation of a phase of ‘‘W’’ turns in the trigonometric circle: we no longer need to know the number of rotations of the unwrapped phase value before applying the gain. Once an integer value gain has been applied, we can unwrap the phase. In fact, the unwrapped frequency spectrogram fs fU unwrap P ym;i ym1;i unwrap P ymr ;i ymr 1;i m;i ðP; mr Þ 2plP for m ¼ 2; 3; 4; . . . ; M ð11Þ allows us to display the frequency relative to the reference slice mr on a continuous scale. The frequency presents the same topological characteristic as the phase in the trigonometric circle, which means that it can be illustrated as superposed on a cyclic scale (Eq. (9)) or unwrapped on a conventional scale (Eq. (11)). The term P ym;i ym1;i in Eqs. (9) and (11) can be replaced by NP X
Pk ymþk;i ymþk1;i
(12)
k¼1
for weighting the computation of the phase derivative. This smoothing distributes the gain P according to an approximately Gaussian distribution of smaller weight PI so that P¼
NP X k¼1
Pk
(13)
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
129
Fig. 1. Power spectrogram and frequency spectrogram of the frequency-modulated signal cosð2p 200 Hz t þ 50 expð2tÞ sinð2p 4 Hz tÞÞ sampled at 2000 samples/s (N ¼ 256 for a 4-term Blackman–Harris window, spectral interpolation 2 , l ¼ 3, P ¼ 33).
with integers Pi. The advantage of the weighting is that we obtain smoothing with a Gaussian distribution whereas increasing the distance l between slices amounts to having a smoothing with a rectangular window where P ¼ Pi N p . The phase values from the STFT fill a matrix measuring M I, where M is the total number of time slices, I ¼ N=2 corresponds to the number of spectral lines. The term ymþk;i ymþk1;i then becomes the simple subtraction of a matrix from the same matrix offset by a time slice, namely M I operations. Applying the operator n to Pk ymþk;i ymþk1;i , like the operator unwrap, requires very little computation time, namely a
multiple of MI operations. On the whole, the computation effort increases in proportion to the dimensions of the time–frequency matrix. Fig. 1 depicts a power spectrogram with, below it, the unwrapped frequency spectrogram of a numerically generated cosine, which is frequency-modulated by a damped sine wave. Since the instantaneous frequency of a sine wave has no significance other than below the main lobe of the corresponding spectral component, only the result under the main lobe is displayed. In order to do this, we fill in the corresponding power spectrogram valleys in black to drown the secondary components with the lowest amplitudes. The height of the filling is adjusted manually by the user, depending on the amount of information to be underlined. The LabViewTM code and runtime is available on request. In computing the STFT, we recommend applying a spectral window having high side lobe rejection. A typical lobe comprises 5–7 spectral lines, which gives a frequency plateau of the same width if the lobe emerges sufficiently from the filling level. Obviously we can enhance the graphic resolution of the plateau by increasing the number of spectral lines beneath it. We need to only perform a spectral interpolation, for example by increasing the width of the spectral window with zero padding. We should clarify that the number of samples N here is the number before zero padding. For example, an interpolation of 2 with N ¼ 256 means that 256 zeros are added and that the TRF will apply to 512 samples, which gives 256 spectral lines. The modulation of phase or frequency generates a special fringe pattern: the fringe appears perpendicular to the path of the sine wave in the time–frequency plane, as illustrated in Fig. 1. Fig. 2, by contrast, shows that a
ARTICLE IN PRESS 130
F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
Fig. 2. Frequency spectrogram of a numerically generated sine wave, modulated in amplitude and sampled at 2 ks/s. (2 Hz/fringe, N ¼ 256 for a 4-term Blackman–Harris window, spectral interpolation 4 , l ¼ 3, P ¼ 333).
Fig. 3. Two different analyses of the same part of signal illustrated in Fig. 1 showing two apparent angles of the fringe alignment with the path of the sine wave.
modulation in amplitude has a different pattern: fringe waves in parallel to the path of the sine wave. Here a 4 spectral interpolation is required to be able to display this fringe across the entire width of the plateau. Caution is necessary, however, when we come to interpret the fringe patterns. As shown in Fig. 3, for a same frequency resolution the perspective given by two parameters sets changes the apparent angle between the fringe pattern and the path followed by the sine wave. Lastly, a recent paper [13] has shown that the variance on the frequency is minimal at the centre of the frequency plateau and is of the same order of magnitude as the variance obtained using the spectral interpolation algorithms (IFFT), namely near the Crame´r–Rao bound on the frequency estimate.
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
131
Fig. 4. Phase spectrogram of a sine wave with a phase modulation of 7p radians, i.e. the signal cosð2p 200 Hz t þ p cosð2p 2Hz tÞÞ with f s ¼ 1 s/s (N ¼ 128 samples, 2 spectral interpolation, l ¼ 5 and P ¼ 1).
4. Phase spectrogram The phase difference is equal to 2p
lP U f ðP; mr Þ þ C f s m;i
(14)
between two successive time slices where C is an unknown corresponding to the number of complete turns in the trigonometric circle between two successive slices of the STFT. The mean frequency of the component thus determines the magnitude of C. The information that C contains is of little use because the focus of our interest is the phase variations around the mean increment in the phase increase, not the slope of the phase increase due to the mean frequency of the component. Therefore, if we remove C from the computation, numerical integration of the type 8 0; m¼1 > < m X Fm;i ðP; mr Þ ¼ 2p l (15) fU m ¼ 2; 3; . . . ; M m;j ðP; mr Þ; > : fs j¼2 yields an unwrapped-phase value starting from the first time slice. In Eq. (15), the phase calculation is referenced to both the first slice for calculating the phase and slice mr for calculating the reference frequency. The latter frequency determines the Dy to be removed between each slice so as to eliminate the constant term of the phase increase. The generalised phase spectrogram defined as yG m;i ðP; mr Þ Fm;i ðP; mr Þ Fmr ;i ðP; mr Þ
(16)
allows us to refer both the phase calculation and the frequency calculation simultaneously to the mr th time slice. For a display with a succession of fringes, we use this presentation yG m;i ðP; mr Þ but we could also write yG m;i ¼ Fm;i ðP; mr Þ
m1 FM1;i ðP; mr Þ M 1
(17)
if we want to impose a null phase value on the two extremities of the measurement. If the distance l and the distribution Pi have an influence on the variance in the frequency spectrogram, they have none on the phase spectrogram [13]. In the latter case, the distance l affects the graphic definition on the time axis while the weight P restrains the phase to the interval ]p=P; p=P], superposing the distribution of phase values for no useful purpose. The weight P is therefore set at unity for calculating the phase spectrogram and the feature proposed at Eq. (12) is not used. The numerical example in Fig. 4 shows the presence of a uniform phase plateau across the width of the main lobe of the component. Also, as in the frequency spectrogram with a phase modulation, the phase spectrogram
ARTICLE IN PRESS 132
F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
Fig. 5. Phase spectrograms of the free response of beams excited by a mechanical impact. The vertical line in each figure shows the position of reference mr (2p radians/fringe, with f s ¼ 25:6 ks=s, N ¼ 256 samples, 8 spectral interpolation, l ¼ 5 and P ¼ 1).
of a phase modulation shows perpendicular fringes in the path of the sine wave while an amplitude modulation shows fringes in the opposite direction on the phase spectrogram (not illustrated. See [13]). Lastly, the use of coloured fringes makes it far easier to read the graphics. With monochrome displays we had to use a different shade of grey for the extremities of the display range to determine the direction of the phase increase. With colours we can have colour continuity when we cross the boundary from p to p. All we need to do is to use three or more colours in the range concerned in order to determine the direction of the phase or frequency progression by counting the fringes.
5. Applications to monitoring and diagnostics The phase and frequency spectrograms were developed in the framework of a crack detection project to reveal slight changes in the frequency of the free damping response of excited modes [14,15]. These frequency changes result in phase shifts, as illustrated in Fig. 5. The modes show a constant phase in the free response of the uncracked beam while a phase drift appears with the other three cracked beams and increases with the surface area of the crack, occurring mainly when the vibration amplitude is no longer enough to completely open the crack. To open a fatigue crack in steel we need to exert a force or it will remain tightly closed. Three vibration states succeed each other on the damped free response of a steel beam when the amplitude of the impact is high enough to open the crack straightaway. Fig. 6 shows where these states are located on a frequency spectrogram of the damped free response of the steel beam that is most badly cracked. In the first state, the crack opens completely and closes alternately with the vibration beat. Since there is less rigidity when the crack is open, the mean frequency of the first vibration is lower. In the second state, the vibration amplitude allows only partial opening of the crack; it is in this second state that the damping reaches a maximum and the frequency varies most. Lastly there comes a point where the vibration amplitude is clearly not capable of
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
133
Fig. 6. Frequency spectrogram of the free response of a badly cracked beam excited by a mechanical impact (f s ¼ 25:6 ks=s).
Fig. 7. Frequency spectrogram of an acoustic signal measured remotely by a microphone (above) and at the emission source by an accelerometer (below).
opening the crack and the third state begins, characterised by the highest modal frequencies and less damping than in the first two states. In the previous example, the amplitude of the transient signal is multiplied by an exponential gain increasing with time in order to keep the amplitude of the vibration components near a constant value. This not only gives a more uniform plateau width but also reduces the effect of free damping on the amplitude modulation. The sound wave propagating in a moving medium, unlike light, does not travel at a constant speed in all directions for all observers. The wind, for example, can change the speed of the propagation of a sound wave and thereby alter the time delay between emission and reception, so here we have phase modulation. To this we have to add a slight amplitude modulation, which is difficult to observe on a power spectrogram. Here again, the frequency spectrogram appears as a useful tool for displaying and quantifying the physical phenomenon as shown in Fig. 7, which presents the vibration of one of the dominant modes of one of the sheds of a high-voltage porcelain insulator in a distribution substation. These modes had been excited by the impact of a small metal ball. The wave distortion here corresponds to wind speed fluctuations of around 1 m/s between the source and the measurement. This phenomenon cannot be seen in the lab measurements but it made us take another look at the algorithms developed in the laboratory for detecting cracked insulators. In the area of vibration diagnostics, we are sometimes called upon to determine whether the intermittent bursts observed on a spectrogram come from a mode excited randomly or from some distant excitation with a stable frequency but transmitted sporadically. Fig. 8 shows the phase coherence from one burst to another in
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
134
Fig. 8. Power spectrogram and phase spectrogram (2p radians/fringe) of an acceleration signal recorded in the cabin of an airline pilot. This measurement was kindly provided by Rene´ Archambault IMS. (f s ¼ 24 Ks=s, N ¼ 128, l ¼ 3, P ¼ 1).
time, revealing the presence of a stable-vibration source that is observed sporadically. Vibration coming from a resonance excited by random shocks would not have this phase continuity. 6. Dual-channel applications Dual-channel applications [16] offer the opportunity to observe the frequency gap or the phase gap between different channels or measurements in the time–frequency plane. The required condition is that the signals have similar fingerprints in the time–frequency plane so that if we want to compare sine waves present on two channels or measurements, these must follow the same trajectories in the time– frequency plane. If need be, we can have a time translation in order to obtain the superposition of the fingerprints. We cannot resort to subtraction of the phase argument of two STFTs to observe the dual-channel phase because this phase difference does not fulfill the phase sampling condition. For example, if we cross the boundary p=p on a channel at a different time than on another channel, a phase jump can occur which can be interpreted in two ways: either as one of the channels crossing the boundary or as a real difference by more than a value p between the two channels without any crossing of this boundary in the trigonometric circle. Unwrapping of the phase on each channel must be done before phase subtraction or else information will be lost, resulting in ambiguous interpretation of the phase increases. Subtraction of two unwrapped phase spectrograms allows us to present the evolution of the phase between sinusoids from different measurements that are superposed in the time–frequency plane. To be more precise, the dual-channel phase spectrogram here gives the pattern of the relative difference in the phase evolution of two juxtaposed sinusoids and not the phase difference itself. If we need to show the phase difference between two sinusoids, the process must include the phase of the reference spectral slice mr common to each channel (A and B) by adding the phase vector P ymr ;i prior to subtracting the phase spectrograms so that yAB dualchannel
1 A y yB , P
(18)
where yA ¼ B
y ¼
P yA mr ;i þ P yBmr ;i
þ
m X
! mr X A A unwrap P yA unwrap P yA , m;i ym1;i m;i ym1;i
j¼2
j¼2
m X
! mr X B B B B unwrap P ym;i ym1;i unwrap P ym;i ym1;i .
j¼2
j¼2
Thus, we obtain a phase difference between components, which is more practical for such purposes as modal analysis of an unsteady frequency response, for example.
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
135
As for the dual-channel frequency spectrogram, this could be generated by the difference A B f AB dualchannel f m;i ðPÞ f m;i ðPÞ
(19)
between the phase derivatives calculated on two STFTs each of which corresponds to a measurement channel. This wrapped instantaneous frequency gap appears significant where the components are superposed. By contrast, the unwrapped frequency gap Unwrapped Unwrapped f AB fA ðP; mr Þ f Bm;i Unwrapped ðP; mr Þ m;i dualchannel
(20)
gives a significant result for the frequency lines included in the frequency plateau of the component at the intersection of the time slice mr .
Fig. 9. Time trace, single- and dual-channel frequency spectrograms and dual-channel phase spectrogram of a frequency sweep test with an imposed displacement on a HV cable.
ARTICLE IN PRESS 136
F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
Fig. 10. Dual-channel frequency spectrogram comparing the phase evolution between the responses to two impacts on a porcelain insulator observed by the same channel.
The concept of the dual channel for phase or frequency spectrograms is wider than the same concept for a dual-channel frequency analyser. In fact, in the latter case, we are comparing the phase frequency distribution of two channels sampled simultaneously, for one time slice. With the spectrogram, however, we have an extra dimension: time. The dual-channel frequency or phase spectrogram allows us to observe the temporal evolution of a frequency or phase relation of time–frequency components. All that is required is that the components compared share the same spectral position and some coherence between them (in other words, a linear physical relationship). The components compared may occupy different temporal positions. For example, the time–frequency components detected by two remote microphones must be juxtaposed in time: note that any positioning error here causes interference, which reduces the surface area where there is coherence between the components. The components may even come from different tests. For example, we can compare two transients present on the same channel to observe whether they have the same frequency or even phase modulation characteristics. While today’s modal analysis tools are not adequate for us to tackle unsteady frequency modes, the dualchannel phase spectrogram makes it possible to show the phase relation between channels without accounting for frequency fluctuations, as long as the latter do not move by more than three spectral lines. In addition, the dual-channel phase spectrogram also allows us to observe the phase behaviour of a mode between the fixed and mobile parts on machines. We can compare two time–frequency components present on channels sampled in parallel using two power spectrograms but also using the dual-channel spectrogram, as depicted in Fig. 9. The time signal seen at the top of this figure corresponds to the force applied to the end of a cable for an imposed displacement of constant magnitude and increasing frequency. The frequency spectrogram has various harmonics on the force measurement channel of this imposed oscillation. With the dual channel, adding information from the frequency spectrogram of the displacement channel means we can eliminate any fringes that are present and show only the frequency difference between channels for a given harmonic. In this example, it is clear that there is no frequency shift between the measurement channels for the first harmonic and the phase shows a 1801 shift through the resonance at 0.04 s. When we compare two time–frequency components, we want to know to what extent these components resemble each other. Take note that these components can be present on two channels acquired in parallel or on a single channel acquired at two different instants. It is useful to compare the time plots when these components are nearly identical. However, the least deviation of a vibration mode or the presence of an ambient noise that differs under these components will mask the correspondence in vibration modes that are identical. The energy density spectrogram can then be used to acknowledge the similarity of the modal amplitudes present. The frequency spectrogram allows us to observe whether there are any fluctuations in the frequency from one component to another. Lastly, the
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
137
dual-channel phase spectrogram gives us the opportunity to show the phase correspondence between components. Fig. 10 gives an example where a steel ball is released from a fixed height to impact on the shed of a porcelain insulator. Comparison of the two vibratory responses here reveals the reproduction quality of the measurement. As can be seen, several modes keep the same phase progression pattern for over 500 cycles. This sharpness of the spectrogram opens the way to new vibratory analysis procedures. 7. Conclusion The information contained in the phase of the STFT is easily readable in the time–frequency plane and can be used for diagnostic purposes. Fringe patterns can be seen allowing us to distinguish between weak frequency modulation and an amplitude modulation. The capacity to distinguish these modulations in the visual interpretation opens the way to new diagnostic approaches. The method proposed here can be summarised as aligning the phase with the vertical, according to the frequency axis in the time–frequency plane, then unwrapping the phase horizontally, on the time axis. The derivative of the phase in the time axis yields the frequency spectrogram whereas the integral of this derivative gives the generalised phase spectrogram. The display is restricted to the phase or frequency values calculated under the spectral lobes corresponding to significant amplitude sinusoids; the phase value is not used anywhere except under these lobes. Development of a dual-channel time–frequency analyser is quite straightforward: it involves simply determining the difference between two unwrapped-phase spectrograms or between two frequency spectrograms. There is but one constraint, namely the presence of similar, superimposable time–frequency components in the two spectrograms. We believe we have found a practical way to present the phase information of a signal comprising several sinusoids. Long ignored in favour of presenting the spectral power in the time–frequency plane, the phase information of a STFT can now assume its rightful place and find appropriate application niches. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
D. Gabor, Theory of communication, Journal of IEE 93 (3) (1946) 429–457. J. Ville, The´orie et applications de la notion de signal analytique, Cables & Transmission 2A (1) (1948) 61–74. L. Cohen, Time–frequency distribution—a review, Proceeding IEEE 77 (7) (1989) 941–981. F. Auger, P. Flandrin, Improving the readability of time–frequency and time-scale representations by using the reassignment method, IEEE Transaction on Signal Processing 43 (5) (1995) 1068–1089. K. Kodera, R. Gendrin, C. de Villedary, Analysis of time-varying signals with small BT values, IEEE Transaction on ASSP 26 (1) (1978) 64–76. P. Flandrin, Time–frequency and time-scale, in: Proceedings of the Fourth Annual Acoustics, Speech and Signal Processing Workshop on Spectrum Estimation and Modeling, 1988, pp. 77–80. A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, New Jersey, 1975. W.B. Hershey, J.M. Kim, Impulse response applications to nondestructive testing, International Advances in Nondestructive testing 15 (1990) 289–311. J.M. Tribolet, A new phase unwrapping algorithm, IEEE Transactions on Acoustics, Speech and Signal Proccessing 25 (2) (1977). B. Bhanu, J.H. McClellan, On the computation of the complex spectrum, IEEE Transactions on Acoustics, Speech and Signal Processing 28 (1980) 583–585. R. McGowan, R. Kuc, A direct relation between a signal time series and its unwrapped phase, IEEE Transactions on Acoustics, Speech and Signal Proccessing 30 (1982) 719–726. H. Nashi, Phase unwrapping of digital signals, IEEE Transactions on Acoustics, Speech and Signal Processing 37 (11) (1989) 1693–1702. F. Le´onard, Spectrogramme de phase et spectrogramme de fre´quence, Revue du Traitement du Signal 17 (4) (2001) 269–286. F. LE´onard, J. Lanteigne, S. Lalonde, Y. Turcotte, Vibration behavior of a cracked cantilever beam, in: Proceedings of the 18th International Modal Analysis Conference, San Antonio, Texas, 2000. F. Le´onard, J. Lanteigne, S. Lalonde, Y. Turcotte, Vibration behavior of a cracked cantilever beam, Mechanical Systems and Signal Processing 15 (3) (2001) 259–548. F. Le´onard, Spectrogramme de phase et spectrogramme de fre´quence, Hydro-Que´bec Internal Report ireq-98-178, 1998.
Phase spectrogram and frequency spectrogram as new diagnostic tools Franc- ois Le´onard Institut de recherche d’Hydro-Que´bec, 1800 boul. Lionel-Boulet, Varennes, Que´bec, Canada J3X 1S1 Received 11 April 2005; received in revised form 2 August 2005; accepted 5 August 2005 Available online 23 September 2005
Abstract The power spectrogram has found wide use in diagnostics based on mechanical vibrations. It makes use of the amplitude values of the Short Time Fourier Transform (STFT) but does not consider the information contained in the phase values. In the phase spectrogram, on the other hand, it is precisely the phase values of the STFT that are processed so that they can be presented graphically and are useful. The time–frequency distribution of the phase of a specific component supplies information about the phase modulations around a reference point, which determines both the reference phase and the reference frequency for this component. The frequency spectrogram is calculated from the phase difference between each time slice of the STFT. The frequency spectrogram shows the drift on the instantaneous frequency of each spectral component. These tools have proven to be useful complements to the power spectrogram. In crack detection, for example, they accentuate any frequency drift that occurs in the damped free-decay response of a cantilever beam, which appears to be a function of the crack depth. When identifying a sporadically pulsed vibration source, the ability to show the phase coherence from one impact to the other in time on the spectrogram makes it possible to draw a conclusion about the presence of a stable-vibration source hidden in the background. A dual-channel phase spectrogram is also presented with two application examples. Published by Elsevier Ltd. Keywords: Spectrogram; Phase spectrogram; Frequency spectrogram; STFT; Phase unwrapping; Modulation; Instantaneous frequency; Dual-channel time-frequency distribution; Crack detection
1. Introduction The time–frequency representation was first mentioned in the work of Gabor [1] and Ville [2] in the mid1940s. Gabor represented the signal in the time–frequency plane as a distribution of ‘‘quanta of information’’ whereas Ville was in search of a means of obtaining an energy density. Since then, many efforts have been invested to combine and consolidate the different approaches under Cohen’s class [3]. However, although certain algorithms make use of the phase in the time–frequency plane [4–6], no author has come up with a way to graphically represent the phase in this plane.
E-mail address: [email protected]. 0888-3270/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.ymssp.2005.08.011
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
126
Nomenclature i index of the spectral line m,k index of the time slice mr index of the reference time slice, fixed value fs sampling frequency, fixed value f m;i ðPÞ instantaneous frequency of the wrapped frequency spectrogram fU m;i ðP; mr Þ instantaneous frequency of the unwrapped frequency spectrogram l number of time samples between two consecutive time slices, fixed value P, Pk integer gain, constant value Fm;i ðP; mr Þ instantaneous unwrapped phase related to the reference time slice mr W; Wm;i integer number of phase rotations in the trigonometric circle Sm;i STFT result, amplitude complex value y; ym;i unwrapped phase y; ym;i wrapped phase yG m;i ðP; mr Þ instantaneous phase of the generalised unwrapped phase spectrogram
To start with, we need a time–frequency distribution of the amplitude and phase of the signal. A fair number of algorithms are able to provide this distribution but, for reasons of popularity and simplicity, we take as our starting point the Short Time Fourier Transform (STFT). Usually the term ‘‘spectrogram’’ is reserved for presentation based on a power distribution of the STFT which, on a time signal s(t), is written: Z 1 F s ðt; nÞ ¼ sðtÞ h ðt tÞej2pnt dt, (1) 1
with F s ðt; nÞ denoting the complex distribution function and h ðt tÞ, the conjugate of the spectral window used as a time–frequency kernel. The resulting product F s ðt; nÞ F s ðt; nÞ gives the spectrogram. In this text, we will refer to the ‘‘power spectrogram’’ when we refer to this product, as opposed to the phase spectrogram and frequency spectrogram, which will be introduced later. The STFT is a reversible transform since we need only a surface integration in the time–frequency plane [6] such as Z Z 1 sðtÞ ¼ F s ðt; nÞ hðt tÞej2pnt dt dn (2) 1
in order to find the time signal s(t). On the other hand, it is much more difficult to reconstitute the signal exactly using just the power spectrogram [7] because the latter does not contain all the information needed to reconstruct the original time signal: what is missing is the phase information contained in the ‘‘phase’’ argument of the STFT. It should be noted that interference between the components present can be observed on the power spectrogram. Actually, it is the phase information that this interference contains that allows us to estimate the time signal from the power spectrogram. Although the latter does contain a certain amount of phase information, it is nevertheless incapable of providing an adequate illustration of this phase. This paper explains how the phase spectrogram, using the argument of the STFT, provides access to a source of information which completes that offered by the power spectrogram. The information in a phase spectrogram can take different forms. The easiest to understand intuitively is the phase variation of a sine wave compared to a reference sine wave, the latter being the same sine wave taken at a given time and not having changed its frequency since. On the other hand, in deriving this phase, we generate a frequency spectrogram that gives us a different perception of the instantaneous frequency.
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
127
2. Phase unwrapping Before we tackle the algorithms of the phase and frequency spectrogram, we should clarify the mathematics involved in phase unwrapping. When we extract a phase value from a complex amplitude in a Fourier domain, we obtain its projection y in the trigonometric circle. The unwrapped phase value y ¼ y þW 2p
with W 2 Z
(3)
differs by W rotations in the trigonometric circle; we know the position of the phase on the circle but not the number of rotations. Phase unwrapping amounts to finding this number, which is done progressively by tracing round the trigonometric circle repetitively, usually starting from a chosen reference point. The unwrapped value y only has any meaning if we can associate it with a reference that has physical significance. For example, in cepstral and modal analysis, phase unwrapping is done on the frequency axis with y ¼ 0 at 0 Hz as the reference point. This is unavoidable because the d.c. component has a null imaginary value and, also, because a time delay applied to the 0 Hz frequency does not result in a phase shift. On the other hand, the choice of reference is more arbitrary for the phase spectrogram where phase unwrapping follows the time axis in the time-frequency plane. The phase here never seems smooth without appropriate unwrapping [8]. Several publications deal with phase unwrapping on the frequency axis [7,9–12] but no author presents the phase unwrapping on the time axis in the time–frequency plane. However, whether it be the frequency or the time axis, phase unwrapping consists in selecting the shortest path on the trigonometric circle. Usually, the unwrapping can only be achieved if the phase does not vary by more than p radians between two successive time slices. The problem is that two successive time slices of a STFT can be several time samples apart. Let l be the step size in number of samples separating two adjacent time slices of the STFT. From Nyquist criterion, the highest frequency that sampling without spectral aliasing can withstand is equal to a half sampling frequency, which produces an increase of p radians between each time sample for this component. The latter will therefore have a phase increase of lp radians between two successive time slices of the STFT. It seems, therefore, that the phase increase exceeds p for frequencies higher than f s =ð2lÞ, where f s is the sampling frequency. At first sight, phase unwrapping appears difficult if not impossible on the time axis of a STFT but, as we will see later, the difficulty can be overcome. In this text, we use underscored notations to mean that the y phase y is contained within the trigonometric circle. The underscored operator n ¼ mod2p fn þ pg p
(4)
therefore removes the 2p multiple value in the operand. We should also make it clear that the phase ym;i comes from the trigonometric function atan2 (*,*) applied to the result of the STFT such that ym;i ¼ atan 2ðReðSm;i Þ; ImðS m;i ÞÞ
(5)
with ym;i 2 p; p, where m is the index of the time slice and i that of the spectral line. In calculating the STFT, application of a spectral window comprising N samples will give N/2 spectral lines such that i 2 f0; . . . ; N=2 1g is the spectral-line index. The phase-unwrapping algorithm in time is recursive and usually starts with the first time slice. A simple form of this algorithm, starting with the first sample, consists in writing ynm;i ¼ ym;i þWm;i 2p
with Wm;i 2 Z,
(6)
where Wm corresponds to the number of rotations or turns completed in the trigonometric circle. The integer Wm;i is established recursively such that 8 W for po ym;i ym1;i pp > < m1;i 1 þ W for ym;i ym1;i p p m1;i and W1;i ¼ 0, Wm;i ¼ (7) > : Wm1;i -1 forpo ym;i ym1;i so that any phase shift with an amplitude exceeding p can be attributed to the fact that the boundary p=p of the trigonometric circle has been crossed. If the phase sampling condition is respected, then ym;i ¼ ym;i . For
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
128
simplification, we write ym;i ¼ unwrap ym;i
(8)
to express the unwrapping of the phase distribution ym;i . The operator unwraps the phase without worrying about there being possibly n extra 2p radians between two consecutive samples. Lastly, in order to refer the unwrapped phase Wm;i to the reference time slice mr , all we need to do is to subtract the latter phase such that ym;i ymr ;i gives the unwrapped phase in relation to the reference slice. 3. Frequency spectrogram The wrapped frequency spectrogram f f m;i ðPÞ s P ym;i ym1;i 2plP
(9)
is obtained from the phase difference between successive time slices of the STFT where P 2 Z. A mathematical demonstration has been developed and can be found in [13]. The gain P multiplies the sensitivity expressed as the fringes per Hz in such a way that fs fs ; f m;i ðPÞ 2 (10) 2lP 2lP defines the frequency axis scaling. All we need to do is display the frequency on a scale of colours where the two extremes f s =2lP and f s =2lP have the same colour, in order to watch the frequency develop without any break in the progression. The gain P is adjusted so that the progression can be monitored and we can easily count the fringes on the screen. The sensitivity f s =lP is expressed in Hz/fringe. that this definition of the frequency spectrogram depends largely on the equality We should emphasise
P ym;i ym1;i ¼ P ym;i ym1;i þW2p
for integers W and P. With an integer value gain, we amplify the
phase difference between two time slices without any disturbance from an error in the estimation of a phase of ‘‘W’’ turns in the trigonometric circle: we no longer need to know the number of rotations of the unwrapped phase value before applying the gain. Once an integer value gain has been applied, we can unwrap the phase. In fact, the unwrapped frequency spectrogram fs fU unwrap P ym;i ym1;i unwrap P ymr ;i ymr 1;i m;i ðP; mr Þ 2plP for m ¼ 2; 3; 4; . . . ; M ð11Þ allows us to display the frequency relative to the reference slice mr on a continuous scale. The frequency presents the same topological characteristic as the phase in the trigonometric circle, which means that it can be illustrated as superposed on a cyclic scale (Eq. (9)) or unwrapped on a conventional scale (Eq. (11)). The term P ym;i ym1;i in Eqs. (9) and (11) can be replaced by NP X
Pk ymþk;i ymþk1;i
(12)
k¼1
for weighting the computation of the phase derivative. This smoothing distributes the gain P according to an approximately Gaussian distribution of smaller weight PI so that P¼
NP X k¼1
Pk
(13)
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
129
Fig. 1. Power spectrogram and frequency spectrogram of the frequency-modulated signal cosð2p 200 Hz t þ 50 expð2tÞ sinð2p 4 Hz tÞÞ sampled at 2000 samples/s (N ¼ 256 for a 4-term Blackman–Harris window, spectral interpolation 2 , l ¼ 3, P ¼ 33).
with integers Pi. The advantage of the weighting is that we obtain smoothing with a Gaussian distribution whereas increasing the distance l between slices amounts to having a smoothing with a rectangular window where P ¼ Pi N p . The phase values from the STFT fill a matrix measuring M I, where M is the total number of time slices, I ¼ N=2 corresponds to the number of spectral lines. The term ymþk;i ymþk1;i then becomes the simple subtraction of a matrix from the same matrix offset by a time slice, namely M I operations. Applying the operator n to Pk ymþk;i ymþk1;i , like the operator unwrap, requires very little computation time, namely a
multiple of MI operations. On the whole, the computation effort increases in proportion to the dimensions of the time–frequency matrix. Fig. 1 depicts a power spectrogram with, below it, the unwrapped frequency spectrogram of a numerically generated cosine, which is frequency-modulated by a damped sine wave. Since the instantaneous frequency of a sine wave has no significance other than below the main lobe of the corresponding spectral component, only the result under the main lobe is displayed. In order to do this, we fill in the corresponding power spectrogram valleys in black to drown the secondary components with the lowest amplitudes. The height of the filling is adjusted manually by the user, depending on the amount of information to be underlined. The LabViewTM code and runtime is available on request. In computing the STFT, we recommend applying a spectral window having high side lobe rejection. A typical lobe comprises 5–7 spectral lines, which gives a frequency plateau of the same width if the lobe emerges sufficiently from the filling level. Obviously we can enhance the graphic resolution of the plateau by increasing the number of spectral lines beneath it. We need to only perform a spectral interpolation, for example by increasing the width of the spectral window with zero padding. We should clarify that the number of samples N here is the number before zero padding. For example, an interpolation of 2 with N ¼ 256 means that 256 zeros are added and that the TRF will apply to 512 samples, which gives 256 spectral lines. The modulation of phase or frequency generates a special fringe pattern: the fringe appears perpendicular to the path of the sine wave in the time–frequency plane, as illustrated in Fig. 1. Fig. 2, by contrast, shows that a
ARTICLE IN PRESS 130
F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
Fig. 2. Frequency spectrogram of a numerically generated sine wave, modulated in amplitude and sampled at 2 ks/s. (2 Hz/fringe, N ¼ 256 for a 4-term Blackman–Harris window, spectral interpolation 4 , l ¼ 3, P ¼ 333).
Fig. 3. Two different analyses of the same part of signal illustrated in Fig. 1 showing two apparent angles of the fringe alignment with the path of the sine wave.
modulation in amplitude has a different pattern: fringe waves in parallel to the path of the sine wave. Here a 4 spectral interpolation is required to be able to display this fringe across the entire width of the plateau. Caution is necessary, however, when we come to interpret the fringe patterns. As shown in Fig. 3, for a same frequency resolution the perspective given by two parameters sets changes the apparent angle between the fringe pattern and the path followed by the sine wave. Lastly, a recent paper [13] has shown that the variance on the frequency is minimal at the centre of the frequency plateau and is of the same order of magnitude as the variance obtained using the spectral interpolation algorithms (IFFT), namely near the Crame´r–Rao bound on the frequency estimate.
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
131
Fig. 4. Phase spectrogram of a sine wave with a phase modulation of 7p radians, i.e. the signal cosð2p 200 Hz t þ p cosð2p 2Hz tÞÞ with f s ¼ 1 s/s (N ¼ 128 samples, 2 spectral interpolation, l ¼ 5 and P ¼ 1).
4. Phase spectrogram The phase difference is equal to 2p
lP U f ðP; mr Þ þ C f s m;i
(14)
between two successive time slices where C is an unknown corresponding to the number of complete turns in the trigonometric circle between two successive slices of the STFT. The mean frequency of the component thus determines the magnitude of C. The information that C contains is of little use because the focus of our interest is the phase variations around the mean increment in the phase increase, not the slope of the phase increase due to the mean frequency of the component. Therefore, if we remove C from the computation, numerical integration of the type 8 0; m¼1 > < m X Fm;i ðP; mr Þ ¼ 2p l (15) fU m ¼ 2; 3; . . . ; M m;j ðP; mr Þ; > : fs j¼2 yields an unwrapped-phase value starting from the first time slice. In Eq. (15), the phase calculation is referenced to both the first slice for calculating the phase and slice mr for calculating the reference frequency. The latter frequency determines the Dy to be removed between each slice so as to eliminate the constant term of the phase increase. The generalised phase spectrogram defined as yG m;i ðP; mr Þ Fm;i ðP; mr Þ Fmr ;i ðP; mr Þ
(16)
allows us to refer both the phase calculation and the frequency calculation simultaneously to the mr th time slice. For a display with a succession of fringes, we use this presentation yG m;i ðP; mr Þ but we could also write yG m;i ¼ Fm;i ðP; mr Þ
m1 FM1;i ðP; mr Þ M 1
(17)
if we want to impose a null phase value on the two extremities of the measurement. If the distance l and the distribution Pi have an influence on the variance in the frequency spectrogram, they have none on the phase spectrogram [13]. In the latter case, the distance l affects the graphic definition on the time axis while the weight P restrains the phase to the interval ]p=P; p=P], superposing the distribution of phase values for no useful purpose. The weight P is therefore set at unity for calculating the phase spectrogram and the feature proposed at Eq. (12) is not used. The numerical example in Fig. 4 shows the presence of a uniform phase plateau across the width of the main lobe of the component. Also, as in the frequency spectrogram with a phase modulation, the phase spectrogram
ARTICLE IN PRESS 132
F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
Fig. 5. Phase spectrograms of the free response of beams excited by a mechanical impact. The vertical line in each figure shows the position of reference mr (2p radians/fringe, with f s ¼ 25:6 ks=s, N ¼ 256 samples, 8 spectral interpolation, l ¼ 5 and P ¼ 1).
of a phase modulation shows perpendicular fringes in the path of the sine wave while an amplitude modulation shows fringes in the opposite direction on the phase spectrogram (not illustrated. See [13]). Lastly, the use of coloured fringes makes it far easier to read the graphics. With monochrome displays we had to use a different shade of grey for the extremities of the display range to determine the direction of the phase increase. With colours we can have colour continuity when we cross the boundary from p to p. All we need to do is to use three or more colours in the range concerned in order to determine the direction of the phase or frequency progression by counting the fringes.
5. Applications to monitoring and diagnostics The phase and frequency spectrograms were developed in the framework of a crack detection project to reveal slight changes in the frequency of the free damping response of excited modes [14,15]. These frequency changes result in phase shifts, as illustrated in Fig. 5. The modes show a constant phase in the free response of the uncracked beam while a phase drift appears with the other three cracked beams and increases with the surface area of the crack, occurring mainly when the vibration amplitude is no longer enough to completely open the crack. To open a fatigue crack in steel we need to exert a force or it will remain tightly closed. Three vibration states succeed each other on the damped free response of a steel beam when the amplitude of the impact is high enough to open the crack straightaway. Fig. 6 shows where these states are located on a frequency spectrogram of the damped free response of the steel beam that is most badly cracked. In the first state, the crack opens completely and closes alternately with the vibration beat. Since there is less rigidity when the crack is open, the mean frequency of the first vibration is lower. In the second state, the vibration amplitude allows only partial opening of the crack; it is in this second state that the damping reaches a maximum and the frequency varies most. Lastly there comes a point where the vibration amplitude is clearly not capable of
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
133
Fig. 6. Frequency spectrogram of the free response of a badly cracked beam excited by a mechanical impact (f s ¼ 25:6 ks=s).
Fig. 7. Frequency spectrogram of an acoustic signal measured remotely by a microphone (above) and at the emission source by an accelerometer (below).
opening the crack and the third state begins, characterised by the highest modal frequencies and less damping than in the first two states. In the previous example, the amplitude of the transient signal is multiplied by an exponential gain increasing with time in order to keep the amplitude of the vibration components near a constant value. This not only gives a more uniform plateau width but also reduces the effect of free damping on the amplitude modulation. The sound wave propagating in a moving medium, unlike light, does not travel at a constant speed in all directions for all observers. The wind, for example, can change the speed of the propagation of a sound wave and thereby alter the time delay between emission and reception, so here we have phase modulation. To this we have to add a slight amplitude modulation, which is difficult to observe on a power spectrogram. Here again, the frequency spectrogram appears as a useful tool for displaying and quantifying the physical phenomenon as shown in Fig. 7, which presents the vibration of one of the dominant modes of one of the sheds of a high-voltage porcelain insulator in a distribution substation. These modes had been excited by the impact of a small metal ball. The wave distortion here corresponds to wind speed fluctuations of around 1 m/s between the source and the measurement. This phenomenon cannot be seen in the lab measurements but it made us take another look at the algorithms developed in the laboratory for detecting cracked insulators. In the area of vibration diagnostics, we are sometimes called upon to determine whether the intermittent bursts observed on a spectrogram come from a mode excited randomly or from some distant excitation with a stable frequency but transmitted sporadically. Fig. 8 shows the phase coherence from one burst to another in
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
134
Fig. 8. Power spectrogram and phase spectrogram (2p radians/fringe) of an acceleration signal recorded in the cabin of an airline pilot. This measurement was kindly provided by Rene´ Archambault IMS. (f s ¼ 24 Ks=s, N ¼ 128, l ¼ 3, P ¼ 1).
time, revealing the presence of a stable-vibration source that is observed sporadically. Vibration coming from a resonance excited by random shocks would not have this phase continuity. 6. Dual-channel applications Dual-channel applications [16] offer the opportunity to observe the frequency gap or the phase gap between different channels or measurements in the time–frequency plane. The required condition is that the signals have similar fingerprints in the time–frequency plane so that if we want to compare sine waves present on two channels or measurements, these must follow the same trajectories in the time– frequency plane. If need be, we can have a time translation in order to obtain the superposition of the fingerprints. We cannot resort to subtraction of the phase argument of two STFTs to observe the dual-channel phase because this phase difference does not fulfill the phase sampling condition. For example, if we cross the boundary p=p on a channel at a different time than on another channel, a phase jump can occur which can be interpreted in two ways: either as one of the channels crossing the boundary or as a real difference by more than a value p between the two channels without any crossing of this boundary in the trigonometric circle. Unwrapping of the phase on each channel must be done before phase subtraction or else information will be lost, resulting in ambiguous interpretation of the phase increases. Subtraction of two unwrapped phase spectrograms allows us to present the evolution of the phase between sinusoids from different measurements that are superposed in the time–frequency plane. To be more precise, the dual-channel phase spectrogram here gives the pattern of the relative difference in the phase evolution of two juxtaposed sinusoids and not the phase difference itself. If we need to show the phase difference between two sinusoids, the process must include the phase of the reference spectral slice mr common to each channel (A and B) by adding the phase vector P ymr ;i prior to subtracting the phase spectrograms so that yAB dualchannel
1 A y yB , P
(18)
where yA ¼ B
y ¼
P yA mr ;i þ P yBmr ;i
þ
m X
! mr X A A unwrap P yA unwrap P yA , m;i ym1;i m;i ym1;i
j¼2
j¼2
m X
! mr X B B B B unwrap P ym;i ym1;i unwrap P ym;i ym1;i .
j¼2
j¼2
Thus, we obtain a phase difference between components, which is more practical for such purposes as modal analysis of an unsteady frequency response, for example.
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
135
As for the dual-channel frequency spectrogram, this could be generated by the difference A B f AB dualchannel f m;i ðPÞ f m;i ðPÞ
(19)
between the phase derivatives calculated on two STFTs each of which corresponds to a measurement channel. This wrapped instantaneous frequency gap appears significant where the components are superposed. By contrast, the unwrapped frequency gap Unwrapped Unwrapped f AB fA ðP; mr Þ f Bm;i Unwrapped ðP; mr Þ m;i dualchannel
(20)
gives a significant result for the frequency lines included in the frequency plateau of the component at the intersection of the time slice mr .
Fig. 9. Time trace, single- and dual-channel frequency spectrograms and dual-channel phase spectrogram of a frequency sweep test with an imposed displacement on a HV cable.
ARTICLE IN PRESS 136
F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
Fig. 10. Dual-channel frequency spectrogram comparing the phase evolution between the responses to two impacts on a porcelain insulator observed by the same channel.
The concept of the dual channel for phase or frequency spectrograms is wider than the same concept for a dual-channel frequency analyser. In fact, in the latter case, we are comparing the phase frequency distribution of two channels sampled simultaneously, for one time slice. With the spectrogram, however, we have an extra dimension: time. The dual-channel frequency or phase spectrogram allows us to observe the temporal evolution of a frequency or phase relation of time–frequency components. All that is required is that the components compared share the same spectral position and some coherence between them (in other words, a linear physical relationship). The components compared may occupy different temporal positions. For example, the time–frequency components detected by two remote microphones must be juxtaposed in time: note that any positioning error here causes interference, which reduces the surface area where there is coherence between the components. The components may even come from different tests. For example, we can compare two transients present on the same channel to observe whether they have the same frequency or even phase modulation characteristics. While today’s modal analysis tools are not adequate for us to tackle unsteady frequency modes, the dualchannel phase spectrogram makes it possible to show the phase relation between channels without accounting for frequency fluctuations, as long as the latter do not move by more than three spectral lines. In addition, the dual-channel phase spectrogram also allows us to observe the phase behaviour of a mode between the fixed and mobile parts on machines. We can compare two time–frequency components present on channels sampled in parallel using two power spectrograms but also using the dual-channel spectrogram, as depicted in Fig. 9. The time signal seen at the top of this figure corresponds to the force applied to the end of a cable for an imposed displacement of constant magnitude and increasing frequency. The frequency spectrogram has various harmonics on the force measurement channel of this imposed oscillation. With the dual channel, adding information from the frequency spectrogram of the displacement channel means we can eliminate any fringes that are present and show only the frequency difference between channels for a given harmonic. In this example, it is clear that there is no frequency shift between the measurement channels for the first harmonic and the phase shows a 1801 shift through the resonance at 0.04 s. When we compare two time–frequency components, we want to know to what extent these components resemble each other. Take note that these components can be present on two channels acquired in parallel or on a single channel acquired at two different instants. It is useful to compare the time plots when these components are nearly identical. However, the least deviation of a vibration mode or the presence of an ambient noise that differs under these components will mask the correspondence in vibration modes that are identical. The energy density spectrogram can then be used to acknowledge the similarity of the modal amplitudes present. The frequency spectrogram allows us to observe whether there are any fluctuations in the frequency from one component to another. Lastly, the
ARTICLE IN PRESS F. Le´onard / Mechanical Systems and Signal Processing 21 (2007) 125–137
137
dual-channel phase spectrogram gives us the opportunity to show the phase correspondence between components. Fig. 10 gives an example where a steel ball is released from a fixed height to impact on the shed of a porcelain insulator. Comparison of the two vibratory responses here reveals the reproduction quality of the measurement. As can be seen, several modes keep the same phase progression pattern for over 500 cycles. This sharpness of the spectrogram opens the way to new vibratory analysis procedures. 7. Conclusion The information contained in the phase of the STFT is easily readable in the time–frequency plane and can be used for diagnostic purposes. Fringe patterns can be seen allowing us to distinguish between weak frequency modulation and an amplitude modulation. The capacity to distinguish these modulations in the visual interpretation opens the way to new diagnostic approaches. The method proposed here can be summarised as aligning the phase with the vertical, according to the frequency axis in the time–frequency plane, then unwrapping the phase horizontally, on the time axis. The derivative of the phase in the time axis yields the frequency spectrogram whereas the integral of this derivative gives the generalised phase spectrogram. The display is restricted to the phase or frequency values calculated under the spectral lobes corresponding to significant amplitude sinusoids; the phase value is not used anywhere except under these lobes. Development of a dual-channel time–frequency analyser is quite straightforward: it involves simply determining the difference between two unwrapped-phase spectrograms or between two frequency spectrograms. There is but one constraint, namely the presence of similar, superimposable time–frequency components in the two spectrograms. We believe we have found a practical way to present the phase information of a signal comprising several sinusoids. Long ignored in favour of presenting the spectral power in the time–frequency plane, the phase information of a STFT can now assume its rightful place and find appropriate application niches. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
D. Gabor, Theory of communication, Journal of IEE 93 (3) (1946) 429–457. J. Ville, The´orie et applications de la notion de signal analytique, Cables & Transmission 2A (1) (1948) 61–74. L. Cohen, Time–frequency distribution—a review, Proceeding IEEE 77 (7) (1989) 941–981. F. Auger, P. Flandrin, Improving the readability of time–frequency and time-scale representations by using the reassignment method, IEEE Transaction on Signal Processing 43 (5) (1995) 1068–1089. K. Kodera, R. Gendrin, C. de Villedary, Analysis of time-varying signals with small BT values, IEEE Transaction on ASSP 26 (1) (1978) 64–76. P. Flandrin, Time–frequency and time-scale, in: Proceedings of the Fourth Annual Acoustics, Speech and Signal Processing Workshop on Spectrum Estimation and Modeling, 1988, pp. 77–80. A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, New Jersey, 1975. W.B. Hershey, J.M. Kim, Impulse response applications to nondestructive testing, International Advances in Nondestructive testing 15 (1990) 289–311. J.M. Tribolet, A new phase unwrapping algorithm, IEEE Transactions on Acoustics, Speech and Signal Proccessing 25 (2) (1977). B. Bhanu, J.H. McClellan, On the computation of the complex spectrum, IEEE Transactions on Acoustics, Speech and Signal Processing 28 (1980) 583–585. R. McGowan, R. Kuc, A direct relation between a signal time series and its unwrapped phase, IEEE Transactions on Acoustics, Speech and Signal Proccessing 30 (1982) 719–726. H. Nashi, Phase unwrapping of digital signals, IEEE Transactions on Acoustics, Speech and Signal Processing 37 (11) (1989) 1693–1702. F. Le´onard, Spectrogramme de phase et spectrogramme de fre´quence, Revue du Traitement du Signal 17 (4) (2001) 269–286. F. LE´onard, J. Lanteigne, S. Lalonde, Y. Turcotte, Vibration behavior of a cracked cantilever beam, in: Proceedings of the 18th International Modal Analysis Conference, San Antonio, Texas, 2000. F. Le´onard, J. Lanteigne, S. Lalonde, Y. Turcotte, Vibration behavior of a cracked cantilever beam, Mechanical Systems and Signal Processing 15 (3) (2001) 259–548. F. Le´onard, Spectrogramme de phase et spectrogramme de fre´quence, Hydro-Que´bec Internal Report ireq-98-178, 1998.