Modern Methods Used in the Complex Analysis of the Phonocardiography Signal

Sometimes we need an accurate or approximate representation of a quantity in a different form; either the quantity is given by an analytical expression or by the finite set of values. Using transformations, we usually measure the similarity of a given function with an entire class of functions depending on one or more parameters (such as the frequency in the Fourier transform) which may change continuously or discretely. The wavelet representation is given in the spacefrequency domain, opposite to the Fourier analysis that gives only a frequency representation. Compact supports of wavelets provide a space, and their oscillatory nature provides a frequency representation of a transformed function. It is clear that such representation is essential for the non stationary signal. The wavelet representation of a function has the multiresolution property, which means that it is given on several resolution scales. Details defined on various refinement levels (fine meshes) are added to the rough approximation determined on a coarse mesh. The dimension of the data set that store information about the function is considerably decreased while the most important information is not lost. This is very important for a good compression that saves storage in a system memory. A data compression is fundamental for an efficient analyzing of a signal with an extremely high noise density and a large variety of shapes (Abbas K.Abbas & Rasha Bassam, 2009).


Introduction
Sometimes we need an accurate or approximate representation of a quantity in a different form; either the quantity is given by an analytical expression or by the finite set of values.Using transformations, we usually measure the similarity of a given function with an entire class of functions depending on one or more parameters (such as the frequency in the Fourier transform) which may change continuously or discretely.The wavelet representation is given in the spacefrequency domain, opposite to the Fourier analysis that gives only a frequency representation.Compact supports of wavelets provide a space, and their oscillatory nature provides a frequency representation of a transformed function.It is clear that such representation is essential for the non stationary signal.The wavelet representation of a function has the multiresolution property, which means that it is given on several resolution scales.Details defined on various refinement levels (fine meshes) are added to the rough approximation determined on a coarse mesh.The dimension of the data set that store information about the function is considerably decreased while the most important information is not lost.This is very important for a good compression that saves storage in a system memory.A data compression is fundamental for an efficient analyzing of a signal with an extremely high noise density and a large variety of shapes (Abbas K.Abbas & Rasha Bassam, 2009).
If a compression scheme is lossless it can be always recover the entire original message by uncompressing then a compressed message that has the same total entropy as the original, but in fewer bits.

The fast wavelet transform
The approximation at the m+1 scale is given by equation 1

 
(5) As a result the equations 4 and 5 are forming the multi-resolution decomposition algorithm used in our paper.

Choosing the optimal wavelet function for the analysis of PGC
Every application which requires the using of a wavelet transform requires the adaptation or finding of a certain type of wavelet function.The first criterion is defined by the correlation factor between the wavelet function and the analyzed signal.Unfortunately the shape of the correlated wavelet function is capable to indicate only the family of wavelet functions but gives no indication regarding the number of the wavelet coefficients.To compute the necessary coefficient numbers, we used the wavelet toolbox of the Matlab program.The PCG signal is decomposed accordingly the algorithm shown in figure 1: A key criterion in determining the type of wavelet function consists in evaluation the overall reconstruction error, for the given type of signals that will appear in the application.So the evaluation must be made for signal reconstruction at the extreme in terms of spectral content.We imposed this condition for the decomposition algorithm to use a minimum number of coefficients, to optimize the convolution calculation time.Therefore the evaluation of the reconstruction error is done by using the functions from the Matlab toolbox which are shown in figure 2 (Gergely S., 2011).

Fig. 2. Evaluation of the reconstruction error
After computing the wavelet function we tested the validity of the selection by taking on all detail components of reconstructed signal and computing the reconstruction error.The goal was to find a maximum error therefore the maximum sensitivity of the wavelet function to the given input signal.The tables are presenting the overall computed errors which must be the lowest possible and the component extraction error which must indicate the highest sensitivity of the wavelet functions.
These two conditions are met by the Daubeschies db2 function.It seemed that this function did answer very well to the PCG signal which is a strongly polynomial type due the massive presence of arbitrary heart murmur signals.Due to the Matlab convention the db2 wavelet function is mathematically defined by four coefficients which are computed by using the orthogonally conditions along with a specially imposed condition which gives the smoothness of the Daubeschies functions The Daubeschies functions family plays an important role in signal analysis because it is most often used to break down signals into sub-band components.Therefore the use of Db4 function requires the computation of the scaling coefficients.The scaling equation for a wavelet function having four coefficients is given by equation ( 6): (Gergely et al. 2011) And the equation for the wavelet function is: Prof. Ingrid Daubeschies imposed a supplementary condition regarding the smoothness of the wavelet function which is given by equation 8: By replacing with the computed coefficients, we get the final form: Equations 9 and 10 are used to compute the partial multiresolution decomposition up to level seven, which is used in the characterization of the PCG signals.The wavelet transform is done by using Daubeschies 4 coefficients along with signal decimation by two.Another benefit of using the wavelet transform is due the filtering effect which reduces the bandwidth of the signal.

