2.1. Acquisition of Respiratory Sounds

Sound signals belonging to 15 individuals who had severe asthma, moderate asthma, mild asthma, chronic obstructive pulmonary disease (COPD), hypersensitivity pneumonitis (HP), and the normal subjects were recorded from the right and left part of chest separately. To prevent sounds coming from ambient, sound recording process was carried out in the laboratories of the College of Medicine at University of Gaziantep. Sony ECM T150 microphone with an air capsule was used to record sounds [8]. The working method of the microphones can be one or all directions. The microphone used in this study is one direction. It was mounted onto the stethoscope head to get standard and valuable sound [8]. Its impedance is 2.2 kΩ. Response bandwidth is between 30-15000 Hz. The sampling rate is 8 kHz and 16 bits. The duration of the recorded sounds varies between 10 and 14 second.

Recording respiratory sounds are not completely distinguishable due to combination of noise sources such as heartbeat, muscle movements, and other sounds not associated with respiration. According to literature studies, heart sounds, muscle sounds and frictional noise constitute low frequency range below 100Hz. Furthermore, sounds which have frequency components above 2000Hz are not carry important information about the respiratory sounds. In accordance with this information, we can say that respiratory sounds have considerable frequency components between 100Hz and 2000Hz frequency bands. Therefore, irrelevant frequency components should be filtered.

In the pre-processing stage of this study, the respiratory sound signals were filtered with high pass filter at 100 Hz (to eliminate muscle sounds, heart sounds and frictional noise), and low pass filter at 2000 Hz (to avoid aliasing).

When we examine the literature studies about the respiratory sounds, we have seen that Bessel high pass filter and Butterworth low pass filter are often used to filter respiratory sounds. For example, Pasterkamp and coworkers (1996) [22] utilized a 6th order Butterworth low pass filter with cut off frequency of 2400Hz. Broersen and de Waele (2000) [23] used a 4th order Bessel band pass filter with cut off frequencies of 100Hz and 1500Hz. Güler et al. (2005) [21] filtered out heart sounds and frictional sounds using a sixth-order Bessel Filter. Similarly eight-order Butterworth filter with flat passband characteristics was used as anti-aliasing filter. Sen and Kahya (2003) [24] used Bessel high pass filter (6th order) and Butterworth low pass filter (8th order). Yeginer et al. (2009) [25] also used sixth order Bessel high-pass filters and eighth order Butterworth low-pass filters in order to minimize heart sound interference and frictional noise with minimal phase distortion and to prevent aliasing. Due to the widespread use of these filters, 6th order Bessel high pass filter and 8th order Butterworth low pass filter were used in this study.

After the filtering process, pre-processing stage was continued with segmentation of respiratory sounds. Since the each of sound samples contains a varying number of respiratory cycles, each sound signal was segmented into a small duration of one complete cycle of respiration. And then every cycle was separated as inhalation and exhalation sound signals according to the start and end times of inhalation and exhalation phases. The start and end times of phases as milliseconds were determined by utilizing audacity tool. Audacity is multi-track audio editor and recorder. This tool can be accessed free on the web. After the start and end times of phases were determined, respiratory sound cycles were divided into inhalation and exhalation using Matlab software package. In this way, each segment contained one inhalation or exhalation phase.

As a result of the segmentation of the sounds belonging to 15 people with different pulmonary conditions, sound segments which of amounts are mentioned in the following list were obtained.

• ✓

  67 right basal inhalation sound segments

• ✓

  67 right basal exhalation sound segments

• ✓

  70 left basal inhalation sound segments

• ✓

  70 left basal exhalation sound segments

Every sound segment was evaluated and processed as a separate pattern.

Distributions of the respiratory sound segments according to 6 pulmonary conditions are as shown Fig. 1.

Fig. 1.

Distribution of the lung sound segments according to 6 pulmonary conditions.

Following the pre-processing stage, frequency spectra of sound segments were obtained by using spectral methods such as FFT, AR, ARMA.

2.2. Used Feature Extraction Methods

Feature extraction is a special form of dimensionality reduction [8]. Along with the feature extraction process, high dimensional vectors are converted into lower dimensional vectors. It is a mapping f: RH→RL, where L<H. In the feature extraction process, characteristics of the signal with a small data rate are captured [4].

