Real-Time Implementation of a Frequency Shifter for Enhancement of Heart Sounds Perception on VLIW DSP Platform

Abstract

Auscultation of heart sounds is important to perform cardiovascular assessment. External noises may limit heart sound perception. In addition, heart sound bandwidth is concentrated at very low frequencies, where the human ear has poor sensitivity. Therefore, the acoustic perception of the operator can be significantly improved by shifting the heart sound spectrum toward higher frequencies. This study proposes a real-time frequency shifter based on the Hilbert transform. Key system components are the Hilbert transformer implemented as a Finite Impulse Response (FIR) filter, and a Direct Digital Frequency Synthesizer (DDFS), which allows agile modification of the frequency shift. The frequency shifter has been implemented on a VLIW Digital Signal Processor (DSP) by devising a novel piecewise quadratic approximation technique for efficient DDFS implementation. The performance has been compared with other DDFS implementations both considering piecewise linear technique and sine/cosine standard library functions of the DSP. Piecewise techniques allow a more than 50% reduction in execution time compared to the DSP library. Piecewise quadratic technique also allows a more than 50% reduction in total required memory size in comparison to the piecewise linear. The theoretical analysis of the dynamic power dissipation exhibits a more than 20% reduction using piecewise techniques with respect to the DSP library. The real-time operation has been also verified on the DSK6713 rapid prototyping board by Texas Instruments C6713 DSP. Audiologic tests have also been performed to assess the actual improvement of heart sound perception. To this aim, heart sound recordings were corrupted by additive white Gaussian noise, crowded street noise, and helicopter noise, with different signal-to-noise ratios. All recordings were collected from public databases. Statistical analyses of the audiological test results confirm that the proposed approach provides a clear improvement in heartbeat perception in noisy environments.

1. Introduction

Cardiac auscultation is still a very important technique for primary cardiovascular assessment of patients. Heart sounds are generated by the flow of blood through the heart chambers as the heart valves open and close during the heart cycle [1,2,3,4]. In healthy subjects, two normal sounds can be perceived for each heartbeat, and are commonly, onomatopoeically described as ‘lub’ and ‘dub’.

Normal body sounds have a very low sound level in the range of 22–30 dB in free space. A simple stethoscope is able to render, for the operator’s ear, a significantly higher sound level of 60–65 dB, while also attenuating the contribution of external noises (thus highly increasing the signal-to-noise ratio). Despite the simplicity of the instrument, efficiency in mechanical transduction is very important in order to render the lower frequencies, at the edge of audibility, of which heart sounds are composed. The spectral content of the main heart sounds (S1 and S2) is confined to low frequencies, spanning up to about 100 Hz. The human ear threshold of audibility in this frequency range is particularly high, then the perceived heart sounds are particularly feeble. Figure 1 schematically shows the spectrum of heart sounds in relation to the human ear threshold of audibility. Heart sounds are too weak to be heard because the human ear is far less sensitive to them with respect to the speech range (1000–2000 Hz).

For this reason, the suppression of external noises from a stethoscope is an important quality parameter along with audio amplification. With the advent of digital stethoscopes, it has been possible to filter and amplify heart sounds to improve their perception [7,8,9,10]. Today it is also possible to use special processing [11] to improve diagnostics, as well as to use devices to visualize signals [12]. Moreover, novel promising wearable devices have been proposed, which could enable continuous, long-term monitoring of heart sounds [13,14,15,16]. Cardiac auscultation is an invaluable tool to assess the health status of a patient in emergency medicine contexts, specifically in extramural settings, where other cardiac monitoring devices, such as those for electrocardiography and photoplethysmography, or for wireless monitoring of pacemakers [17], are more difficult or impossible to use.

1.1. Methods for Noise Suppression

Many researchers have proposed a variety of methods [18] for heart sound denoising, either based on hardware or software solutions. Among software solutions, a first class of algorithms is based on thresholding applied in transform domains, including Short Time Fourier Transform [19,20,21,22], Fast Fourier Transform [23,24,25], Wavelet transform [19,26,27,28,29,30,31,32,33,34,35], and Convolution Wavelet Packet Transform [36,37,38,39]. However, the choice of the threshold function [35] strongly affects the denoising effect. In addition, with STFT it is impossible to extract some useful clinical information, such as the time distance between aortic valve closure and pulmonic valve closure [19]. A second class of algorithms is based on mode decomposition methods, such as Empirical Mode Decomposition [40,41,42].

Ensemble Empirical Mode Decomposition [43,44], Complementary Ensemble Empirical Mode Decomposition [45,46], and Variational Mode Decomposition [47,48,49] are mode decomposition methods that decompose the heart sound signal into a series of intrinsic mode functions (IMFs) and specific algorithms must be designed to select and retain only the IMFs that capture the main information content of the signal, during the reconstruction of the signal, while discarding those related to noise. However, this class of algorithms, particularly the EMD, is affected by the mode mixing issue, which is caused by the intermittency of signals and noises [42]. In addition, mode decomposition algorithms usually require the a priori selection of the number of IMFs, which affects the decomposition results and determines the impact of mode mixing. EEMD was first conceived to mitigate the mode mixing observed in EMD; however, the reconstructed signal is usually affected by residual noise [49]. CEEMD was proposed precisely to eliminate such residual noise. VMD performs better than EMD and its variants, by realizing an adaptive decomposition of each mode to overcome the mode mixing problem [47]. It has been used to separate heart sounds and lung sounds and to improve the detection of the subtle third heart sound (S3) [48]. However, an accurate selection of the VMD parameters is instrumental in achieving a good decomposition effect [49]. Additional algorithms for heart sound denoizing are Total Variation filtering [50,51] and Blind Source Separation [52,53,54].

Hardware solutions are generally based on the use of additional microphones or PCG sensors, which are used to obtain independent measurements of external noises, and allow enhancing the quality of noisy heart sound recordings via subtractive techniques. This is the case of the spectral subtraction algorithm described in [55], whose performance, however, turned out to be heavily dependent on the characteristics of the additional microphone used to pick up the external noises. A different approach was described in [56], where an array of microphones, appropriately positioned in specific measurement points on the body, and a frequency-domain Independent Component Analysis are used to implement Direction of Arrival and Beamforming techniques for noise suppression in heart sound recordings. The main limitation of such hardware solutions is precisely the need for additional sensors that have to be placed on a subject’s body, possibly increasing the subject’s discomfort (especially in the case of sensor arrays), and also requiring accurate positioning to ensure effective results.

It is important to highlight that all the methods described above have not been designed for real-time operation, nor to enhance human audibility, but rather to improve computer analysis. For this reason, the ability of these methods to improve human audibility of heart sounds has not been tested on cohorts of human subjects.

1.2. Pitch Shift Approaches