The time-frequency representation of PCG signals
A common situation which occurs in the analysis of the PCG signals is that the domains of signal frequency spectra are almost identical but the temporal distribution of spectral components is totally different.For this reason the frequency analysis requires a different approach to traditional methods.The main parameter that leads to differentiation of the two signals is the energy distribution in time domain, which is well evidenced by Parseval's relation given in equation 12: So the best approach is a non-stationary analysis of signal properties coupled with a timefrequency representation type.Wavelet representation is made in time-frequency domain opposed to Fourier analysis as only be effective in the frequency domain representation.The compact support offered by the wavelet transform allows the analysis by space and the oscillating character of a signal.This type of analysis is best suited for non-stationary signal type.Although apparently PCG signal is characterized by time periodicity of cardiac activity, the high frequency components of the signal spectrum are strongly marked by acoustic noise which is created by blood flow through vessels and cardiac valves.Wavelet representation of a function is characterized by multi-resolution property which means that it can be decomposed in several scales.By performing wavelet transform on a signal it results a considerably signal vector size decreases due to decimation by 2, at each level of transform.This is essential as a good signal compression to store transform in a low amount of available memory and having limited resources which is specific to the embedded type of microcomputers.The simulation software packages of the Matlab programming environment, offers a wide range of functions as a representation of time-frequency analysis but their export is generally possible only at the graphic level.The correct design of the time-frequency scalogram algorithm requires more detailed study of how the spectrum inversion occurs at the decimation by two of the wavelet high pass filter.Identifying the areas in which the frequency spectrum inversion occurs can be done by the representation of wavelet packet on several levels.In figure 4 it can be observed the inversion algorithm then takes place.
The time-frequency scalogram decomposition is actually a total multi-resolution analysis, so that the representation of time-frequency scalogram requires the correct reordering of the resulted spectrum fields.So it can be seen how the spectrum fields results from filtering G (high pass), following the decimation by 2 produces spectral inversion.Given the above observations, in figure 4 we can see another interesting phenomenon in the last decomposition, where the order of occurrence of the frequency domains follows a pattern similar to the Gray code (A.Jensen & A. la Cour-Harbo, 2001).Spectral range as natural frequency ordering is done by interpreting the initial results of permutations that are written in binary 000, 001, 011, 010, 110, 111, 101, 100 (the string of numbers marked in green).An autonomous medical device witch is using a graphical color display, showing a timefrequency scalograme as a result of multiple wavelets transforms, is a very complex problem.To solve the implementation on a microcontroller, the specific problems associated with time-frequency representation have following steps: a. Correct scaling of the calculated results of multi-resolution wavelet transforms on both x, y coordinates.b.Calculation of the samples power average coefficient scaling.c. Construction of time-frequency plane by allocating the appropriate memory addresses for each levels of decomposition.d.Correction of the natural order of frequency representation.Due to decimation by 2 at each level of decomposition, the filtered frequency domain appears in the mirror.e. Setting the calculated energy values to the equivalent RGB colors.f.Converting data from the memory allocated time-frequency representation to a bxp type file.We assigned this type of graphics file to the DEM128160 display.The graphic conversion program between bmp files and bxp files were done by using the graphics library of the Embarcadero C + + program.
The time-frequency plane is used to describe the way in which the signal energy is distributed along the signal.The time-frequency scalogram is actually a full decomposition when analyzing multi-resolution, so that proper representation requires reordering the domain of the resulted spectrum.Spectral areas ordering by following the natural frequency, is done by interpreting the initial permutations results which are written in binary 000, 001, 011, 010, 110, 111, 101, 100.The new frequency order follows a Gray code pattern.It may be noted that Gray property type permutations is changing a single bit to change from one value to another.The computation of the time-frequency plane requires reiterating the wavelet transform routines for 254 times.The samples numbers of the original signal are decreasing to half after each transformation.Therefore the total amounts of required operations are 32768*7=229376 MAC instructions.At the level of the dsPIC30F6012 microcontroller the overall computing time is approximately 0.96 seconds.
The displaying of the scalogram requires the construction of a special algorithm which is capable to calculate the addresses allocated to each pixel.The address of the first pixel is programmed for the upper left corner.This reference address position is convenient for standard graphics display and alphanumeric character set implementation. In Adding a further level of decomposition, increases the frequency resolution by increasing the number of bands, but also it occurs a lower temporal resolution by reducing the number of samples.This is a consequence of the principle that states that a time-frequency representation has an intrinsic limitation, as the product in time and frequency resolution is limited by Heisenberg's uncertainty principle which expresses the equation 14: Where t  and f  are given by equations 15 and 16: Where  is the Fourier transform of base wavelet function  In figure 6 it is showed a few time-frequency scalograms obtained by using of the multiresolution algorithm presented in figure 5.
Fig. 6.Time-frequency scalograms of the acquired PCG signals

