Time-Frequency Analysis of Clinical Percussion Signals Using Matrix Pencil Method

This paper discusses time-frequency analysis of clinical percussion signals produced by tapping over human chest or abdomen with a neurological hammer and recorded with an air microphone. The analysis of short, highly damped percussion signals using conventional time-frequency distributions (TFDs) meets certain difficulties, such as poor time-frequency localization, cross terms, and masking of the lower energy features by the higher energy ones. The above shortcomings lead to inaccurate and ambiguous representation of the signal behavior in the time-frequency plane. This work describes an attempt to construct a TF representation specifically tailored to clinical percussion signals to achieve better resolution of individual components corresponding to physical oscillation modes. Matrix Pencil Method (MPM) is used to decompose the signal into a set of exponentially damped sinusoids, which are then plotted in the time-frequency plane. Such representation provides better visualization of the signal structure than the commonly used frequency-amplitude plots and facilitates tracking subtle changes in the signal for diagnostic purposes. The performance of our approach has been verified on both ideal and real percussion signals.TheMPM-based time-frequency analysis appears to be a better choice for clinical percussion signals than conventional TFDs, while its ability to visualize damping has immediate practical applications.


Introduction
Clinical percussion is a centuries-old bedside diagnostic technique used to identify various conditions of the thorax and abdomen in health and disease [1,2].It is a method of eliciting sounds by striking body parts with fingertips or a percussion hammer [1,2].The purpose of percussion is to provide the physician with qualitative information about the size, consistency, and borders of the vital organ or pathology under the percussed area.Trained physicians are able to recognize many different kinds of percussion sounds, the main three of which are widely known as "resonant, " "tympanic, " and "dull" [1,2].Resonant sounds are low-pitch, hollow sounds heard over normal lung tissue.Abnormal lungs may be hyperresonant, dull, or stony dull.Dullness is expected over the liver and over the heart [2].Typical tympanic sound is drum-like and is usually heard over the air cavities in the abdomen as well as over pneumothorax [2,3].In the latter case, percussion response of the chest changes from resonant to tympanic due to acoustic impedance mismatch at the air cavity boundary [1].Due to disrupted acoustic coupling between the chest wall and the lung parenchyma (normally provided by the wet contact between visceral and parietal pleurae), the transfer of the acoustical energy to the lungs (where it subsequently dissipates) becomes less efficient [1].Hence the energy dissipates less rapidly and stays within the chest wall for a longer period of time, leading to the extended signal duration and decreased damping [1].The goal of this work is to come up with a practical way to visualize this effect of changed damping of certain signal components in the time-frequency plane.
The majority of preceding work attempting the objective classification of percussion sounds is based on the timedomain and Fourier spectral analysis [2,4].Spectral peaks correspond to individual oscillation modes, while their width is a direct measure of damping.When closely spaced oscillation modes have high damping, their broad spectral peaks overlap and often cannot be resolved with required accuracy.Time-frequency analysis combines the advantages of the pure time-domain and spectral methods and brings important additional benefits, such as the ability to trace the evolution of the spectral content with time, thus providing a convenient tool to dissect, analyze, and interpret signals [5,6].
A clinical percussion signal represents free transient response of the human body to a brief mechanical impact.Due to damping, the amplitude of these signals changes with time, and hence they can be considered nonstationary in the broad sense of this definition [7].Since the impact usually excites oscillations of more than one anatomical subsystem, percussion signals usually represent more than one independent oscillation mode with different frequencies and hence can be considered multicomponent, again, in the broad sense of this term [7].Representing such signals in the time-frequency domain is not always easy due to their short duration and possible coexistence of multiple frequency components with very different amplitudes.When dealing with multicomponent signals, most conventional TFDs demonstrate poor time-frequency localization, lack of accuracy, insufficient cross terms suppression, and masking of the lower energy transient components by the higher energy ones [3,8,9].Such drawbacks severely distort the time-frequency portrait of percussion signals, limiting practical applications.The most extensively exploited linear time-frequency analysis methods are the Short Time Fourier Transform (STFT) and Wavelet Transforms [5,10].Although STFT is immune to cross terms, its biggest limitation is the tradeoff between the time and frequency resolution due to the Heisenberg uncertainty principle [8,11].In case of the Wavelet Transform, the frequency resolution decreases with increasing frequency [10].If more than one frequency component exists in the signal, satisfactory resolution for all modes cannot be obtained.This is an inherent characteristic of the wavelet analysis [10].On the other hand, the quadratic (Cohen's class) time-frequency distributions are based on estimation of instantaneous energy using bilinear operation on the signal [12].Among the most commonly used timefrequency distributions from Cohen's class are Wigner-Ville Distribution (WVD), Choi-Williams Distribution (CWD), and Zhao-Atlas Marks (ZAM) Distribution [5,6,8].WVD has excellent time-frequency resolution due to the absence of averaging over any finite time interval [11,13].While performing well on the single-component signals, WVD's performance degrades dramatically as the number of components increases, because of the spurious peaks arising from the cross terms and confusing the visual interpretation of its timefrequency spectrum [8,14].CWD is an improvement of the WVD intended to suppress the cross terms by multiplying it by a kernel function [6].Unfortunately, increased suppression of the cross terms invariably leads to smearing or loss of resolution of the autoterms in the time-frequency plane.
Another possible way to construct a time-frequency representation of a signal is by representing it as a linear superposition of elementary functions with well-behaved TFDs and then assembling such elementary TFDs on a compound plot.The majority of methods for performing linear signal decomposition involve overcomplete waveform dictionaries.By selecting the optimal (in certain sense) set of available waveforms from the dictionary, one can obtain a sparse model of the signal.In this work, we chose to represent clinical percussion signals as a linear combination of exponentially damped sinusoids (EDS).Such decomposition seems compatible with the underlying physics.Indeed, the percussed human body represents a passive viscoelastic system with multiple degrees of freedom.A percussion impact excites several such degrees at once, and then the oscillations fade out with their respective rates.If the impact is weak, the transient response can be considered linear, and the damping of all normal modes can be considered exponential [4].
In this work, the decomposition of medical percussion signals into damped harmonics was accomplished using Matrix Pencil Method (MPM).The MPM was chosen because it is reportedly the best performing of all damped harmonic analysis methods, especially in the presence of noise.The MPM decomposes the signal into a number of EDS components [15], detecting the exact location of the dominant frequencies and accurately distinguishing closely spaced frequency components [15,16].Each EDS is fully characterized by its amplitude, frequency, damping, and initial phase.The first two parameters can be visualized in the frequency domain as vertical line segments of certain height (amplitude).This representation allows for direct comparison with the Fourier amplitude spectrum [17].However, without seeing the damping parameter it is difficult to judge the relative contribution of individual EDS to the signal.Visualizing the damping parameter would help better capture the behavior of each EDS and show their time evolution in the time-frequency plane.

Materials and Methods
The percussion signals analyzed in this paper were collected from normal volunteers at the Detroit Medical Center by trained medical personnel.The testing of human subjects was performed according to Protocol # 0710005340 (Portable Pulmonary Injury Diagnostic Device) approved by the Human Investigation Committee for the Wayne State University International Review Board (M1) for the period from 25 September 2008 to 24 September 2009.The signals were produced by gently tapping with a neurological hammer over a mediator plate (plessimeter) placed on the subject's chest or abdomen.The signals were received with a tripod-based omnidirectional electret condenser air microphone.In each test, the microphone was placed 150-300 mm away from the percussion spot.Then the signals were amplified and digitized at 48 KHz sampling rate using a 24-bit computer sound card.Additional descriptions of our testing setup can be found in [4,17,18].Two examples of percussion signal waveforms are shown in Figure 1.The "resonant" signal (Figure 1(a)) was acquired over the upper chest in the subclavicular area and the "tympanic" signal (Figure 1(b)) was acquired over the left portion of the abdomen.The signals recorded in our experiments are quite repeatable, hardware-independent, and, in fact, similar in shape to those recorded by Murray and Neilson in 1975 [2].The amplitude spectra of both signals are also shown in Figure 1.Although the fine spectral structure may differ significantly from signal to signal, in general, "tympanic" spectra tend to be narrower and with fewer major subpeaks than "resonant" ones.

Damped Harmonic Model
In this model, the observed percussion response () = ()+ () is represented as a sum of noise () and the actual signal composed of  EDS with frequencies   , amplitudes   , initial phases   , and damping factors   (total 4 parameters).In the discrete case, Here,  = 0, . . ., −1,  = 1, . . ., ,   is the sampling period,   =    −  are complex amplitudes, and   =  (2  −  )  are signal poles [16,19].Decomposition (1) is carried out using MPM, which is described in more detail in [15][16][17][18][19][20].One of the major MPM advantages is that it does not require any prior knowledge or guessing about the EDS present in the signal.Another advantage is the built-in signal denoising process.After performing singular value decomposition (SVD) of the Hankel data matrix constructed from signal samples, the individual singular values   are arranged on the main diagonal of the resulting matrix in the descending order:  1 ≥  2 ≥ ⋅ ⋅ ⋅ ≥  min .Some noise can be filtered out by discarding singular values   smaller than the empirically chosen threshold , for example, those for which (  / max ) <  (in this work  = 10 −3 ) [15].This built-in thresholding process is very useful, but it does not always completely eradicate the noise.The remaining noise can be further reduced using physical reasons.Indeed, in typical situations, some of the EDS components produced by the MPM may have too low amplitude, zero or negative damping, or unrealistically high frequency.Such EDS can be discarded, as they are artificially introduced by the MPM to better fit the noise or discontinuities in the signal.If the noise level is too high, the number of such artificial EDS can grow uncontrollably, as the signal no longer obeys the damped harmonic model.Therefore, prior to applying MPM, it is necessary to estimate the Signal to Noise Ratio (SNR), which in clinical percussion signals is typically not constant.It is highest at the onset and then decreases as the signal fades out due to damping.A detailed analysis of the effects of the noise on the MPM performance can be found in [17].The most optimal conditions for MPM application to clinical percussion signals estimated in [17] suggest that the SNR at the end of the signal window should not exceed 20 dB. 1 with a 50-1000 Hz bandpass Butterworth filter and digitally finding the approximate location of the signal onset, 1500 samples (∼30 ms duration) of each signal were processed with the MPM algorithm without downsampling.The Butterworth filter was chosen as it has a reasonably sharp frequency response but does not induce substantial passband ripple.The four parameters of all resulting EDS (only those with positive frequency and damping values were selected) are listed in Tables 1 and 2. The EDS are sorted by the normalized amplitude.

Comparison of the Conventional and MPM-Based TFD.
The time-frequency analysis of resonant and tympanic clinical percussion signals of Figure 1 is presented in Figure 2. Figures 2(a) and 2(b) show the results of applying WVD and CWD (with the value of the cross term suppression parameter  = 0.8).These two distributions were chosen because they are the most popular quadratic TFDs and because one of them uses a special kernel function for reducing the cross terms and the other one does not.WWD (Figure 2(a)) is the most blurred and the most difficult to interpret of all TFDs presented.Multiple cross term artifacts (not suppressed by the WVD) indirectly confirm the multicomponent character of clinical percussion signals.
The CWD (Figure 2(b)) appears to be much cleaner than WVD.Beyond the blurred 0-6 ms interval for the resonant signal and 0-12 ms for the tympanic signal, CWD shows several components (three for the resonant signal and two for the tympanic signal) with different frequencies, which do not change with time.The amplitudes of these components appear to decay with time.The frequency values of these components are in qualitative agreement with the frequencies of the high-amplitude EDS from Tables 1 and 2. However, not all of the EDS can be identified from the CWD plot.For example, the 299.1 Hz mode of the resonant signal does not extend beyond the blurred 0-6 ms interval due to relatively high damping and therefore cannot be detected.The blurring within these time intervals can be explained by the high intensity of cross terms, which CWD suppresses at the expense of the resolution.The frequency resolution of the STFT, another conventional non-model-based TFD, can be adjusted by choosing appropriate length of the sliding window, apodization, zero padding, and other common Fourier domain techniques.Figure 2(c) shows the results of applying a reassigned spectrogram (STFT) with a 1400-point apodization window (almost as long as the signal itself).This method is able to resolve same frequency components as CWD.Also, similar to CWD, the separation is not satisfactory in the beginning ( < 5-9 ms).This can be explained by noting that the spectrogram is essentially a Fourier transform that fits each damped oscillation mode with infinite sine waves.Hence, each damped mode appears to be broadband and contains multiple Fourier harmonics at the same time, making the regular STFT (not shown here) highly blurred and not practically useful.The reassignment procedure improves the resolution of the spectrogram by artificially highlighting its ridges.This helps visualize the Fourier bandwidth of each EDS and even single out the main three frequency components.However, even after reassignment, the STFT of percussion signals still contains too many artifacts to be practically useful.Figure 2(d) shows the EDS-based time-frequency representation of the chest and abdominal percussion signals.Since the EDS decomposition is performed using MPM, this representation is named MPM TFR.Each EDS (defined by its frequency , amplitude , damping , and phase ) is plotted as a line parallel to the time axis, having ordinate value of  and colored according to the value of  − .This representation is very intuitive and readily reveals the main features of the signal.The damping of each EDS component can be estimated from its visible duration, while the color indicates the instantaneous amplitude (darker color means higher amplitude).Each component is associated with a particular frequency and experiences exponential decay with time.Components with high starting amplitude and low damping are dominant and determine the overall character of the signal.Weak components or those with very high damping are less important, as they are often artifacts produced by the MPM to fit the noise.It should be mentioned that the graphical representation of the MPM results is not completely equivalent to that of the results obtained by the other methods.The frequencies on the vertical axis correspond to frequency of the EDS components and not to Fourier components.Even though they are represented by the same parameter, the frequency, the set of EDS represent a different basis for the decomposition of the signals.Damping of a harmonic component leads to widening of its spectral peak.This may result in overlapping of various components in the spectrum.In the EDS representation, the damping influences the apparent length of the horizontal lines.There is no overlapping along the frequency axis in the EDS representation, no matter how close in frequency and how damped the components are.This feature should be considered when one decides which representation is the most suitable for a specific application.
The MPM TFR images have no cross term artifacts, infinite frequency resolution, and provide convenient visualization of damping.However, MPM TFR follows from the EDS-based model, which is only an empirical suggestion and needs to be justified.Indeed, decomposition into a number of EDS is just one possible type of representing the signal.Other basis functions can be selected to fit the same signal, and the resulting TFR may look different than the MPM-based one.One argument in favor of the EDSbased model comes from physical reasons, as described in the Introduction.Another independent proof could be obtained by demonstrating similarity of the time-frequency portraits produced by the MPM TFR and by non-modelbased distributions.Indeed, the time-frequency portrait is an intrinsic property of the signal, and it should not depend on the choice of any particular analysis tool.Various TFDs can visualize this portrait only approximately, while introducing their specific artifacts.If certain TFDs introduce less severe artifacts than others for a given signal type, they are expected to produce more or less similar looking time-frequency portraits (approaching the actual portrait).On the other hand, showing that a model-based distribution produces time-frequency portrait that looks similar to that produced by a non-model-based distribution could be considered as an argument in favor of the model.
The STFT and CWD belong to different classes (Fourier transform-based and bilinear) yet show very similar timefrequency portraits of the same clinical percussion signal.Notably, all dominant features present in the spectrogram and CWD are also present in the MPM TFR plot.Based on the above reasoning, this could be considered an independent proof of the validity of the EDS-based model of clinical percussion signals.On the other hand, MPM TFR reveals additional components, not resolved by conventional distributions.Those components do not play critical role in the signal reconstruction as they have relatively low amplitudes and/or significantly higher damping compared to the dominant ones, determining the signal behavior.The physical origin of these minor components is unclear.The EDS with low amplitude, too high frequency, negative or excessive damping, and so forth could be artificially introduced by MPM to fit the noise or edge discontinuities in the signal.Alternatively, they may represent real (but weak) oscillation modes of the human body and be of practical importance.On the other hand, the noise (which is always present in percussion signals) could be the major cause of the poor resolution of the WVD, CWD, and reassigned spectrogram.Without the noise, these distributions could possibly show much better resolution.In order to answer those questions and to better understand the effect of noise on the behavior of the MPM TFR and other conventional distributions, we constructed two idealized, noise-free signals closely resembling the actual resonant and tympanic percussion waveforms.These simulated signals are analyzed below.

Analysis of Simulated Percussion Signals.
Two simulated noise-free percussion signals were constructed by adding four highest-amplitude EDS from Table 1 and three such EDS from Table 2, respectively, to model resonant and tympanic signals.The resulting signals (1500 points sampled at 48000 Hz) and their spectra are shown in Figure 3.
The simulated "resonant" signal contains three EDS with frequencies of 75.9, 142.5, 197.4,and 299.1 Hz.The first three EDS have well-resolved spectral peaks, while the fourth one is masked by the broad (highly damped) and higher-amplitude third peak.The simulated "tympanic" signal contains three EDS with frequencies of 159.6, 210.9, and 277.0 Hz.All three spectral peaks are very closely spaced, but still resolvable, as the first (and the strongest) two of them have relatively low damping.
The time-frequency analysis of both simulated signals is presented in Figure 4. Figures 4(a) and 4(b) show the results of applying WVD and CWD (with  = 0.8).Visual examination of the WVD plots reveals many additional components, which were not present in the original signal.These are undesirable cross term artifacts produced by the WVD and are making the image difficult to interpret.Still, the WVDs of ideal signals are cleaner than those of real signals (Figure 2).As expected, the level of cross term artifacts in the CWD images is significantly reduced.CWD allows resolving, respectively, three and two dominant components in both images and noticing that their frequencies remain constant while the amplitudes decay with time.The only issue with CWD is its poor frequency resolution in the beginning ( < 6 and  < 12 ms), where all components have high amplitudes and the strongest interference.This happens because CWD suppresses the cross terms at the expense of the resolution.Therefore, even for these simplistic simulated signals neither WVD nor CWD can simultaneously resolve the cross terms problem and increase the resolution.This problem is exacerbated in practical applications, as the real percussion signals always contain noise.Figure 4(c) shows the results of applying a reassigned spectrogram (STFT) with a 1400-point apodization window.This method shows resolution similar to CWD, only with slightly reduced lowresolution zone ( < 5-9 ms).2, but with frequency values of the second and third EDS changed to 310 and 580 Hz.The resulting CWD and reassigned spectrogram (not shown here) successfully resolved these frequencies starting from 3 ms.Further spreading of the EDS frequencies or reducing the damping (less spectral overlapping) leads to further shrinkage of the blurred zone.Therefore, the blurred areas are caused by the spectral overlapping and not by the noise.The observed "blurring" is then a real feature of the noiseless signal.It is due to the attempt to decompose a signal made from EDS into harmonic modes (analogous to decomposition of a vector into a nonorthogonal basis).For these cases, the MPM TFR correctly identified all constituent EDS including weak and highly damped modes, not resolvable by conventional distributions, and eliminates the overlapping along the frequency axis.

Using MPM TFR to Detect Changes in Physical Condition.
Up to this point, we have demonstrated the superiority of the MPM-based TFR over WVD, CWD, and spectrogram when applied to clinical percussion signals.Below we discuss some of its practical applications.The rates of decay (damping), the frequencies, and the relative amplitudes of dominant EDS are important parameters, possibly reflecting pathophysiological properties of the percussed organ.
An illustrative example of how the MPM-based timefrequency analysis can be practically deployed to detect changes in the signal is shown in Figure 5.In this experiment, a normal volunteer did not eat for 15 hours and then consumed a carbonated drink.The percussion signals were recorded over his abdomen and chest before and after consuming the drink.The MPM TFRs of these signals were constructed and compared.It can be seen that before consuming the drink, the least damped EDS at ∼240 Hz responsible for the overall tympanic character of the abdominal signal lasted for only ∼20 ms (see Figure 5, top row).After taking the drink, the intestinal gas pattern changed, so that the frequency of the least damped EDS shifted to ∼190 Hz, and its duration increased beyond 30 ms.Changes in the modal composition of the chest signal can be observed as well (see Figure 5, bottom row).For example, an additional strong but highly damped mode has appeared at ∼300 Hz, as marked by the dashed oval.Other signal modes have remained almost unchanged.One possible explanation of this result is that gas-filled stomach has freed one of the chest vibration modes indirectly through diaphragm.
It should be noted that none of the observed effects can be reliably detected using conventional TFRs examined in this paper due to their inferior resolution in the 0-5 or 0-12 ms interval and due to their insensitivity to slight changes in the damping.

Conclusions
The proposed MPM TFR yields a clear and efficient graphical representation of both resonant and tympanic percussion signal types.The ability to resolve oscillation modes with closely spaced frequencies and high damping makes MPM TFR a preferable method for the analysis of clinical percussion signals.Its superiority becomes especially evident when analyzing the 5 ms time interval following the signal onset, where such commonly used TFDs as STFT, WVD, and CWD all fail to resolve individual frequency components.
We have successfully demonstrated the usefulness of MPM TFR in detecting subtle percussion signal changes caused by physiological processes.We were able to detect the appearance of a new oscillation mode or slight changes in

Figure 1 :
Figure 1: Two examples of audible percussion signals collected from healthy volunteers: (a) Signal from the left subclavicular area ("resonant" character); (b) signal from the left abdominal area ("tympanic" character).The normalized spectra of both signals are shown in (c).

Figure 3 :
Figure 3: Noise-free "resonant" and "tympanic" percussion signals synthesized with four and three highest-amplitude EDS extracted from the real signals shown in Figure 1.The spectra of both signals are shown in (b).

Figure 4 (
d) contains timefrequency representations of the original EDS components used to construct both synthetic signals (MPM TFR).Both images have no cross term artifacts, infinite frequency resolution, and provide convenient visualization of damping.They adequately represent the noise-free resonant and tympanic signals and are clearly superior to WVD, CWD, and spectrogram.By comparing the noise-free plots of Figure4with their real-world counterparts of Figure2it can be observed that the extent of the low-resolution zone in CWD and reassigned spectrogram is approximately the same for the real and noisefree percussion signals.Therefore, the low-resolution zone in the CWD and STFT cannot be explained by the noise in the signal.In another simulated experiment, a signal was synthesized from the same top three EDS of Table

Figure 5 :
Figure 5: MPM TFR of percussion signals from a normal volunteer before and after intake of a carbonated drink.The percussion signals were recorded both over the abdomen and the chest.

Table 1 :
Parameters of all EDS extracted by the MPM algorithm from the resonant percussion signal of Figure1(a).

Table 2 :
Parameters of all EDS extracted by the MPM algorithm from the tympanic percussion signal of Figure 1(b).