Another possibility for improving the perception of heart sounds can be achieved by shifting the heart sounds to higher frequencies. In fact, as the frequency increases, the threshold of audibility lowers, and sounds are perceived even better. Methods based on this approach have proven effective in situations characterized by considerable external noise, such as aboard helicopters [57,58], aircraft [59,60], ambulances [61] or other evacuation vehicles, emergency rooms [62,63], rural clinics [64], and crowded places such as stadiums or crowded streets. Frequency translation can be achieved through digital signal processing techniques. Signal modulation techniques, such as those used for radio or telephone transmissions, can be used for this purpose. Some techniques for translating sounds to the best frequencies perceived by the human ear realize a conspicuous frequency translation of the heart sounds, either to the thousands of Hertz [65] or with a translation of 500 Hz [66], achieving considerable improvements in the acoustic perception of the operator. However, if the pitch is increased significantly, the sounds are perceived as artifactual or unreal. Another study [67] has proposed an adjustable frequency modulator with a significantly reduced frequency shift, realized using Hilbert transform. In this case, the Hilbert transformer was implemented as a simple time delay, based on the assumption that heart sounds have a narrow bandwidth. In addition to stethoscopes with frequency-shifting capabilities, hybrid stethoscopes have been proposed, which embed ultrasonic sensors that make the Doppler signal audible, so as to enable cardiac diagnostics in scenarios of intense external noise [68,69]. The enhancements of human audibility of heart sounds provided by the methods proposed in [66,67] have not been tested on an actual cohort of human subjects. The method proposed in [65] appears to have been tested on a cohort of human subjects, but the kinds and levels of noise adopted have not been specified. The methods proposed in [68,69] require very peculiar stethoscopes with ultrasonic sensors, which are not very common.

In summary, there is a gap in the literature regarding an accurate, real-time method for frequency shifting of heart sounds, where the ability to improve the human audibility of heart sounds has been evaluated by experimental testing on a real cohort of human subjects.

The present study aimed to fill this gap and describes the design, implementation, and experimental testing of a novel, real-time, frequency shifter for heart sounds, implemented on a Very-Long-Instruction-Word (VLIW) Digital Signal Processor (DSP). The proposed method is based on the Hilbert transform and provides two improvements with respect to the solution of [67]: (i) the Hilbert transformer is implemented using an FIR filter to improve signal quality; (ii) a Direct Digital Frequency Synthesizer (DDFS) is used to allow the agile modification of the frequency shift. The proposed implementation allows translating the spectrum of heart sounds even by only a few tens of Hertz. A novel piecewise quadratic technique is also proposed for the implementation of the DDFS and demonstrates the best trade-off between execution time, memory resources, and power consumption. Auscultation tests were performed on 41 volunteers to assess the actual performance of the proposed method. The results confirm the capability of the proposed method to enhance the audibility of heart sounds in the presence of different kinds of noise. The frequency shifter here described could be easily embedded in modern electronic stethoscopes at very modest costs, thus providing valuable support in in emergency medicine contexts, such as those described in [58,59,60,61,62,63,64,65], where heart sounds auscultation is fundamental to assess the health status of a patient and could be seriously hindered by environmental noises.

2. Materials and Methods

2.1. Top-Level Algorithm

The top-level architecture of the algorithm used to implement the frequency shift operation is shown in Figure 2. The operation is achieved using a single-sideband modulator realized using the Hilbert transform. The behavior of this solution can be briefly explained as follows. The Hilbert transform, which can be considered to be a linear system with a frequency response:

$$H(f)=\{\begin{array}{cc}-j\hfill & \hfill f>0\\ j\hfill & \hfill f<0\end{array}$$

(1)

is used to form the analytical signal x_c(n) associated with x(n) as follows:

$$x_c(n)=x(n)+j x_H(n)$$

(2)

x_H(n) being the Hilbert transform of x(n). As is well-known, the Fourier transform of the associated analytical signal includes only the positive-frequency axis components of X(n), as follows:

$$X_c(f)=\{\begin{array}{cc}2X(f)\hfill & \hfill f>0\\ 0\hfill & \hfill f<0\end{array}$$

(3)

The frequency shift operation, therefore, can be computed by multiplying the complex signal x_c(t) by $e^{j\hspace{0.17em}2\pi \hspace{0.17em}{f}_{shift}/f_s\hspace{0.17em}n}$, f_shift being the desired frequency shift and f_s being the sampling frequency, and taking the real part, in order to re-obtain a real signal with Hermitian Fourier Transform, obtaining the following:

$$Y(f)=\{\begin{array}{cc}X(f-f_{shift})\hfill & f>f_{shift}\hfill \\ 0\hfill & 0\le f<f_{shift}\hfill \end{array}$$

(4)

where clearly:

$$y(n)=\mathrm{Re}\left[x_c(n)\cdot e^{j\hspace{0.17em}2\pi \hspace{0.17em}{f}_{shift}/f_s\hspace{0.17em}n}\right]$$

(5)

The architecture of Figure 2 implements (5) using a Direct Digital Frequency Synthesizer (DDFS) to compute the two signals cos(2π f_shift/f_s n) and sin(2π f_shift/f_s n).

Please note that a similar technique has been devised in [67], where, however, oscillators are used to obtain sine/cosine signals and the Hilbert transform is approximated using a simple delay, chosen in order to realize the required π/2 phase shift (see (1)) only at a given frequency, f_mid, chosen in the middle of the input signal bandwidth, i.e., by calling δ the implemented delay, 2π∙f_mid∙δ = π/2. This unfortunately leads to obtaining a not negligible output noise power at frequencies lower than f_shift (low-sideband noise power) since in the architecture of Figure 2, the final adder which produces y(n) actually implements a cancellation between the two adder inputs for f < f_shift and, when the Hilbert transform is implemented with a simple delay, the cancellation is perfectly realized only at frequency f_shift − f_mid. As an example, Figure 3a shows the spectrum of a cardiac audio signal sampled with a frequency of 2 kHz and used as an input signal for the system of Figure 2, and Figure 3b shows the corresponding elaborated signal when the Hilbert transform is implemented with a simple delay, using f_mid = 50 Hz (this value has been chosen in order to minimize the low-sideband noise power). As it can be observed, the output spectrum nullifies only at frequency f_shift − f_mid = 50 Hz, and no negligible frequency components can be seen in the low-sideband. As reported in

Table 1. Ratio between low-sideband noise power and output power by considering different techniques to implement the Hilbert transform.

Technique	Low-Sideband Noise Power/Output Power
Simple delay approx. [67]	−11.9 dB
Proposed with 20-th order FIR filter	−21.2 dB
Proposed with 40-th order FIR filter	−23.3 dB
Proposed with 60-th order FIR filter	−30.3 dB
Proposed with 80-th order FIR filter	−37.7 dB
Proposed with 100-th order FIR filter	−42.9 dB

, the overall ratio between the low-sideband noise power and the output power is −11.9 dB.

To reduce the low-sideband noise, in this paper, we implemented the Hilbert transform using a linear-phase FIR filter, whose coefficients have been obtained using the Parks–McClellan optimal FIR filter design algorithm implemented in Matlab. Please note that in this case, the FIR filter implements the Hilbert transform with a delay given by half of the filter order. To compensate for this delay, a delay Δ also equal to one-half the filter order is added in Figure 2 at the input of the multiplier-driven by x(t). Figure 3c shows the obtained output spectrum in the case of a 40th-order FIR filter. Table 1 summarizes the achieved ratio between low-sideband noise power and output power in the case of FIR filter orders 20, 40, 60, 80, and 100. It can be observed that a significant reduction of the low-sideband noise power can be achieved with respect to the previous approach of [67]. In the following sections, we will investigate the implementation with a 40th-order filter.

2.2. Implementation

In this paper, we will focus on an implementation using a Texas Instruments VLIW DSP. In the following sections, we will describe in detail the implementation of the two main components of the algorithm of Figure 2 that is the Hilbert transform, implemented using a FIR filter, and the DDFS.

2.2.1. FIR Filter Implementation (Hilbert Transform)

As well-known DSPs are specifically designed to implement in an efficient way FIR filters. The exploited implementation (shown in Figure 4) uses a simple direct form implementation of the filter, using an array managed as a linear buffer to store the input samples x(n), x(n − 1), …, x(n − M), M being the filter order. Please note that the delay Δ shown in Figure 2 can be easily realized by taking the sample x(n − M/2) from the array.

In the figure, the notation Q_m,n is used to indicate a signed 2’s complement representation with a most-significant bit of weight 2^m and a less-significant bit of weight 2⁻ⁿ. As can be observed, the input is assumed contained in the range [−1, 1) and is represented in Q_0,15 using 16 bits. This size fits the data-size of the inputs of the multiplier available in the DSP. The filter coefficient, h_i, is also assumed represented in Q_0,15, using 16 bits. To manage the overflow, the sum-of-products is computed using the guard-bits approach and the available data-size provided by the DSP is computing the sum using 40 bits in Q_9,30. The filter output is obtained in Q_0,15 (note that the 16-bits representation is also needed to compute the multiplication between x_H(n) and sin(2π f_shift/f_s n) in Figure 2) using a rounding operation on the less-significant bits and a saturation operation to reduce the most-significant bits. The employ of saturation helps in mitigating the limited overflow cases that could happen whether a full-scale, large bandwidth input signal is applied.

2.2.2. DDFS Implementation

The top-level DDFS implementation is shown in Figure 5. A phase accumulator is used to compute instantaneous phase φ/2π encoded in an unsigned N-bits representation (in our implementation N = 32 to fit DSP register size). This phase is then truncated to P-bits. The heart of the DDFS is the two blocks that compute sine and cosine functions in [0, π/4) interval that is the two functions:

$$y_s(z)=\mathrm{Sin}\left[\frac{\pi }{4}\cdot z+{\phi }_{LSB}\right]$$

(6)

$$y_c(z)=\mathrm{Cos}\left[\frac{\pi }{4}\cdot z+{\phi }_{LSB}\right]$$

(7)

for z ∈ [0, 1). The 1’s complement, the swapper, and the two inverters (inv.), driven by the mapper, are used to reconstruct the output sine and cosine functions in the whole [0, 2π) interval. Please note that in (6) and (7), φ_LSB represents the angle corresponding to the one-half LSB of signal φ/2π after P-bits truncation, which is: φ_LSB = 2π∙2^−P/2.

In the literature, many approaches have been proposed to implement (6) and (7), including simple Lookup-Table (LUT) methods, LUT-compression techniques, piecewise approximations, angle rotation approaches (CORDIC algorithm), and mixed LUT and angle rotation techniques. In this paper, we introduce an extension of the piecewise linear technique proposed in [70] to the piecewise quadratic case.

In piecewise approximation techniques the interval [0, 1) of z is subdivided into 2^H sub-intervals, each of length h = 2^−H. Since we consider a piecewise quadratic approximation, in the generic i-th sub-interval, the two functions y_s(z) and y_c(z) are approximated as follows:

$$y_s(z)\simeq a_{0\hspace{0.17em}i}+a_{1\hspace{0.17em}i}\cdot z_i+a_{2\hspace{0.17em}i}\cdot z_i^2+{\epsilon }_{s\hspace{0.17em}2\hspace{0.17em}i}(z_i)$$

(8)

$$y_c(z)\simeq b_{0\hspace{0.17em}i}+b_{1\hspace{0.17em}i}\cdot z_i+b_{2\hspace{0.17em}i}\cdot z_i^2+{\epsilon }_{c\hspace{0.17em}2\hspace{0.17em}i}(z_i)$$

(9)

where z_i represents the distance between z and the starting point of the sub-interval: z_i = z − i · h.

By following the approach of [70], we consider a truncated Chebyshev expansion in order to obtain the coefficients a_ji and b_ji of (8) and (9) which approximate (6) and (7) in the quasi-minimum maximum error sense. To that purpose the variable z_i is substituted by variable w so that the interval [i∙h, (i + 1)∙h) of z_i corresponds to the interval [−1, +1) of w that is:

$$w=\frac{2}{h}\cdot z_i-1$$

(10)

Using this substitution, the functions to be approximated in the i-the sub-interval becomes:

$$y_{s\hspace{0.17em}i}(w)=\mathrm{Sin}\left[\frac{\pi }{4}\cdot \left(\frac{h}{2}\cdot w+\frac{h}{2}+i\cdot h\right)+{\phi }_{LSB}\right]$$

(11)

$$y_{c\hspace{0.17em}i}(w)=\mathrm{Cos}\left[\frac{\pi }{4}\cdot \left(\frac{h}{2}\cdot w+\frac{h}{2}+i\cdot h\right)+{\phi }_{LSB}\right]$$

(12)

These two functions are expanded in a series of Chebyshev polynomials of the first kind (T_n(w)):

$$y_{s\hspace{0.17em}i}(w)\simeq {\alpha }_{0\hspace{0.17em}i}\cdot T_0(w)+{\alpha }_{1\hspace{0.17em}i}\cdot T_1(w)+{\alpha }_{2\hspace{0.17em}i}\cdot T_2(w)+{\alpha }_{3\hspace{0.17em}i}\cdot T_3(w)+\dots$$

(13)

$$y_{c\hspace{0.17em}i}(w)\simeq {\beta }_{0\hspace{0.17em}i}\cdot T_0(w)+{\beta }_{1\hspace{0.17em}i}\cdot T_1(w)+{\beta }_{2\hspace{0.17em}i}\cdot T_2(w)+{\beta }_{3\hspace{0.17em}i}\cdot T_3(w)+\dots$$

(14)

where we recall that Chebyshev polynomials of the first kind T_n(w) can be written using the following recursion:

$$\{\begin{array}{l}{T}_0(w)=1\\ T_1(w)=w\\ T_{n+1}(w)=2w\cdot T_n(w)-T_{n-1}(w)\end{array}$$

(15)

Since Chebyshev polynomials are orthogonal in the general sense in [−1, +1] interval, the coefficient in (13) and (14) can be obtained as follows:

$${\alpha }_{j\hspace{0.17em}i}=\frac{{\displaystyle {\int }_{-1}^{+1}\frac{1}{\sqrt{1-w^2}}\cdot y_{s\hspace{0.17em}i}(w)\cdot T_j(w)\hspace{0.17em}dw}}{{\displaystyle {\int }_{-1}^{+1}\frac{1}{\sqrt{1-w^2}}\cdot T_j^2(w)\hspace{0.17em}dw}}$$

(16)

$${\beta }_{j\hspace{0.17em}i}=\frac{{\displaystyle {\int }_{-1}^{+1}\frac{1}{\sqrt{1-w^2}}\cdot y_{c\hspace{0.17em}i}(w)\cdot T_j(w)\hspace{0.17em}dw}}{{\displaystyle {\int }_{-1}^{+1}\frac{1}{\sqrt{1-w^2}}\cdot T_j^2(w)\hspace{0.17em}dw}}$$

(17)

By solving the integrals, we obtain:

$${\alpha }_{0\hspace{0.17em}i}=J_0\left(\frac{\pi }{8}\cdot h\right)\cdot \mathrm{Sin}\left[\frac{\pi }{4}\cdot \left(i\cdot h+\frac{h}{2}\right)+{\phi }_{LSB}\right]$$

(18)

$${\alpha }_{1\hspace{0.17em}i}=2\cdot J_1\left(\frac{\pi }{8}\cdot h\right)\cdot \mathrm{Cos}\left[\frac{\pi }{4}\cdot \left(i\cdot h+\frac{h}{2}\right)+{\phi }_{LSB}\right]$$

(19)

$${\alpha }_{2\hspace{0.17em}i}=-2\cdot J_2\left(\frac{\pi }{8}\cdot h\right)\cdot \mathrm{Sin}\left[\frac{\pi }{4}\cdot \left(i\cdot h+\frac{h}{2}\right)+{\phi }_{LSB}\right]$$

(20)

$${\alpha }_{3\hspace{0.17em}i}=-2\cdot J_3\left(\frac{\pi }{8}\cdot h\right)\cdot \mathrm{Cos}\left[\frac{\pi }{4}\cdot \left(i\cdot h+\frac{h}{2}\right)+{\phi }_{LSB}\right]$$

(21)

$${\beta }_{0\hspace{0.17em}i}=J_0\left(\frac{\pi }{8}\cdot h\right)\cdot \mathrm{Cos}\left[\frac{\pi }{4}\cdot \left(i\cdot h+\frac{h}{2}\right)+{\phi }_{LSB}\right]$$

(22)

$${\beta }_{1\hspace{0.17em}i}=-2\cdot J_1\left(\frac{\pi }{8}\cdot h\right)\cdot \mathrm{Sin}\left[\frac{\pi }{4}\cdot \left(i\cdot h+\frac{h}{2}\right)+{\phi }_{LSB}\right]$$

(23)

$${\beta }_{2\hspace{0.17em}i}=-2\cdot J_2\left(\frac{\pi }{8}\cdot h\right)\cdot \mathrm{Cos}\left[\frac{\pi }{4}\cdot \left(i\cdot h+\frac{h}{2}\right)+{\phi }_{LSB}\right]$$

(24)

$${\beta }_{3\hspace{0.17em}i}=2\cdot J_3\left(\frac{\pi }{8}\cdot h\right)\cdot \mathrm{Sin}\left[\frac{\pi }{4}\cdot \left(i\cdot h+\frac{h}{2}\right)+{\phi }_{LSB}\right]$$

(25)

where J_n(x) represents the n-th order Bessel function of the first kind.

The quasi-minimum maximum error quadratic approximation of y_s and y_c can be obtained by truncating (13) and (14) to the second order that is:

$$y_s(z)\simeq {{\alpha }_{0\hspace{0.17em}i}\cdot T_0(w)+{\alpha }_{1\hspace{0.17em}i}\cdot T_1(w)+{\alpha }_{2\hspace{0.17em}i}\cdot T_2(w)|}_{w=\frac{2}{h}{z}_i-1}$$

(26)

$$y_c(z)\simeq {{\beta }_{0\hspace{0.17em}i}\cdot T_0(w)+{\beta }_{1\hspace{0.17em}i}\cdot T_1(w)+{\beta }_{2\hspace{0.17em}i}\cdot T_2(w)|}_{w=\frac{2}{h}{z}_i-1}$$

(27)

By comparing (26) and (27) with (8) and (9), and using (15) we obtain:

$$a_{0\hspace{0.17em}i}={\alpha }_{0\hspace{0.17em}i}-{\alpha }_{1\hspace{0.17em}i}+{\alpha }_{2\hspace{0.17em}i}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}_{1\hspace{0.17em}i}=\frac{2}{h}\cdot \left({\alpha }_{1\hspace{0.17em}i}-4\hspace{0.17em}{\alpha }_{2\hspace{0.17em}i}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{a}_{2\hspace{0.17em}i}=\frac{8}{h}\cdot {\alpha }_{2\hspace{0.17em}i}$$

(28)

$$b_{0\hspace{0.17em}i}={\beta }_{0\hspace{0.17em}i}-{\beta }_{1\hspace{0.17em}i}+{\beta }_{2\hspace{0.17em}i}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{b}_{1\hspace{0.17em}i}=\frac{2}{h}\cdot \left({\beta }_{1\hspace{0.17em}i}-4\hspace{0.17em}{\beta }_{2\hspace{0.17em}i}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{b}_{2\hspace{0.17em}i}=\frac{8}{h}\cdot {\beta }_{2\hspace{0.17em}i}$$

(29)

It is worth noting that the first-order approximation of [70] can be easily re-obtained by imposing α_{2 i} = β_{2 i} = 0 in (28) and (29).

The approximation error in (8) and (9) can be approximated as the first term dropped from the series (13) and (14):

$${\epsilon }_{s\hspace{0.17em}2\hspace{0.17em}i}(z_i)\simeq {{\alpha }_{3\hspace{0.17em}i}\cdot T_3(w)|}_{w=\frac{2}{h}{z}_i-1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_{c\hspace{0.17em}2\hspace{0.17em}i}(z_i)\simeq {{\beta }_{3\hspace{0.17em}i}\cdot T_3(w)|}_{w=\frac{2}{h}{z}_i-1}$$

(30)

Since every Chebyshev Polynomial T_j(ω) in contained in the interval [−1, +1], the maximum errors are given by (see also (21) and (25)):

$$\max \left|{\epsilon }_{s\hspace{0.17em}2\hspace{0.17em}i}\right|\simeq \left|{\alpha }_{3\hspace{0.17em}i}\right|=2\cdot J_3\left(\frac{\pi }{8}\cdot h\right)\cdot \mathrm{Cos}\left[\frac{\pi }{4}\cdot \left(i\cdot h+\frac{h}{2}\right)+{\phi }_{LSB}\right]$$

(31)

$$\max \left|{\epsilon }_{c\hspace{0.17em}2\hspace{0.17em}i}\right|\simeq \left|{\beta }_{3\hspace{0.17em}i}\right|=2\cdot J_3\left(\frac{\pi }{8}\cdot h\right)\cdot \mathrm{Sin}\left[\frac{\pi }{4}\cdot \left(i\cdot h+\frac{h}{2}\right)+{\phi }_{LSB}\right]$$

(32)

In this case, the worst-case error is obtained for the approximation of the sine function in the first sub-interval (i = 0). Therefore, the overall maximum approximation error due to piecewise quadratic approximation (ε_alg) is given by:

$${\epsilon }_{alg}=2\cdot J_3\left(\frac{\pi }{8}\cdot h\right)\simeq \frac{1}{192}{\left(\frac{\pi }{4}\cdot h\right)}^3$$

(33)

where we have used the approximation:

$$J_n(x)\simeq \frac{1}{n!}\cdot {\left(\frac{x}{2}\right)}^2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}x\ll 1$$

(34)

By proceeding in a similar way for a piecewise linear approximation, we find:

$${\epsilon }_{alg}=2\cdot J_2\left(\frac{\pi }{8}\cdot h\right)\simeq \frac{1}{16}{\left(\frac{\pi }{4}\cdot h\right)}^2$$

(35)

Table 2. Achieved ε_alg by varying the number of sub-intervals (LSB_y = 2⁻¹⁵). The values with grey background are the ones that were used for the actual implementation.

Number of Sub-Intervals (2H)	Piecewise Linear Approx	Piecewise Quadratic Approx
	εalg/LSBy	εalg/LSBy
4	78.893	1.291
8	19.735	0.161
16	4.935	0.020
32	1.234	2.52 × 10−3
64	0.308	3.15 × 10−4
128	0.077	3.94 × 10−5

reports the achieved ε_alg by varying the number of sub-intervals (2^H). Please note that according to DSP data-sizes, we plan to implement a DDFS with 16-bit outputs with a weight of the less-significant bit LSB_y = 2⁻¹⁵. For this reason, the error values reported in the table are normalized to LSB_y. By looking at reported values, it can be observed that if we would like to achieve a ε_alg lower than 1/4 LSB_y, the piecewise linear approximation requires 128 sub-intervals. On the other hand, using the piecewise quadratic approximation developed in this paper only eight sub-intervals can be used; this reduces significantly the size of the arrays needed to store the coefficients (28) and (29).

Figure 6 shows the implementation scheme for the sine computation block of the DDFS shown in Figure 5 (a similar scheme applies also to the cosine block), where Figure 6a considers the piecewise linear technique and Figure 6b focuses on the piecewise quadratic approach. Signal representation is reported using the notation U_m,n and Q_m,n for unsigned and signed representations, respectively. It is assumed that in both cases, 2^m represents the weight of the most-significant bit, while 2⁻ⁿ represents the weight of the less-significant bit. Please note that the number of bits of all signals is chosen in order to meet the fixed-point hardware operators of the DSP (that is 16-bit input and 32-bit output for the multipliers and 32-bit input and output for the adders). If, for the sake of simplicity, we focus on the piecewise quadratic implementation of Figure 6b, since the number of sub-intervals is chosen equal to 2^H, the H = 3 most-significant bits of z can be used to obtain the index i which identifies the sub-interval; the remaining bits of z, represents z_i = z − i ∙ h. The data-size of z is chosen in order to meet the 16-bit size of z_i, since this signal is the input of the following multipliers. The coefficients a_{1 i} and a_{2 i} are stored on 16 bits, resulting in inputs of the multipliers. In this case, the unsigned representation is chosen since a_{1 i} and –a_{2 i} results always positive. The coefficients a_{0 i}, being required only at the input of an adder, can be stored on 32 bits, in Q_0,31 representation. A final rounding operation is used to provide the output in Q_0,15 representation.

The results of the complete error analysis of the two DDFS implementations (using piecewise linear or piecewise quadratic approximations), are reported in

Table 3. DDFS error analysis and simulation. The total errors columns are highlighted in grey.

	Maximum Absolute Errors Normalized to LSBy
	Piecewise Approx. (εalg)	a0 i Coeff.	a1 i Coeff.	a2 i Coeff.	z Signal	z2 Signal	Inputs of the Adder in Figure 6	Adder Output in Figure 6	Total Error	Total Error (Simulat.)
LSB Reduction Method	-	Rounding	Rounding	Rounding	Rounding	Trunc.	Trunc.	Rounding	-	-
Piecewise linear technique	0.077	7.63 × 10−6	0.002	-	0.002	-	1.53 × 10−5	0.500	0.581	0.576
Piecewise quadratic technique	0.161	7.63 × 10−6	0.031	5.96 × 10−8	0.025	0.005	3.05 × 10−5	0.500	0.722	0.684

. It is worth noting that rounding on z signal is implemented by initializing the phase accumulator in Figure 2 to the rounding constant 0.5 ∙ 2^−P. It can be observed that some error components are almost negligible, and this allows us to employ the simpler approach of truncation in place of rounding both at the output of the multiplier computing z² in Figure 6b and at the inputs of the adder in Figure 6a,b. The theoretical overall maximum absolute error of piecewise linear DDFS is 0.581 LSB_y, while in the case of piecewise quadratic approximation the error is slightly higher (0.722 LSB_y). In both cases, the overall error close to 0.5 LSB_y confirms the accurate sizing described in Figure 6. We have also simulated the two DDFS implementations and the simulated values are reported in the last column of Table 3. It can be observed that in the case of piecewise quadratic DDFS, the simulated maximum error is slightly lower than the theoretical one; this is due to the fact that the theoretical error components considered in Table 3 tend to compensate for each other: as an example the error of a_{1 i} coefficient is maximum for maximum z_i and is positive while in the same condition, the error due to piecewise quadratic approximation is maximum in modulus but negative. In the following section, we will consider the simulated values of the DDFS errors.

If we call ε_DDFS the total maximum absolute value of the DDFS error (reported in the last column of Table 3), the overall maximum absolute error of the whole elaboration of Figure 2 can be written as:

$${\epsilon }_{elab\hspace{0.17em}total}=2\cdot {\epsilon }_{DDFS}\cdot \max \left|x(n)\right|+{\epsilon }_{FIR}+0.5\cdot {\mathrm{LSB}}_{\mathrm{y}}$$

(36)

where ε_FIR represents the total error of the FIR filter used to implement the Hilbert transform and the final contribution of 0.5 LSB_y is due to the final rounding operation needed to obtain the final output y(n) in Q_0,15 representation. Please note, however, that this simple relationship describes a worst-case behavior: as an example, if we look at (31) and (32) the sum of the absolute errors due to piecewise approximation at the two outputs of the DDFS is equal to $\sqrt{2}$ times the error of the single W output, not to two times this error.

Since the only approximation in Figure 4 is the final rounding, we have that also ε_FIR is equal to 0.5 LSB_y. As a consequence, ε_{elab total} can be written as:

$${\epsilon }_{elab\hspace{0.17em}total}={\mathrm{LSB}}_{\mathrm{y}}+2\cdot {\epsilon }_{DDFS}\cdot \max \left|x(n)\right|$$

(37)

For a worst-case input signal, with maximum value equal to 1, ε_{elab total} is therefore equal to 2.15 LSB_y and 2.37 LSB_y for piecewise linear and quadratic DDFS, respectively. The error becomes lower for an input signal that is not at full-scale. As an example, for a maximum absolute input value of 1/3, the errors for piecewise linear and quadratic DDFS implementations becomes 1.38 LSB_y and 1.46 LSB_y, respectively. These theoretical values describe again a worst-case condition since error components in the real system tend not to be at the maximum at the same time or tend to cancel each other out. As an example, by simulating the whole algorithm using a random input signal with 1/3 maximum absolute value, the simulated maximum absolute error results were 1.13 LSB_y and 1.18 LSB_y, respectively.

2.3. Experimental Auscultation Tests

The actual performance of the proposed frequency shift enhancement was assessed by carrying out specific auscultation tests. To this aim, 41 participants (23 males, 18 females, age 34.3 ± 11.6) not skilled in medical auscultation and with no history of hearing impairment were enrolled. Each participant was asked to listen to short audio clips consisting of a heart sound recording (from the public database “Classification of Heart Sound Recordings: The PhysioNet/Computing in Cardiology Challenge 2016” [71]) corrupted by noise, either unprocessed, or processed with the proposed frequency shift enhancement. Three noise sources were considered, namely additive white Gaussian noise (AWGN), crowded avenue noise (publicly available at [72]), and helicopter noise (publicly available at [73]). For each noise source, several noisy audio clips were created by corrupting each heart sound recording, either unprocessed or processed via the proposed algorithm, with different noise levels, so as to obtain values of signal-to-noise ratio (SNR) ranging from 0.025 to 1 (steps of 0.025 up to 0.3, and steps of 0.1 from 0.3 up to 1). In the case of processed heart sounds, two frequency shift (f_shift) values were considered, namely 50 Hz and 100 Hz. For each noise source, three sets of 19 audio clips were obtained, namely unprocessed (referred to as “no_elab”), processed with 50 Hz shift (referred to as “elab_50”), and processed with 100 Hz shift (referred to as “elab_100”). Therefore, each participant listened to 171 audio clips, divided into three sessions (one for each noise source), with 5-min breaks in between. In each session, the sets of unprocessed and processed audio clips were administered separately, with a random sequence of SNR values. The participants, who were blind to the specific set being administered, were asked to acknowledge their ability to reliably perceive the heart sounds; the lowest SNR (SNR_min) value for which a participant confirmed the ability to hear the heart sounds was considered to be the minimum hearing level for the related set and noise source. At the end of the experimental auscultation tests, for each participant, the SNR_min values for no_elab, elab_50, and elab_100 sets of clips were obtained for each noise source (i.e., AWGN, crowded avenue, helicopter).

2.4. Statistical Analyses

Statistical analyses of the results of the experimental auscultation tests were carried out to verify if the proposed method enhanced the participants’ ability to perceive heart sounds in the presence of different kinds of noise. An enhancement in hearing ability for a specific noise source would result in a reduction in the SNR_min obtained for the processed clips (elab_50, elab_100) with respect to the unprocessed clips (no_elab). Considering that the effects of the proposed frequency shift could change between the participants, specific tests were performed to assess if the differences between the SNR_min values of processed and unprocessed clips were statistically lower than 0. In particular, for each noise source, the Wilcoxon Signed-Rank Test was applied to SNR_min matched pairs of elab_50 and no_elab clips, and to SNR_min matched pairs of elab_100 and no_elab clips.

3. Results

3.1. DSP Implementation Results

The algorithm presented in the previous section has been described in the C programming language and optimized for implementation on Texas Instruments C67xx DSPs (Texas Instruments, Dallas, TX, USA).

Table 4. DDFS implementation results (C67xx DSP—ideal memory model).

Sin/Cos Computation Technique	LUT Storage Tech.	LUT Size (Bytes)	Code Size (Bytes)	Total Memory Size (Bytes)	Percentage BSS Section Occupation	Execution Time (Clock Cycles)
Piecewise Lineax Approx.	far alloc.	1536	224	1760	-	25.0
Piecewise Lineax Approx.	near alloc.	1536	224	1760	4.7%	19.0
Piecewise Quadratic Approx.	near alloc.	128	256	384	0.4%	21.0
Library functions—Single-Precision	-	-	928	928	-	177.0
Library functions—Double-Precision	-	-	1120	1120	-	346.5

reports the achieved performance by considering the DDFS implementation. In this case, the execution time was obtained by simulating the code on a DSP simulator with an ideal memory model (i.e., without considering cache). It can be observed that a total of five implementations were considered. The first implementation used the piecewise linear approximation technique described in the previous section and exploited a “far allocation” approach to implement the Lookup Table (LUT) needed to store a_ji and b_ji coefficients. The second implementation still used a piecewise linear approximation technique, but in this case, the LUT employed a “near allocation” technique, i.e., it was stored in the BSS data section (as well-known, this approach reduces the time needed to access the coefficients). The third implementation used the piecewise quadratic approach described in the previous section; in this case, only “near allocation” was considered since the required LUT size was very small (only 128 bytes). The last two implementations used the DDFS architecture of Figure 5 and computed the required sine and cosine functions using the optimized floating-point functions available in the FastRTS Library provided with the DSP. In this case, two implementations were considered, one using single-precision sine/cosine functions, and one using double-precision functions. Considering the total memory size required to implement the DDFS, it can be noted that the best solution is the piecewise quadratic implementation, which requires only 384 bytes. The piecewise linear implementation resulted in the largest required memory size, which is dominated by the required LUT size (1536 bytes). In addition, the piecewise linear DDFS using near allocation required also a non-negligible 4.7% of the total BSS section size, which was used precisely to store variables using “near allocation”. By considering, on the other hand, the execution time, the fastest implementation was the piecewise linear DDFS using “near allocation”, which was executed in only 19 clock cycles. Also remarkable is the performance of the piecewise quadratic implementation, which was executed in 21 clock cycles and turned out to be also faster than the piecewise linear DDFS with “far allocation”. The two implementations using library functions, with single and double precision, were much slower than other implementations, requiring 177 and 346.5 clock cycles, respectively. Please note that this last value is not an integer because the execution time turned out to also be data-dependent. It is also interesting to observe that the proposed fixed-point implementations outperform standard floating-point implementations using FastRTS Library, despite the DSP being able to implement floating-point arithmetic in hardware.

Table 5. Full elaboration algorithm implementation results (C67xx DSP—ideal memory model). The columns related to total memory size and execution time are highlighted in grey.

Sin/Cos Computation Technique	LUT Storage Tech.	Memory Size (Bytes)						Execution Time (Clock Cycles)
		DDFS LUT	FIR LUT	DDFS Code	FIR Code	Code Other	Total	DDFS	FIR Sum-of-Prod.	FIR Buffer Update	Other	Total
Piecewise Lineax Approx.	far alloc.	1536	84	224	64	640	2548	25.0	28.0	46.0	29.0	128.0
Piecewise Lineax Approx.	near alloc.	1536	84	224	64	640	2548	19.0	28.0	46.0	29.0	122.0
Piecewise Quadratic Approx.	near alloc.	128	84	256	64	640	1172	21.0	28.0	46.0	29.0	124.0
Library functions—Single-Precision	-	-	84	928	64	640	1716	177.0	28.0	46.0	29.0	280.0
Library functions—Double-Precision	-	-	84	1120	64	640	1908	342.0	28.0	46.0	29.0	445.0

summarizes the performances achieved by considering the whole elaboration algorithm described in the previous sections, by considering the five possible DDFS implementations, also outlined in Table 4. As in Table 4, a DSP simulator with an ideal memory model was considered. By looking at the execution time of the FIR filter, it can be observed that the performances achieved are coherent with the approach used to realize and optimize the filter, which is using a linear buffer algorithm for filter input and implementing the sum-of-products with a software-pipelined loop using 2× unrolling (the considered filter order is 40). Overall, the data presented confirm that the piecewise quadratic implementation resulted in the best tradeoff between memory size and execution time, with a required memory size of about 1.1 KB and a total execution time of only 124 clock cycles for the whole algorithm. The fastest implementation is again the one using the piecewise linear technique with “near allocation”, which achieves an execution time of 122 clock cycles, but requires a total memory size of about 2.5 KB. Implementations using FastRTS Library functions, with single- or double-precision arithmetic, are the slowest, with execution times of 280 and 445 clock cycles, respectively.

To test the developed algorithm, we used a DSK6713 rapid prototyping board, shown in Figure 7. The board is equipped with a Texas Instruments C6713 DSP, which uses a 225 MHz internal core clock frequency and provides a 256 KB internal memory. The board also includes a 16 MB external SDRAM, 512 KB flash memory, and a TLV320AIC23B (AIC23) Audio CODEC. In our experiments, the analog audio input is provided to the Codec using the “line in” 3.5 mm input jack, and the elaborated audio output is obtained using either the “line out” or the “headphone” 3.5 mm output jacks.

The minimum sampling frequency that can be configured in the Codec is 8 kHz; however, in the considered application, a sampling frequency of 2 kHz is more than adequate to sample cardiac audio signal (in fact f_s = 2 kHz is the design value considered in the previous section). Therefore, in the realized tests, the Codec is configured for an 8 kHz sampling frequency, and digital down and up sampling is implemented on the input and output signals, respectively, in order to keep a sampling frequency of 2 kHz for the core elaboration algorithm (which implement the frequency shift operation). Anti-aliasing and reconstruction filtering of digital down and up conversions is realized using a 39th-order FIR filter.

To realize the experiments, the code has been implemented using Code Composer Studio developing environment of Texas Instruments. The communication between the DSP and the AIC23 Codec of the board is realized using the McBSP interfaces of the DSP and employing a block processing approach with a block size of four samples. The EDMA unit of the DSP is used to transfer the blocks of data between the McBSP interfaces and the internal memory of the DSP, using a ping-pong buffering approach. The block size (equal to 4) has been chosen in order to match the ratio between the Codec sampling frequency (8 kHz) and core elaboration sampling frequency (2 kHz): this results in keeping a stream processing elaboration approach for core elaboration. Since the DSK6713 has on-board JTAG Emulation, the board is connected to the PC running Code Composer Studio by simply using a USB cable. Code Composer Studio is used for debugging and profiling. Realtime audio input and output are achieved using the 3.5 mm jacks of the board.

The results of the profiling achieved using the DSK6713 board are reported in

Table 6. Experimental profiling results obtained using the DSK6713 board. The column related to total execution time is highlighted in grey.

Sin/Cos Computation Technique	LUT Storage Tech.	Execution Time (Clock Cycles)
		Core Elaboration (Frequency Shift)	Down/Up Sampling	ISRs	Total
Piecewise Lineax Approx.	far alloc.	128.2	221.0	63.0	412.2
Piecewise Lineax Approx.	near alloc.	122.2	215.0	62.3	399.5
Piecewise Quadratic Approx.	near alloc.	124.1	218.2	62.3	404.7
Library functions—Single-Precision	-	285.2	214.1	62.6	561.9
Library functions—Double-Precision	-	510.2	228.2	69.4	807.8

. If we compare the execution time of the core elaboration with the data obtained in Table 5, we can observe that the execution time obtained by the profiling is very close to the value obtained by the simulator in the case of piecewise DDFS implementations. In this case, in fact, we have verified that the number of cache misses results is very low. A larger number of cache misses was observed in the implementations using single- and double-precision library functions, specifically considering program memory, and this corresponds to an evident increase in the execution time, with respect to the results of Table 5. This behavior can be explained by considering the code size reported in

Table 7. Required memory size in DSK6713 board implementation. The columns related to total code size and data size are highlighted in grey.

Sin/Cos Computation Technique	LUT Storage Tech.	Code Size (Bytes)				Data Size (Bytes)
		Core Elaboration (Frequency Shift)	Down/Up Sampling	ISRs	Total	Core Elaboration (Frequency Shift)	Down/Up Sampling	ISRs	Total
Piecewise Lineax Approx.	far alloc.	928	1280	1184	3392	1620	280	70	1970
Piecewise Lineax Approx.	near alloc.	928	1280	1184	3392	1620	280	70	1970
Piecewise Quadratic Approx.	near alloc.	960	1280	1184	3424	212	280	70	562
Library functions—Single-Precision	-	1632	1280	1184	4096	84	280	70	434
Library functions—Double-Precision	-	1824	1280	1184	4288	84	280	70	434

, and observing that in the case of implementations using single- and double-precision library functions, the code size of frequency shift elaboration (1632 and 1824 bytes, respectively) represents a consistent portion of the DSP Program Cache memory (4 KB). In this case, moreover, the total code size results are also equal to or larger than the Program Cache memory size.

Please note that in this application the total execution times reported in Table 6 results are much lower than the available time for real-time elaboration, considering the DSP clock frequency (225 MHz) and sampling frequency (2 kHz) is about 225 MHz/2 kHz = 112,500 clock cycles. This calls for a possible sensible reduction in DSP clock frequency that could be beneficial to reduce power dissipation. Since the DSP clock input frequency cannot be changed in the DSK6713 board, we performed this analysis using the DSP power model provided by Texas Instruments [74]. The power analysis results are reported in

Table 8. Power Analysis results (DSP TI C6713, V_DD = 1.2 V). The column related to total dynamic power is highlighted in grey.

Sin/Cos Computation Technique	LUT Storage Tech.	Minimum DSP Clock Freq. (MHz)	Power estimation
			Avg. Istr. Executed per Clock Cycle	DSP Core Utilization (%)	Dynamic Power (mW)
					Activity Power	Clock Tree	Total Power
Piecewise Lineax Approx.	far alloc.	0.824	4.02	50.2%	0.58	1.76	2.35
Piecewise Lineax Approx.	near alloc.	0.799	4.09	51.1%	0.57	1.71	2.28
Piecewise Quadratic Approx.	near alloc.	0.809	4.11	51.4%	0.58	1.73	2.31
Library functions—Single-Precision	-	1.124	2.95	36.8%	0.61	2.40	3.01
Library functions—Double-Precision	-	1.616	2.20	27.5%	0.69	3.46	4.14

. In this table, the first column reports the minimum possible DSP core clock frequency which allows respecting the real-time constraint (C/f_clock = 1/f_s, C being the total execution time in clock cycles reported in Table 6). The power analysis is performed by assuming that in each case, the DSP is clocked using the minimum possible clock frequency. To conduct an accurate power analysis, the DSP percentage core utilization has been evaluated using the profiler of the cycle-accurate simulator of C6713 DSP; in this context, the percentage core utilization has been obtained as the ratio between the average number of instructions executed per clock cycle (computed as total number of instructions executed divided by the total cycle count, value also reported in Table 8) and 8, which represents the number of functional units of the DSP (maximum number of instructions per execution packet). The obtained dynamic power dissipation is reported in the last three columns of the table. It can be observed that piecewise techniques allow achieving a remarkable power dissipation reduction with respect to single- and double-precision implementations, with a power saving that is close to 25% with respect to single precision and 44% with respect to double precision. The lowest power dissipation is obtained using piecewise linear implementation with near allocation of the LUT.

3.2. Auscultation Tests Results

Table 9. Descriptive statistics of SNR_min values obtained for unprocessed (no_elab) and processed (elab_50, elab_100) audio clips, for each noise source.

Noise Source	Processing	Min	1st Quartile	Median	Mean	3rd Quartile	Max	SD
AWGN	no_elab	0.125	0.150	0.175	0.195	0.225	0.300	0.05038
	elab_50	0.075	0.100	0.125	0.124	0.150	0.200	0.02766
	elab_100	0.075	0.100	0.125	0.124	0.150	0.175	0.02468
CROWDED STREET	no_elab	0.200	0.300	0.400	0.423	0.500	0.700	0.15236
	elab_50	0.075	0.125	0.150	0.158	0.175	0.250	0.04087
	elab_100	0.150	0.225	0.275	0.282	0.300	0.500	0.09087
HELICOPTER	no_elab	0.100	0.175	0.200	0.259	0.275	0.700	0.15550
	elab_50	0.050	0.100	0.100	0.114	0.125	0.225	0.03309
	elab_100	0.075	0.125	0.150	0.161	0.200	0.300	0.05392

shows some descriptive statistics of the SNR_min values obtained for no_elab, elab_50, and elab_100 clips, for each noise source. For each experimental setting, Figure 8 displays the pairwise violin plots of the SNR_min values, by comparing the results obtained for processed and unprocessed audio clips. In particular, the panels of the left column compare the results related to elab_50 and no_elab clips, while the panels of the right column compare the results related to elab_100 and no_elab clips; the panels of each row show the results obtained for a specific noise source. The plots provide evidence that the distributions of the SNR_min values related to the signals elaborated with the proposed frequency shifter are different, always showing lower central location values and overall, less variability.

The results of the Wilcoxon Signed-Rank tests are outlined in

Table 10. Results of the Wilcoxon Signed-Rank Test performed on matched pairs of SNR_min values related to processed (elab_50 or elab_100) vs. unprocessed (no_elab) audio clips. The z-statistic used for computing approximate p-value (z-val) and the matched pairs rank biserial correlation coefficient (r_rb) are also reported.

Noise Source	Comparison	W	z-val	p-Value	rrb	# Pairs
AWGN	elab_50 Hz vs. no_elab	861	5.61	1.005 × 10−8	1.00	41
	elab_100 Hz vs. no_elab	861	5.60	1.057 × 10−8	1.00	41
CROWDED STREET	elab_50 Hz vs. no_elab	861	5.58	1.212 × 10−8	1.00	41
	elab_100 Hz vs. no_elab	861	5.59	1.107 × 10−8	1.00	41
HELICOPTER	elab_50 Hz vs. no_elab	861	5.60	1.066 × 10−8	1.00	41
	elab_100 Hz vs. no_elab	780	5.47	2.210 × 10−8	1.00	39

. As it can be observed, all tests show very small p-values (p < 0.001) supporting the conclusion that there is large statistical evidence against the null hypothesis that there is no difference between the median SNR_min of the elaborated and not elaborated signals in favor of the alternative that the SNR_min is higher in the case of the not elaborated signals. As higher SNR_min values indicate worse audibility, the tests confirm that the elaboration algorithm reduces the SNR_min, thus improving the audibility of heart sounds. In addition, the effect size has been evaluated by means of the matched pairs rank biserial correlation coefficient (r_rb). It always reached a value of 1, thus highlighting the high practical relevance of the research outcomes. The analyses confirm the usefulness of the proposed method, which proved capable of enhancing the audibility of heart sounds in the presence of noise.

4. Discussion

This study presented a novel algorithm to enhance the acoustic perception of heart sounds in noisy environments in real time. The proposed approach shifts the frequency spectrum of heart sounds to slightly higher frequencies (from tens of Hz up to 100 Hz) in order to achieve two conditions at the same time: moving the heart sounds spectrum in a frequency range corresponding to higher human audibility, and also moving such spectrum away from environmental low-frequency noises (e.g., vehicles engines, crowded streets), which are the ones that most affect the perception of heart sounds. The frequency shift is performed via a single-sideband modulation based on the Hilbert transform. The proposed algorithm was implemented on a VLIW DSP platform, which allowed real-time operation. At the core of the DSP implementation is a novel piecewise quadratic approximation approach for the implementation of the DDFS. This approach was compared with a piecewise linear approximation method and with library functions, in terms of required memory size, execution time, and power dissipation.

Both piecewise approximations provided a remarkable decrease of about 90% in the execution time of the DDFS as compared to library functions with single and double precision. The piecewise quadratic approximation also provided a remarkable reduction in required memory size, namely about 60% with respect to the library functions, and about 80% with respect to the piecewise linear approximation, which in fact, required higher memory resources than the library functions. Considering the whole algorithm implementation, the piecewise quadratic approximation ensured a 30% to 40% reduction in memory size, a 55% to 72% reduction in execution time, and a 25% to 44% reduction in power consumption, as compared to library functions. It also provided a more than 50% reduction in memory size, as compared to the piecewise linear approximation, as well as similar execution time and power consumption. According to these results, the proposed novel piecewise quadratic approximation for the DDFS provided the best tradeoff between memory size, execution time, and power consumption.

The actual performance of the proposed algorithm in enhancing the acoustic perception of heart sounds for a human operator in noisy environments was assessed by carrying out specific auscultation tests on 41 volunteers. Several audio clips consisting of heart sound recordings (either unprocessed or processed via the proposed frequency shifter) corrupted with different kinds of noise, at different signal-to-noise ratios, were administered to the volunteers. In the case of processed signals, two shift frequencies, namely 50 Hz and 100 Hz, were considered. The volunteers were asked to acknowledge their ability to reliably perceive the heart sounds, in order to determine the minimum SNR values for reliable audibility. The minimum SNR values annotated for all the volunteers were analyzed and Wilcoxon Signed-Rank tests were performed to verify that the minimum SNR values obtained for the audio clips processed with the proposed algorithm were statistically lower than the minimum SNR values obtained for the unprocessed audio clips. The tests confirmed that the proposed algorithm, both with 50 Hz and 100 Hz frequency shifts, enhanced the audibility of heart sounds, as it allowed the volunteers to reliably perceive the heart sounds for lower SNR values, i.e., in the presence of higher noise levels.

This study has some limitations, essentially regarding the cohort of volunteers. Indeed, the enrollment did not account for the age of the volunteers, which has some relationship with the general hearing ability of a subject and may have an impact on the actual enhancement that the proposed method could provide. Future studies will focus on additional experiments on a larger cohort of subjects, with a good representation of sex, age, and other characteristics that may impact the hearing ability of a human subject. Further tests are also envisioned to be performed on a large number of cardiologists skilled in cardiac auscultation, in order to evaluate the ability of the proposed method to maintain essential diagnostic features of heart sounds.

5. Conclusions

The frequency shift of heart sounds via single-sideband modulation based on the Hilbert transform enhances the acoustic perception of heart sounds corrupted by different levels of AWGN, crowded avenues, and helicopter noises, according to the results of auscultation tests carried out on 41 volunteers. The algorithm was implemented in real time on a Very-Long-Instruction-Word Digital Signal Processor with hardware floating-point operations support. A novel fixed-point piecewise quadratic approximation approach was proposed for the implementation of the Direct Digital Frequency Synthesizer at the core of the single-sideband modulator, demonstrating to provide the best tradeoff between required memory size, execution time, and power consumption, as compared to a fixed-point piecewise linear approximation, and optimized floating-point library functions provided with the DSP. The proposed frequency shifter could be equipped in modern electronic stethoscopes at very modest costs, and could provide valuable support in all those situations where heart sounds auscultation is fundamental to assess the health status of a patient, and has to be performed in noisy environments, as an example, in emergency medicine contexts such as roadside or helicopter rescue.