The Shannon multirate entropy in computing of the Euclidian distance of the PCG pathology
In information theory, entropy is a measure of the uncertainty associated with a random variable.In this context, the term usually refers to the Shannon entropy, which quantifies the expected value of the information contained in a message, usually in units such as bits.
Equivalently, the Shannon entropy is a measure of the average information content one is missing when one does not know the value of the random variable.In the wavelet transform theory the Shannon entropy is the measure for a cost function as how much decomposition has to be made to attain a full decomposition.After the signal decomposition, the level iterations contain a specific narrow signature reflected in the information content in a specific frequency band.At the end of the decomposition the signal is characterized by a number of Shannon's entropy coefficients.The first step in the characterization of the PCG signals is to compute the multi-resolution decomposition of the signal by using of the wavelet transform.The level at which the decomposition stops is given by the optimization is the number of calculations and the content of information of on levels below a certain threshold where it no longer pays a further decomposition.The multi-resolution wavelet decomposition shown in figure 7 presents a PCG signal with pathology.These representations are the result of simulations made under the National Instruments -LabWindows CVI platform without using any embedded software packages for the wavelet transform.
Such a signal can be represented and stored using only 14 coefficients by replacing each line in the matrix with 2 coefficients obtained by computing the Shannon entropy of each individual level of decomposition.These two coefficients are taken from the above figure.
By simulation of the calculation algorithms directly in C language, it was significantly reduced the time required to implement the DSP programs.For this reason, the code generation for the DSP microcontroller was used in another C compiler but with the same source code optimized for a different level of precision.Quantization errors which are inherent in floating point calculations do not alter the results because all evaluations were made by identifying only the minimum values of Euclidian distance.The Euclidian distance is computed by using equation 18: The signal coefficients are stored in the SD memory card in specific vectors (address) which represent the comparison reference for the recognizing process of the specific pathology.By using the Shannon entropy coefficients there is not necessary to store the entire pathology signals.These coefficients represent the compressed form of the pathology signals.The pathology signals have therefore specific spectral fingerprints which are used to compute the Euclidian distance between the stored reference vectors and the currently analyzed signal.By using the wavelet transform in the analyzing of PCG signals it is possible to compress and to preserve all time-frequency characteristics of the signals.On the other hand either time domain or frequency domain analysis does not fully describe the nature of non-stationary signal.A pathological PCG signal is dominated by the high frequency components along with the low frequency S type pulses.The statistical characterization method is usable only for a primary signal classification.A precise PCG signal characterization done with the intention of pathology recognizing, is possible only if a reference signal is in fact compared with the fully or partially acquired signal.The frequency content in the multilevel wavelet transform may well be evaluated by the information content of each level defined by the Shannon entropy which is presented in the following equation: The estimation of the signals envelope as a final characterization is a difficult task that involves intensive computing resources.Therefore the algorithm was designed to be implemented on a device using DSP engine for the signal processing.143