Frequency components, their distribution and power information at various frequencies are numerous characteristic features for sound signals. To represent the frequency characteristics of sound signals, generally power spectra of signals are estimated [26]. For a given signal, the power spectra give the distribution of the signal power among various frequencies.

The most common way of generating power spectra is to use spectral analyses techniques. For this reason, we used different spectrum analysis methods namely AR, ARMA and FFT in this study. Detailed information about these techniques is given in this section.

Welch method is a FFT based technique used to estimate the power of signal at various frequencies. Spectral estimation through Welch based on the definition of periodogram method [12]. Calculation of PSD with this method is carried out averaging modified periodograms from overlapped and windowed segments [27].

In this method, a signal s(n), of length N samples, is divided into K data segments of length M, overlapping by D points. The ith segment s_i is defined as following Eq. (1) [27].

If the finite-length signal such as signal with N samples is considered, the spectrum of this signal generally exhibits side lobes and spectral leakage due to discontinuities at the end points [27]. Well, what are the leakage and discontinuities? More importantly, how they are prevented?

The FFT-Welch method assumes that the signal is infinite. That is, signal is exactly repeated throughout all time [28]. As a result of this assumption, signal is periodic at intervals that correspond to the length of the time record [28]. For example, we suppose that our input is a single cycle of a sine wave. FFT-Welch method will actually see it as an endless repetition of that sine wave. The resulting spectrum contains just a single spike or main lobe at the signal’s frequency. Fig. 2 illustrates this situation.

Fig. 2.

Single cycle of sine wave, infinite cycles of a sine wave and its spectrum

But we assume that our signal is a partial cycle of a sine wave. In this case, signal is not periodic. The FFT-Welch method will see this as a signal that jumps straight from the maximum value to zero, repeatedly. In this case, discontinuities occur in the signal and spectrum of this signal is different than single spike. In the spectrum of it, the peak gets smeared out across the frequency domain and noise floor rises up. Even, noise may obscure other frequencies that are present in the signal. Fig. 3 shows partial cycle of sine wave, sine wave with discontinuity and its spectrum.

Fig. 3.

Partial cycle of sine wave, sine wave with discontinuity and its spectrum

For the coherent nice spectrum as Fig. 2, measured signal is should be periodic. Also this signal should have integer number of periods [28]. However, generally, the measured signal isn’t an integer number of periods. For the FFT-Welch, it means that this signal is discontinuous. In the spectrum of discontinuous signal, spectral leakage and side lobes are seen.

Spectral leakage distorts the measurement in such a way that energy from a given frequency component is spread over adjacent frequency lines or bins [28]. It causes of the fine spectral lines to spread into wider signals.

The effects of performing an FFT over a non-integer number of cycles can be minimized by using a technique called windowing [28]. The discontinuities are alleviated and the spectral leakage are minimized by the window function w(n) [27]. A window function alters signal, tapering it to nearly zero at the beginning and end. Therefore, in the welch method; each segment is windowed prior to computing the modified periodogram [27]. That is, we select data window w(n), n=0,…,M-1, and form the sequences s₁(n)w(n),…,s_K-1(n)w(n).

After the windowing operation, for the ith windowed segment, the modified periodogram is calculated using Eq. (2) [26].

In Eq. (2), P(w, i) is the modified periodogram of ith windowed segment, w(n) is the window function and U is the power in the window function. U is defined as Eq. (3) [27].

As shown in Eq. (4) [27], the Welch power spectrum is the average of the K-modified P(w, i) periodograms that were obtained from the overlapped and windowed segments of a signal [27].

Power spectrum estimation using the FFT method is relatively simple, well understood and easy to compute. This method is also fast [29].

The parametric (model-based) methods are also commonly used for spectral estimation. AR and ARMA are the best known of these for spectral analysis of breath sounds. They provide good frequency resolution. Also, they avoid the problems caused by FFT method [29]. However, selection of order is very important in these methods [29]. If a model with too low order is chosen, highly smoothed spectrum is obtained. If order of model is selected too high, we faced with the risk of spurious low-level peaks in the spectrum [29].

Basic idea of the parametric spectral estimation is shown Fig. 4.

Fig. 4.

Block diagram of parametric spectral estimation

The AR model is an infinite impulse response filter or an all pole filter [29] and is named according to the number of observations in past periods. AR (p) model can be defined by AR parameters {a (1), a (2)… a (p)} and white noise pw. There are several algorithms to estimate parameters of AR such as burg, yule walker.