Correlation of the PCG signals
The correlation algorithm is using only the approximation coefficients of the wavelet transform.So that the question is which decomposition level is to be used for the signal correlation algorithm?The envelope of the signal must contain as much as possible amount of samples to get a specific pathology characterization.The property of the multilevel wavelet transform is that, each computed approximation coefficients do contain both the approximation and detail coefficients of the next level.The envelope of the PCG signal is computed by using the third level of signal decomposition which is shown in figure 8: After all the issue is to test the correlation between the acquired signal and all reference signals that are stored in the SD card.Our studies showed that the heart beat rate does not affect the overall spectral density of the PCG signal.The Shannon sampling theorem has been extended to allow for sampling times which are not uniformly spaced (Gergely S., 2011;Fliege N.J., 1994).Several slightly different versions of the non-uniform sampling theorem have arisen.The differences lie in the spaces of functions being considered and the different classes of sampling times which are permitted.The theorem essentially says that a band-limited signal x(t) is uniquely determined by knowledge of its samples {an = x(tn)} as long as the sampling times {tn} occur at a rate which is on average higher than the Nyquist rate.By compressing or expanding the signal the envelope remains shift able and also the resulted frequency variations are not important.The frequency content is previously analyzed by the using of the Shannon entropy which classifies the PCG signal at each wavelet decomposition level.The complex overall algorithm is shown in figure: The above presented algorithm is capable to analyze the extremely complex structure of the PCG signal which is characterized by a signal envelope which is amplitude modulated and in the same time it is frequency modulated having specific and distributed "carrier" frequencies.The frequency band of the original PCG signal is limited extremely sharp by using a 512 coefficients FIR filter which follows the input anti-alias analog filter.The useful frequency band of the PCG signal is situated between 62-800 Hz.
The main difficulty in doing the correlation task is that, the acquired signal is never at the same HBR (Heart Beat Rate) like the reference signals; as a result the time shift will corrupt the peak value of the cross-correlation.All signal envelopes references are recorded at a known HBR, information which is stored in the file of each signal.In signal processing, auto-correlation is a measure of similarity of two waveforms as a function of a time-lag applied to one of them.This is also known as a sliding dot product or inner-product.It is commonly used to search a long duration signal for a shorter, known feature.The auto-correlation (Broersen Piet M. T., 2006). of the third level decomposition is done by using the equation: The reason for which we have used the third level of decomposition instead of the original signal is because of the necessity to reduce the number of computation operations.By using the third level of decomposition, the overall computations are reduced by a factor of 64 due to the squaring of the samplings reduction which is at a factor of 8.The using of the crosscorrelation is the only practical method to extract a certain sequence from a noisy PCG signal which is dominated by the murmur of the heart (Abdallah et al.,1988)..The resulted auto-correlation is shown in figure 10: The correct cross-correlation is complete after the re-sampling of the acquired signal in accordance to the reference signal HBR.After re-sampling, the signal is interpreted as having the same original sampling rate.The necessary re-sampling rate is given by the ratio between the HBR parameter of the analyzed signal and the HBR parameter of the reference: Where L is the interpolation factor which means that between two consecutive signal samples there are inserted a number of L-1 zeros.M is the decimation factor which is always constantly equal with 100.Extensive simulations have proved that the algorithm is permissive to an R RAtio from 0.43 to 0.99.The reference signal is always converted to a HBR equal to 100 before it is stored in the pathology data base.As a consequence the above algorithm has the drawback of using a lot of memory space in the systems memory.
The sensitive part of the overall algorithm is the interpolator used in the multirate sampling module.It turned out that the final correlation operation is obviously sensitive to the unmatched envelope.The issue is the lag over the real envelope of the computed values in case of a under filtering or a diminishing and distorted wave form in case of an over filtering.Therefore it was necessary to calculate the interpolation filter.The first step of using a linear interpolator has the ability of good low frequency attenuation but on the other hand the high frequency components of the envelope signal are slightly distorted.Finally we have preferred a FIR filter.The usually used two filters for the re-sampling process is reduced to a single filter which is chosen by using equation 23; (Fliege N.J., 1994) min , At the third level of signal decomposition the new sampling rate is interpreted as 1000 Hz The frequency domain which must be covered by the filters is shown in figure 11: The impulse response of the low pass filter is given by equation 24: The filter coefficients are weighted with the Hamming window: 0.54The number of filtering coefficients is computed with equation 6.8 The cutting frequency does not depend on the interpolation factor L and is given by: 300 2 For HBR between 43 -99 Where HBR ref =100 It yields for L interpolation factor between 43 to 99 a number of 355 to 817 coefficients.
For the above mentioned R Ratio interval we got a number of coefficients between 355 and 817.The result of the re-sampling process is shown in the figure 12.These values aren't large in comparison to the permanently used 512 coefficients in the overall band pass filter.All these coefficients are stored in the program memory of the DSP microcontroller.

Computing of the PCG signals envelope
Hilbert transform is a process which is capable to extract the precise envelope of a given signal.Generating a signal in the complex domain and out of phase by 90 degrees from the real signal provides a series of numerical processing such as quadrature modulation and demodulation of signals, and the implementation of automatic gain control systems.Time and frequency representation shows how the Hilbert transform cosine signal spectral component rotates with -j and the negative component by +j.An important property of the Hilbert transform is that it is a theoretical system with frequency response magnitude one and phase equal to 90 degrees for all frequencies.This means that a signal passing through a Hilbert system will be weighted by 2 the amplitude and phase will be modified by ¼ the period T. If a real signal x r (t) is modulated in amplitude, then the modulated signal envelope is that it contains the useful information.So the instantaneous envelope E(t) becomes: The above equation is the envelope of the modulator signal and is used to calculate also the PCG signal envelope.Traditional signal demodulation of amplitude modulated signal consists of previously rectifying the signal by squaring and applying a smoothing low pass filter for the carrier frequency.This type of demodulation has been tested for PCG signal envelope but a comparison shows that the voltage ripple obtained by Hilbert transform is clearly in favor.The filter ripple for both methods is shown in figure 13: In terms of the number of calculations required, the Hilbert transform method is faster because it involves the using of a filter with a lower number of coefficients.
Implementation of Hilbert transform involves the calculating of the impulse response of a system considered linear.So that for arbitrary signal:   For a discrete signal the above result has to be adapted, therefore the final form is given by equation 34:  Due to the complex structure of the PCG signal which is characterized by the presence of a large number of distributed peak values, it was necessary to introduce supplementary information regarding the amplitude of the signal.During the researches we applied a new method to analyze the PCG signal by introducing an additional coordinate to the signal vector which is visible in picture 17.The idea is that we switched to two-dimensional analysis similar to pattern recognition located in an image (Andrew K. Chan & Steve J.Liu, 1998).Thus the points defined in terms of measured values X0Y are introducing a second dimension of the signal and can thus be interpreted as an image.Reducing the number of calculations required to scale the amplitude vector x[t] for a number of dots (pixels).This type of resulted image can be processed by using advanced image processing methods.On a closer analysis, there is an important feature which allows to limit one of the coordinates of the vector w[y][x].It can be observed that both signals involved in the algorithm are scaled to the same number of values on the Y coordinate, which provides a significantly higher computing speed.This was imposed because the simple convolution of the reference signal with the analyzed signal in order to get the pathology correlation is not efficient.The solution to the problem was to convert the one dimensional signal vector to a twodimensional vector which at the end represents a graphical image of the PCG signal (Anke Meyer-Base, 2004).The conversion starts with the dividing of the sample amplitude to a given number of desired vertical pixels.If the value is in the computed domain it will be replaced by the value of 1, or if the value is outside of the domain it becomes zero.The signal conversion process is shown in figure 17.

Fig. 17. Two-dimensional conversion of the PCG signal
By extending the method it is possible to give a different magnitude to different image segments by inserting different values instead of ones.As a result the image performs as a picture having a new z coordinate.The modified method makes possible the transformation from absolute pathology detection to a probable type of detection which is useful for primary pathology classification.This way the PCG signal image of the envelope is processed by using conventional image processing algorithms.
The extraction of the PCG signal references is done using a custom designed program which is fully interactive and provides the necessary support for the archiving of the results.All extracted reference vectors are stored in the SD card of a hand held device.The amount of necessary memory to store a reference image does not exceed 2KB.Some reference images are shown in figure 18: Under these conditions the reference image is represented in a matrix of 50x350 points.To test the recognition algorithm, the routine uses computer generated pixels of the image without intensity parameter Z (pixel has value 1) and no interpolation therefore a minimum number of points.The lowest number of points ( 42) is contained by the reference of normal heart rhythm.Analogous to the procedure for the reference signal, it applies the same above procedure for the analyzed signal.As a pattern recognition method, we have used the computation of the Euclidian distance between the analyzed signal and the stored references.