PSD of AR methods is obtained basically from following Eq. (5) and Eq. (6).

In Eq. (5), pw refers to white noise variance to scale the returned PSD. In the same equation, T means that interval in seconds to scale the returned PSD.

In Eq. (6), A(p) is referred AR model of degree p.

The power spectrum estimation is carried out by the Burg algorithm as in Eq. (7).

Where êp represents the total least square error and it is the sum of the forward and backward prediction errors.

ARMA model is a combination of both AR and moving average (MA) model. ARMA(p,q) method can be defined ARMA parameters { a(1), a(2), …, a(p),b(1),b(2), …,b(q)} and white noise pw. PSD of this model is obtained by filtering unit variance noise with a filter having both poles and zeros [29]. PSD of ARMA model is expressed by Eq. (8) and Eq. (9). B(q) is the MA model of degree q.

AR model is suitable for signals that have sharp peaks. MA model is suitable for signals that have not sharp peaks. The ARMA model is suitable for both types of signals. Sounds have sharp peaks. Therefore, use of AR and ARMA models are more appropriate than MA model for sound signal analysis. AR model is more advantageous than ARMA model in terms of mathematical computing density [29].

One of the most important aspects of the use of AR and ARMA methods is the selection of the model orders. Since inhalation and exhalation phases show different characteristics, separate AR-Burg models were applied to each phase in this study. Changing order values were tried for AR and ARMA models and the order values that produce the highest classification accuracy were chosen. While the order of the AR-Burg model was set as 9 for inhalation phases of right basal, it was set as 14 for exhalation phases of same basal. Orders of AR-Burg which are applied to left basal inhalation and exhalation phases were determined respectively as 22 and 16. Orders of ARMA model used in this study are shown following expressions.

p, q order values of ARMA model;

• ✓

  Inhalation sounds of right basal: 21, 20

• ✓

  Exhalation sounds of right basal: 14, 13

• ✓

  Inhalation sounds of left basal: 25, 24

• ✓

  Exhalation sounds of left basal: 27, 26

Power spectrum densities of partitioned sound signals were found through AR-Burg’s method, ARMA method and FFT-Welch method in a Matlab software package.

nfft parameter is crucial for FFT-Welch, AR-Burg and ARMA matlab functions because its value determines the size of the feature vectors. Well, what is the nfft parameter? This parameter refers to the number of FFT points used to calculate the PSD estimate. That is, nfft is the length of FFT and it determines size of PSD. Assuming that the value of length of FFT (nfft) is L, there are (L/2)+1 points in the PSD.

So, how this parameter determines size of feature vectors? In this study PSDs were considered feature vectors of sounds. If we decide that how many points are there in PSD, we determine size of our feature vectors. According to results of experimental studies, we decided that PSDs with 513 points should be obtained as feature vectors in this study. For this reason, value of nfft parameter was determined as 1024.

2.3. Feature Extraction of Respiratory Sounds

In this paper, AR-Burg, ARMA, FFT-Welch methods were applied to each of the inhalation and exhalation sound segments for power spectral density (PSD) estimation. Vectors which consist of value of these PSD were obtained as features of sounds. Fig. 5 shows PSDs of normal healthy person’s inhalation and exhalation sound segments.

Fig. 5.

PSDs of inhalation and exhalations sound segments belonging to healthy person

Fig. 6 indicates PSDs of asthmatic patient’s inhalation and exhalation sound segments.

Fig. 6.

PSDs of inhalation and exhalation sound segments belonging to asthmatic patient

PSDs of patient who had COPD disease is seen as in Fig. 7.

Fig. 7.

PSDs of inhalation and exhalation sound segments belonging to COPD patient

PSDs of inhalation and exhalation sounds relate to patient with HP are shown in Fig. 8. PSDs in all figures are calculated Ar-Burg method.

Fig. 8.

PSDs of inhalation and exhalation sound segments belonging to HP patient

As it is seen from the Fig. 5 and Fig. 6, PSDs of normal and asthmatic respiratory sounds are accumulated between 100 and 600 frequency range. It is understood from these figures that PSDs of exhalation sounds belonging to healthy person and asthmatic patient show more diversity than PSDs of inhalation sounds. Fig.7 and Fig.8 show that PSDs of respiratory sounds related to COPD and HP patients are scattered up to 1000 Hz.

PSDs of inhalation sounds of COPD and HP patients show differences according to that of healthy person. In addition, PSD of exhalation sounds of HP patient are different from that of healthy person.