Construction of a PCG signal analyzer device
The equipment's operating system (Gergely S., 2011) is programmed to handle internal data flow between its component microcontrollers, communication with the outside through the USB port and graphics system.An important component of the program management module for distributing data is between the two microcontrollers.The operating system was designed to allow on board test routines running for the speed performance evaluation of the development system.The hardware system connects to a PC via USB2 interface which is managed by its own program, designed specifically for this application.The host program was implemented in LabWindows CVI which is capable to load into RAM or SD card any routine to emulate the original data in the DSP microcontroller.The resulting data are then transmitted back to the PC for evaluation.This was required because between the programs which are running on your PC and development board are differences in the computing precision.While DSP microcontroller does run routines with the precision given by organizing the data in Q1.15 format, all written routines in LabWindows CVI for the initial simulation are in double precision.The hardware uses a dsPIC30F6012 microcontroller with built-in DSP unit and a PIC18F4550 microcontroller equipped with built-in USB (Gergely S., 2011;Axelson Jan, 2006) full speed mode.
The hardware set of elements which are entering in the composition of the DSP module are specialized for rapid calculation of the amount of products.This equation is fundamental in digital signal processing and mathematical point of view is based on the convolution operation.The development is designed so that it carries out all processes signal analysis using equation 35: Recursive FIR filter implementation can be done easily by using modulo addressing type programming which after computing the values; it provides the return address used for the filter coefficients.Hilbert transform calculation assumes a similar algorithm but to simplify modulo addressing type, due to the small amount of Hilbert coefficients the zero values are inserted in the looping list.
The DSP microcontroller has two direct hardware implemented instructions to compute the Euclidian distance with or without accumulation (ED and EDAC).The minimum value for the distance measured to all references represents the match with the analyzed signal.The data stored on the SD card is downloaded via the USB.
The experiments were made using the internal 12 bit ADC converter.The final objective was to construct an algorithm capable of discriminate between certain pathological situations.The measured frequency band of the phonocardiography signal is between 62Hz to 800Hz.The chosen sampling frequency is a standard 8 KHz.As a result the sampling frequency is much higher than the required Nyquist criterion.This way the aliases are pushed far from the used frequency band.The over-sampling is the best way to design a low order anti-alias analog filter for the input.After the acquisition process the signal is normalized to an amplitude of ±1V or in other words is converted to the Q15 binary signed fraction representation which means 0x7FFF (32767) to 0x8000 (-32768).This way any multiplying with the normalized filtering coefficients does not exceed the value of 1. Obviously because the ADC has only 12 bits the over range risk is less likely.
During the acquisition of the PCG signal the samples are previously filtered with a 512 taps FIR band pass filter to cut under 60db the low frequency noises generated by the patient moving along with the operators hand moving.The presence of the FIR filter is crucial due to the necessity of a band limited signal required by the wavelet transform.The DSP engine is well suitable for the fast computing of the required convolution between the acquired signal and the filtering coefficients.This is a similar process to any FIR or IIR filter which is a familiar in digital signal processing (Emmanuel C. Ifeachor & Barrie W. Jervis, 2002).The wavelet filter uses the hardware implemented MAC (Multiply Accumulate) instruction.Usually by design the FIR filter coefficients are symmetrical, having one or two middle values identical.In the case of the wavelet filter a special attention requires the correct order of indexing the filtered signal coefficients, because of the non-symmetrical coefficient values.In case of a wrong order or array flipping, the process of the convolution ends with a crosscorrelation with unexpected output values.This is happening because the two DSP processes are mathematically related.Thus the memory address of the hardware DSP circular buffer of the coefficient data memory space has to be initiated properly.The wavelet filtering process starts with the saving of the stored data to a SD memory card, card preferable FAT formatted for direct file writing-reading.Each filtering iteration is done first for the approximation coefficients and secondly for the detail coefficients.After the iteration is done the sample number is half of the original signal or the previous approximation coefficients.The results are stored back to the memory card and this way the DSP processes may possibly be tested through the systems USB connectivity by an external PC.The overall necessary processing time for a scale index equal with 7 but only on the wavelet tree is equal to 12.5s, including al the read and write time to the SD memory card.In case of using a large static RAM instead of a SD memory card the overall time goes down to 25ms.The system is equipped with a color LCD display therefore is capable of displaying the signals coefficient scalogram in a complex time-frequency representation.The Daubeschies 4 coefficients are used in a similar manner to compute the wavelet transforms, but the circular buffer is addressed in a way that the resulted signal coefficients are decimated by two.
For a much shorter name using, the two microcontrollers are called "north" and "south".These two units are controlled by an operating system which is embedded in the south microcontroller.The south microcontroller is slower but it has the USB communication module.This module makes possible the firmware update via the USB port or a fast data exchange between the device and the host PC program.The structure of the operating system is presented in figure 20: The above mentioned operating system is an ideal structure which makes possible a separate development of the DSP computing routines.At the end the final program structure can be embedded in one single DSP microcontroller.In that case the system requires an external USB hardware.All resulted data and graphics are displayed on a 1.8" color LCD.This 128x160 pixels display is sufficient for a small portable unit but if the user wishes all graphical data may be transferred to a host PC for a much better image resolution.The display uses a Himax HX8345-A controller which has to be programmed through the 18 bit parallel port.It turned out that the initial high color depth is useless in the application; therefore this was downgraded to 65536 color levels by reducing the bit number of the parallel port.The hardware bus configuration is shown in figure 21: For data displaying it was built a custom designed character set which is structured on 13xW pixels where the W parameter may vary from 4 to 8. By using a dynamically allocated character width, the displayed text does have a much realistic appearance.The device uses also some graphics which uses a BMP structure.The conversion of the images from 18 bit to 16 bit RGB565, required the design of a auxiliary conversion program which makes also possible the transferring of the raw bmp pictures to the SD card of the device.For intermediary data storage the device is using a 4Mb static RAM addressed in a 16 bit structure.The memory can be addressed by both microcontrollers but the data flow is managed by the south master microcontroller.In order to avoid bus conflicts the system uses two buffered data transceivers.The south microcontroller can address the RAM only indirectly because of the 19 bit port which is generated by the DSP north microcontroller.The structure of how the RAM is accessed is shown in figure 22: As it can be seen in the above picture, the DSP microcontroller may well access RAM directly, but entering into this routine is subject to the interrupt generated by the master microcontroller.Meaning the data is controlled by setting the direction of the data access switches SSW and NSW.For optimal control of data flow and to avoid transfer conflicts, the master microcontroller will set that directly or indirectly the access to RAM, by the DSP microcontroller through SW-NSW.To increase the transmission speed we preferred the use of a parallel transmission protocol instead of the SPI transmission.Of course this mode of transmission does increase complexity of the designed hardware.NORDSEL signal selects DSP microcontroller which answers after interruption by turning on SSW that generates parallel port command words transfer by activating RS.If the sent command to the DSP is to acquire an acquisition, than it can communicate directly to RAM.After transmission, the DSP microcontroller sends an interrupt to the south POE2 confirmation signal.The answer will be change south direction by activating the DSP to send other commands.Data transfer between the south microcontroller and the RAM memory is done indirectly by DSP microcontroller.
Data obtained during the execution of signal processing routines and all graphical elements are stored on a SD memory card or MMC card.In this way the available system memory is practically unlimited given the fact that the 32-bit addressing way allows access to 4GB of memory.The memory card can have access to FAT16 or as access to physical memory sectors.Because the FAT protocol is a copyright of Microsoft Corporation, yet we preferred to access the SD card memory as physical memory.Of course in this way, only the host PC can access the card by using a proprietary memory structure.The advantage of this mode of operation is having a full control of stored data without the possibility of unauthorized copying of the content.Bringing the SD card in active state and data transferring to memory is done according to a predetermined by the manufacturer protocol, which includes a series of command words.Also the speed of accessing the SPI communication protocol is critically dependent on the state in which the SD card is accessed.A series of data related to memory capacity, internal organization and memory speed are programmed in a sometimes inaccessible memory area, which is equivalent to hard disks MBR.The organization by predefined memory sector structure is similar to a hard disk.The major difference from a hard disk access speed is very limited on some models.The Best SD cards reach a speed of 20MB/sec but the slower may be under 250Kb/sec.For this reason the SD memory card is the slowest element in the operating system which was taking into account at the design of software routines.The SD card memory of the designed prototype is organized as shown in figure 23:  After sending a control packet to the card, the microcontroller will receive a response R1, R2 or R3; because data transfer is ensured by generating a clock signal, the host must always send a false 0xFF byte type to keep the communication channel active.The allocated response time after sending a command packet (NCR) of the SD card has a length of 0-8 bytes.Ultimately the CRC field is optional to check the validity of data and becomes critical only when using high-speed transmission between host and SD card.To ensure compatibility with SD cards having limited capacity, the application design does not use CRC field because the transmission rate was reduced to a rate of 250Kb.Normally, this rate of transmission can use a minimum speed 400Kb.The necessary condition for the activation of the SD card is to exceed the power supply voltage value of 2 volts then that CS and DI is set to high level for at least 170 SCLK transmitted clock periods.Immediately after activating the card, it can receive native commands.After activating the card it will receive automatically a soft reset.This sets the SPI port system clock to 100 KHz and it is send CMD0 with CS at low level for entry into reset mode.After entry into SPI mode the host will switch off the SD card verification CRC.After acceptance of commands CMD0 the card goes in idle mode.Switching to Idle mode takes up to 400ms while the card is not accessible to the host microcontroller.
The hardware of the device is divided in two parts which are separated by the independent microcontrollers.The DSP microcontroller uses an 80MHz clock frequency which makes possible a real time FIR filtering of the input signal by using of 512 coefficients in about 18μs for every sample.Obviously, that the FIR filtering is used following the anti-alias filtering of the input signal.The supplementary FIR filtering is necessary to avoid the occurrence of frequency artifacts during the wavelet filtering.The hardware structure of the device is shown in figure 25:

Conclusion
Smart use of digital processing tools helps in a new approach as concerning the classical cardiac auscultation, the innovative techniques with the ultimate goal of improvement in prevention and cardiology care.Early identification of cardiac pathologies is closely related to the physician's capacity to perceive and correctly interpret the heart sounds.By using a medical electronic device for PCG signal analysis, it can be increased the rate to identify heart abnormalities, compared with classical method of auscultation.The absolute novelty of the work consists of providing a series of original algorithms that can be used to accurately identify cardiac pathologies.This opens a new road in the PCG signal analysis said to be almost unapproachable.So far, most noninvasive type tests are extracted from the ECG signal.Transformation of PCG signal into a signal with a precise and predictable envelope classified after a pattern archive, together with a known spectral behavior makes possible a full characterization of cardiac pathology obviously within the limits set by a cardiologist professional.
Algorithms presented in this paper refer to those which can be implemented immediately in embedded analysis software.Development of new algorithms for analysis is directly linked to the hardware capabilities of the systems that are to implement, therefore the project started to show weakness when we decided to implement the convolution algorithm of the computed PCG images.The immediate problem that arose was insufficient RAM for the large number of two-dimensional vectors which had to be temporarily stored.Temporary solution was to transfer the algorithm towards the development system by computing the interpolation routine segments with intermediate storage of data on SD memory card.High memory access time of the SD card increases the calculation times to impractical values.So at this moment, the prototype system is capable at this development stage only for the Shannon entropy calculations which are necessary for the display of signal's time-frequency scalogram.The introduction of the additional intelligent decision-making algorithms to the operating system presented in the paper can be a powerful tool for characterization of cardiac pathologies by acquiring the PCG signal.Although this new medical imaging device suffers from a low spatial and temporal resolution; it could be proved to be a good choice for low-cost and mobility strategy in cardiac imaging, rather than the expensive ultrasound imaging devices.The classical auscultation technique benefits from a great quality improvement by using a device which is capable to offer a time-frequency representation.

Future directions
The expected research direction is to be guided to improve the image quality on a larger display and increase the number of wavelet decomposing levels to avoid the necessary interpolations used in present.In terms of hardware we can go two ways of solving the lack of resources.The first way is redesigning part of the current development system to increase the RAM capacity, or alternatively to satisfy also the condition of low energy demand which consist of a total redesign by migrating to the Texas Instruments microprocessors.The change of the microprocessor version implies full rewrite of the operating system.TMS320C5000 series microprocessors are low cost with fixed-point DSP capabilities and having a very low consumption in terms of supply using a voltage of 1.8 volts and an integrated memory of 320KB/16biţi.The move to Texas Instruments microprocessors is supported also by an additional argument; the TI components do have accreditation for standard medical use.Development of electronic equipment is linked to building a database of reference PCG signal images for a large variety of pathologies.The cross-correlation algorithm is although under construction but the proved interest for this device among the cardiologists, makes us confident for the utility of the newly designed features.

Fig. 3 .
Fig. 3. Equations system for solving the Daubeschies 4 coefficients Therefore in the case of Daubeschies multiresolution algorithm we get:

Fig. 4 .
Fig. 4. (a) Wavelet packet order (b) natural order Fig. 5. Time-frequency plane algorithmAfter the total decomposition the obtained frequency resolution is given by equation 13:

Fig. 7 .
Fig. 7. Shannon entropy of PCG signal decompositionEquation 17(Gergely S., 2011) represents the matrix of a PCG signal which is decomposed on up to L levels and containing N samples for the original signal.

Fig. 10 .
Fig. 10.Auto-correlation of the third level decomposed PCG signal After the computing of the autocorrelation the HBR parameter is obtained with the equation 21:

Fig. 12 .
Fig. 12.Time compression of the PCG signal by using the third level of decomposition

Fig. 13 .
Fig. 13.Voltage ripple (a) Squaring method (b) Hilbert transform Fig. 14.(a) Impulse response of the Hilbert transform (b) using of a Blackman window It is possible to use a signal envelope in the evaluation process of the PCG signal because the algorithm does not care about the frequency spectrum of the original signal.The frequency spectrum of the signal was previously evaluated in a probabilistic manner by the Shannon entropy.The Hilbert transform is software implemented by computing the filtering coefficients using the finite impulse response of the transform.Therefore the FIR coefficients are computed by replacing the n index in the above equation.Choosing an odd number of coefficients is an essential criterion because of the center tap of the signal phase condition.To improve the low frequency response of the Hilbert transform we used a supplementary Blackman window.The Hilbert transform implementation is shown in figure15:

Fig. 18 .
Fig. 18.PCG signals reference images Pathology reference images were extracted from the real signals and shows how easy it can be identified with the pathology of origin.The obtained signals does resemble with the well known EKG signals.All noises under the signals envelope were eliminated for a better pathology exploration.By computing the Euclidian distance in order to evaluate the correlation between these reference signals and the analyzed signals, finally we got a nearly perfect match which is shown in figure 19:

Fig. 20 .
Fig. 20.The structure of the operating system

Fig. 23 .
Fig. 23.Data structure of the SD card memory

Fig. 25 .
Fig. 25.Hardware structure of the PCG signal analyzerThe preconditioning of the PCG signal is completed by using a 6 th order Sallen-Key filter which has to cut all unwanted frequency harmonics below 72db at the half of the sampling frequency.That way the 12bit analog-digital converter is acquiring a clean non-aliased signal.The frequency response of the combined anti-alias and FIR filter is presented in figure26.It can be observed the outstandingly sharp, band-pass characteristic of the overall filter.

Fig. 26 .
Fig. 26.Frequency response of the signal preconditioning module

Table 1 .
Coefficients of the wavelet function 0.46 cos