Because the nfft values were chosen as 1024 for all methods, PSDs with 1024/2+1 (513) points were obtained as feature vectors. The dimension of feature vectors is too big for an effective classification. For the successful and short time classification process, dimension of features vectors should be reduced. In order to reduce the dimensionality of the extracted feature vectors, we divided 0-2000Hz frequency band of PSDs into intervals as presented in Table 1.

Division of the frequency bands into intervals was carried out using decomposition rule of discrete wavelet transform (DWT). In discrete wavelet transform, signal is passed through both the high and low pass filters simultaneously. Cut-off frequencies of these filters are designed as one-fourth of the sampling frequency. For the 8000 Hz sampling frequency, DWT produce frequency intervals like Fig. 9. H represents high pass filter and L represents low pass filter in Fig. 9.

Fig. 9.

Ranges of frequency bands in discrete wavelet decomposition

In many studies, wavelet coefficients of frequency intervals produced by DWT have been used as features of respiratory sounds. But, in this study, decomposition rule of DWT was only used to determine the useful frequency intervals.

For our respiratory sounds with 8000 Hz sampling frequency, computed PSDs are scattered 0Hz-2000Hz frequency bands. Therefore, for each PSD, our frequency intervals are F1, F2, F3, F4, F5, F6 and F7. These intervals were generated using decomposition rule of DWT. Since the respiratory sounds do not have any useful frequency and power components in low frequencies, we discarded F6 and F7 frequency intervals. First five frequency intervals were selected to be used to create feature vectors.

Statistical features were extracted from the PSD values belonging to first five frequency intervals. These statistical features;

1. 1.

   Mean of the absolute values of the PSD values in each selected frequency interval.

2. 2.

   Maximum of the PSD values in each selected frequency interval.

3. 3.

   Standard deviation of the PSD values in each selected frequency interval.

4. 4.

   Ratio of the absolute mean values of adjacent frequency intervals.

Feature 1 (mean) and feature 2 (maximum) features represent the frequency distribution of the sound signal. Feature 3 (standard deviation) and feature 4 (ratio of the absolute mean values) describe the amount of changes in frequency distribution [1].

Kandaswamy et al. (2004) [1] extracted these features from the wavelet coefficients and obtained high performance. So, we used the same statistical features for classification.

Feature vectors of the sound segments were created using these statistical features as seen in Table 2. 5 mean values, 5 max values, 5 standard deviation values and 4 ratio of the absolute mean values are the features of our sound segments.

After the feature extraction phase, feature vectors consist of these 19 values applied as input parameters into ANN.

2.4. Classification Using Artificial Neural Network

ANN is an interconnected assembly of simple parallel processing elements, units, nodes or neurons, whose functionality is loosely based on animal neuron [30]. ANN has parallel processing in these neurons and connections with weights between neurons. ANN has to be trained for updating the weights and biases in order to produce the desired mapping. ANN consists of two phases as training and testing. At the training phases, the feature vectors are applied as input to the input layer of network and the network updating the weights and biases, to capture the relationship between the input patterns and desired outputs. The combined inputs are then modified by a transfer functions. This transfer functions can be a threshold function with only passes information if the combined input reaches a certain level, or it can be a continuous function of the combined input. This training phase often takes many hours. The ANN testing phase is real-time which refers to how the network globally processes at its inputs buffer and creates a respond at the output buffer. So, ANN is mostly used to solve complex problems of pattern classification and optimization tasks [1, 31]. ANN is widely used in biomedical field for modelling, data analysis, and diagnostic classification [1, 32–34]. The most frequently used training algorithm in classification problems is the back-propagation (BP) algorithm [1].

In this paper, ANN was used for classification of normal and different pulmonary disorders. İnput layer contains one layer and the numbers of node in this layer are 19. The hidden layer had one layer and 15 nodes for inhalation sounds, 18 nodes for exhalation sounds which gave maximum classification accuracy, the output layer consisted of one layer and 6 nodes for the classification. Used ANN structer is shown in Fig. 10.

Fig. 10.

Used ANN structure

Orange canvas toolbox [35] was used for ANN. This toolbox is available for free at www.ailab.si/orange. Version of used toolbox is 2.7.

Q1, Q2, Q3, Q4, Q5 and Q6 outputs represent mild asthma, moderate asthma, severe asthma, normal, COPD and HP.