Perceptual Echo Control and Delay Estimation

111

Fig. 13. Sparseness measure and its impaction on length of adaptive filter
The quantity in Eq. (56) represents an energy measure in [dBm] within an estimated impulse
response. Fig. 14 demonstrates the curve for this parameter for two frequency domain
algorithms. It can be observed that after 0.5 seconds estimation of IR energy measure for
IPMDF stops fluctuating. Consequently, this fact can be used for switching between
different adaptation schemes. Fig. 14. Estimated misalignment and energy measure for MDF and IPMDF algorithms

Adaptive Filtering Applications

112
Before representing the proposed MDF scheme, declare the following statements:
I. Utilization of sparseness measure, 0<ξ(m)<1 (for real IR)
a. ξ(m)<0.7
i. IR is considered dispersive
ii.
M
t
= L, if fully updated scheme is chosen (mostly during initial period)
iii.
L/K ≤ M
t
< L, if partial updated scheme is chosen
b.


2
,7.0










mkmkLmkmk
ttk
,,,1, φGww 












mgmgmgdiagmk
NkNkNkNt 11

lkN
mw
mw
L
g




const
k



const
LNmif




,












112
,7.0



,0,mod

Tmif









otherwise
LlmofmaximaMtoscorrespondlif
mp
lf
l
,0
12, ,0,,1





,
~
ofhalffirst,










mkmkLmkmk
ttk
,,,1, φGww








const
k





(58)
where
y(m) is a desired signal (echo) and d(m) is adaptive filter’s output. Note that any
reliable adaptive filter with disabled residual echo suppressor has to achieve ERLE of -15dB
within 1 second after starting convergence process (ITU-T G.131, 2003). Fig. 15, 16, 17
illustrate an application of
μ-based selection metric. This kind of metric is used for sub-filter
selection, when estimated sparse measure parameter, ξ, equals or larger than 0.7. This value
was defined experimentally during multiple trials for numerous types of echo path. We
suggest using
μ-based selection metric as an individual block step-size parameter. It helps
accelerating a speed of convergence by allocating larger step-size values for currently
updated sub-filters. If you look at the diagram illustrated in Fig. 17 carefully, you will notice
that the energy, which is available for adaptation, is concentrated around the sparse region
of the echo path. Thus, this fact can be used for selecting sub-filters to be updated along with
setting the step-size parameter for these sub-filters. When estimation of sparse measure, ξ, is
smaller than 0.7, we suggest switching to the
χ – based selection metric. Fig. 15. Sparse impulse response and estimated
μ-based selection metrics

Fig. 16. Dispersive impulse response and estimated
μ-based selection metrics
Fig. 18 demonstrates the normalized misalignment and ERLE parameters obtained for real
speech signals. The proposed partially updated scheme for MDF shows the similar
performance comparing to the other three fully updated frequency domain algorithms.
During our future work we are going to enhance the above described algorithm and
propose a new class of partial sparse-controlled robust algorithms, which will work reliably,

According to the ITU-T Recommendation G.168, this period should not last more than one
second. The Multi-Delay block Frequency domain (MDF) adaptive algorithm can easily
outperform all existing time domain algorithms. Moreover, taking into the account the fact
that the generalized cross-correlation algorithms operate in the frequency domain and use
advantages of the fast Fourier transform, further computational savings for the adaptive
filters are achieved in the frequency domain. Therefore, the fourth section deals with partial,
proportionate, and sparse-controlled adaptive filtering algorithms working in the frequency
domain. What we claimed, within this section, is: a new metric for performing partial
updating; a new approach for designating transitions between MDF and IPMDF-based
updating schemas; a method for estimating step-size control parameter; a new partially
updated sparseness-controlled improved proportionate multi-delay filter; all the approaches
are suitable for implementation whether in time or frequency domains. The proposed
algorithm has both a performance compared to the IPMDF and SC-MDF algorithms and
reduced computational complexity along with the adjustable step-size parameter. Although
the preferred embodiments of the proposed algorithm have been described, it will be
understood by those skilled in the art that various changes may be made thereto without
departing from the main scope of the invention or the appended claims.
6. Acknowledgment
This work was supported by the Grant Agency of the Czech Technical University in Prague,
grant No. SGS 10/275/OHK3/3T/13 and by Grant The Ministry of Education, Youth and
Sports No. MSM6840770014.
7. References
Choi, B K.; Moon, S.; Zhi-Li, Z. (2004) Analysis of Point-To-Point Packet Delay In an
Operational Network, Proceedings of INFOCOM 2004, Twenty-third Annual Joint
Conference of the IEEE Computer and Communications Societies, vol. 3, pp.1797-1807,
ISBN 0-7803-8355-9, Hong Kong, March 7-11, 2004
Gordy, J.D.; Goubran, R.A.
(2006) On the Perceptual Performance Limitations of Echo
Cancellers in Wideband Telephony
IEEE Transactions on Audio, Speech, and

Carter, G.C.
(1976) Time Delay Estimation Ph.D. dissertation, University of Connecticut,
Storrs, CT, pp.67-70
Mueller, M. (1975) Signal Delay,
IEEE Transactions on Communications,
Buchner, H.; Benesty, J.; Gansler, T.; Kellermann, W. (2006) Robust Extended Multidelay
Filter and Double-talk Detector for Acoustic Echo Cancellation,
IEEE Transactions
on Audio, Speech, and Language Processing, vol. 14, issue 5, pp.1633-1644, ISSN 1558-
7916
Youn, D.H.; Ahmed, N.; Carter, G.C. (1983) On the Roth and SCOTH Algorithms: Time-
Domain Implementations, In Proceedings of the IEEE, vol. 71, issue 4, pp.536-538,
ISSN 0018-9219
Zetterberg, V.; Pettersson, M.I.; Claesson, I. (2005) Comparison Between Whitened
Generalized Cross-correlation and Adaptive Filter for Time Delay Estimation,
Proceedings of MTS/IEEE, OCEANS, vol. 3, ISBN 0-933957-34-3, Washington D.C.,
September 17-23, 2005
Hertz, D. (1986) Time Delay Estimation by Combining Efficient Algorithms and Generalized
Cross-correlation Methods,
IEEE Transactions on Acoustics, Speech and Signal
Processing,
vol. 34, issue 1, pp.1-7, ISSN 0096-3518
Knapp, C.; Carter, G.C. (1976) The Generalized Correlation Method for Estimation of Time
Delay,
IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 24, issue 4,
pp.320-327, ISSN 0096-3518
Wilson, K.W.; Darrell, T. (2006) Learning a Precedence Effect-Like Weighting Function for
the Generalized Cross-Correlation Framework, IEEE Transactions on Audio, Speech,
and Language Processing, vol. 14, issue 6, pp.2156-2164, ISSN 1558-7916
Tianshuang, Q.; Hongyu, W. (1996) An Eckart-weighted adaptive time delay estimation

vol. 1, pp.394-398, ISBN 0-7803-5148-7, Pacific Grove, CA, November 1-4, 1998
Benesty, J.; Gay, S.L. (2002) An improved PNLMS algorithm, Proceedings of ICASSP '02, 2002
IEEE International Conference on Acoustics Speech and Signal Processing, vol. 2,
pp.1881-1884, ISSN 1520-6149, Minneapolis, MN, USA, April 27-30, 2002
Fevrier, I.J.; Gelfand, S.B.; Fitz, M.P. (1999) Reduced complexity decision feedback
equalization for multipath channels with large delay spreads, IEEE Transactions on
Communications, vol. 47, issue 6, pp.927-937, ISSN 0090-6778
Douglas, S.C. (1997) Adaptive filters employing partial updates, IEEE Transactions on
Circuits and Systems II, vol. 44, issue 3, pp.209-216, ISSN 1057-7130
Aboulnasr, T.; Mayyas, K. (1999) Complexity reduction of the NLMS algorithm via selective
coefficient update, IEEE Transactions on Signal Processing, vol. 47, issue 5, pp.1421-
1424, ISSN 1053-587X
Aboulnasr, T.; Mayyas, K. (1998) MSE analysis of the M-Max NLMS adaptive algorithm,
Proceedings of IEEE International Conference on Acoustics Speech and Signal Processing,
vol. 3, article ID 10.1109/ICASSP.1998.681776, pp.1669-1672, ISBN 0-7803-4428-6,
Seattle, WA, May 12-15, 1998
Schertler, T (1998) Selective block update of NLMS type algorithms, In Proceedings of IEEE
International Conference on Acoustics Speech and Signal Processing, vol. 3, pp.1717-
1720, ISBN 0-7803-4428-6, Seattle, WA, May 12-15, 1998
Dogancay, K.; Tanrikulu, O. (2001) Adaptive filtering algorithms with selective partial
updates, IEEE Transactions on Circuits and Systems II, vol. 48, issue 8, pp.762-769,
ISSN 1057-7130
Naylor, P.A.; Sherliker, W. (2003) A short-sort M-Max NLMS partial-update adaptive filter
with applications to echo cancellation, Proceedings of ICASSP '03, 2003 IEEE
International Conference on Acoustics Speech and Signal Processing, vol. 5, pp.373-376,
ISBN 0-7803-7663-3, Hong Kong, April 7-10, 2003
Jinhong, W.; Doroslovacki, M. (2008) Partial update NLMS algorithm for sparse system
identification with switching between coefficient-based and input-based selection,
Proceedings of CISS 2008, 42nd Annual Conference on Information Sciences and Systems,
vol.3, pp.237-240, ISBN 978-1-4244-2246-3, Princeton, NJ, March 19-21, 2008

1. Introduction
The electrocardiogram (ECG) records the electrical activity of the heart，which is a
noninvasively recording produced by an electrocardiographic device and collected by skin
electrodes placed at designated locations on the body. The ECG signal is characterized by
six peaks and valleys, which are traditionally labeled P, Q, R, S, T, and U, shown in figure 1. Fig. 1. ECG signal
It has been used extensively for detection of heart disease. ECG is non-stationary
bioelectrical signal including valuable clinical information, but frequently the valuable
clinical information is corrupted by various kinds of noise. The main sources of noise are:
power-line interference from 50–60 Hz pickup and harmonics from the power mains;
baseline wanders caused by variable contact between the electrode and the skin and
respiration; muscle contraction form electromyogram (EMG) mixed with the ECG signals;
electromagnetic interference from other electronic devices and noise coupled from other
electronic devices, usually at high frequencies. The noise degrades the accuracy and
precision of an analysis. Obtaining true ECG signal from noisy observations can be
formulated as the problem of signal estimation or signal denoising. So denoising is the
method of estimating the unknown signal from available noisy data. Generally, excellent

Adaptive Filtering Applications

124
ECG denoising algorithms should have the following properties: Ameliorate signal-to-noise
ratio (SNR) for obtaining clean and readily observable signals; Preserve the original
characteristic waveform and especially the sharp Q, R, and S peaks, without distorting the P
and T waves.
A lot of methods have been proposed for ECG denoising. In general both linear and
nonlinear filters are presented, such as elliptic filter, median filter, Wiener filter and
wavelet transform etc. These methods have some drawbacks. They remove not only noise

has been shown to have asymptotic near-optimality properties over a wide class of
functions. The crucial points are the selections of threshold value and thresholding
function. The generalized threshold function is build. Computationally exact formulas of
bias 、variance and risk of generalized threshold function are derived. Section 5
concentrates on adaptive threshold values based on EEMD.Noisy signal is decomposed
into a series of IMFs, and then the threshold values are derived by the noise energies of
each IMFs. To evaluate the performance of the algorithm, Test signal and Clinic noisy
ECG signals are processed in section 6. The results show that the novel adaptive threshold
denoising method can achieve the optimal denoising of the ECG signal. Conclusions are
presented in section7.

Adaptive Noise Removal of ECG Signal Based On Ensemble Empirical Mode Decomposition

125
2. Empirical mode decomposition
EMD has recently been proposed by N.E.Huang in 1998 which is developed as a data-driven
tool for nonlinear and non-stationary signal processing. EMD can decompose signal into a
series of IMFs subjected to the following two conditions:
1. In the whole dataset, the number of extrema and the number of zero-crossing must
either be equal or differ at most by one.
2. At any time, the mean value of the envelope of the local maxima and the envelope of
the local minima must be zero.
Figure.2 shows a classical IMF. The IMFs represent the oscillatory modes embedded in
signal. Each IMF actually is a zero mean monocomponent AM-FM signal with the following
form:

() ()cos ()xt at t


(1)

Regarding the
1
()htas the new data and repeating steps 1 and 2 until the resulting
signal meets the two criteria of an IMF, defined as
1
()ct. The first IMF
1
()ct contains
the highest frequency component of the signal. The residual signal
1
()rt is given
by
11
() () ()rt xt ct .
4.
Regarding
1
()rt as new data and repeating steps (1) (2) (3) until extracting all the IMFs.
The sifting procedure is terminated until the Mth residue
()
M
rtbecomes less than a
predetermined small number or becomes monotonic.
The original signal x (t) can thus be expressed as following:

1
() () ()
M
jM
j

0.03
0.04

Fig. 2. A classical IMF
The major disadvantage of EMD is the so-called mode mixing effect. For example, the
simulated signal is defined as follows:

() sin(2 ) 10 ()* ( ) ( , 2, 1,0,1,2, )
0.2 0.015 , 0.2 0.03 0.215 0.03
()
0.215 0.015 , 0.215 0.03 0.23 0.03
0,1,2,3
st t wt t n n
tm mt m
wt
mt mt m
m



    

 




    


1
()
n
j
n
j
Xt c r




(4)
3.
Repeat step 1 and step 2 N times, but with different white noise serried w
i
(t) each time,
so

1
()
n
iijin
j
Xt c r



(5)

xz



 (7)
The vector
x
represents noisy signal and

is an unknown original clean signal.
z
is
independent identity distribution Gaussian white noise with mean zero and unit variance .
For simplicity, we assume intensity of noise is one. The step of wavelet shrinkage is defined
as follows:
1.
Apply discrete wavelet transform to observed noisy signal.
2.
Estimate noise and threshold value, thresholding the wavelet coefficients of observed
signal.
3.
Apply the inverse discrete wavelet transform to reconstruct the signal.
The wavelet shrinkage method relies on the basic idea that the energy of signal will often
be concentrated in a few coefficients in wavelet domain while the energy of noise is
spread among all coefficients in wavelet domain. Therefore, the nonlinear shrinkage
function in wavelet domain will tend to keep a few larger coefficients over threshold
value that represent signal, while noise coefficients down threshold value will tend to
reduce to zero.
In the wavelet shrinkage, how to select the threshold function and how to select the
threshold value are most crucial. Donoho introduced two kinds of thresholding functions:







 





(9)

Adaptive Noise Removal of ECG Signal Based On Ensemble Empirical Mode Decomposition

129
Hard threshold function (8) results in larger variance and can be unstable because of
discontinuous function. Soft threshold function (9) results in unnecessary bias due to
shrinkage the large coefficients to zero. We build the generalized threshold function:

1
()
m
m
m
xx
x



m
xxxIx Ix signx
x




    (12)
When m=1, it is soft threshold function; when m=

, it is hard threshold function. When
m=2 it is Non-Negative Garrote threshold function. We show slope signal as an example,
Figure.6 graphically shows generalized threshold functions for different m. It can be clearly
seen that when the coefficient is small, the smaller m is, the closer the generalized function is
to the soft threshold function; when the coefficient is big, the bigger m is, the closer the
generalized function is to the hard threshold function. As 1 m

, generalized threshold
function achieves a compromise between hard and soft threshold function. With careful
selection of m, we can achieve better denoising performance. Fig. 6. Generalized threshold function
We derived the exact formula of mean, bias, variance and
2
l risk for generalized threshold
function.

Adaptive Filtering Applications


 




and  are density and probability function of standard Gaussian random variable
respectively. Then:
Mean:

1
(,) (,) ()
mHm
m
MM A



 (13)
Bias:

2
(,) ( (,) )
mm
SB M
  
 (14)
Variance:
22 2
2122 1
( , ) ( , ) 2 ( ) ( ) ( ) 2 ( , ) ( )

( , ) ( 1)(2 ( ) ( )] ( )( ) ( )( ) ( , )
H H
V M

      
   

2
()1(1)(()( ))()()()()
H


     

  
(,)
m
M


,
(,)
m
SB


,
(,)
m
V

variance is reduced but the bias is increased. The optimal threshold value is the best
compromise between variance and bias and it should minimize the risk of the results as
compared with noise-free data.
Several methods have been proposed for the determinations of threshold values. The
universal threshold, proposed by Donoho and Johnstone, uses the fixed form threshold
equal to the square root of two times the logarithm of the length of the signal. LDT, the level
dependent threshold, proposed by I.M.Johnstone, and B.W.Silverman, uses a different
threshold for each of the levels based on a single formula. Stein Unbiased Risk Estimate
(SURE) is an adaptive threshold selection rule. It is data driven and the threshold value
minimizes an estimate of the risk. Other threshold values include minimaxi threshold etc. In
this paper, an adaptive threshold method is proposed based on EEMD. The threshold values
directly relate to the energy of noise on each IMFs. Next, the derivation of adaptive
threshold values is initiated by the characteristic of Fractional Gaussian noise (fGn).
fGn is a generalization of white noise. The statistical properties of fGn are controlled by a
single parameter H, and the autocorrelation sequence

,,
() ( )
HHiHik
rk EX X


(17)
This can also be defined as:

2
22 2
[] ( 1 2 1 )
2
HH H

In the decomposing of a given fGn, EMD is worked as a dyadic filter. Restricting to the
band-pass IMFs, self-similarity would mean that

(' )
'
', ,
() ( ) ' 2
kk
kk
kH kH H
H
Sf S f kk





 (20)
Given the self-similar relation (6) for PSDs for band-pass IMFs we can deduce how the
variance should evolve as a function of k:

(1)(')
['] [] ' 2
kk
HH
H
Vk Vk k k




H
n
Wcn



(23)
1
c represents the first IMF coefficients.
According to (21)

2( 1)
[] 2
Hk
HH
H
Vk C k




(24)
ˆ
[1]/
HH H
CW

 ，the parameter H and
H


(25)
For white noise,

1
2
H

, 0.719
H

 (26)

2
11
[2.01 0.2( ) 0.12( ) ] 2.01
22
H
HH

  
(27)
The energies of each IMFs can be defined as:

2
2.01 , 2,3,4
0.719
k
n
k
Vk

k
T from
(29); shrink the coefficients using the Non-Negative Garrote threshold function.
4.
Reconstruct the signal by the shrunken IMFs, obtain the denoised signal.

Adaptive Noise Removal of ECG Signal Based On Ensemble Empirical Mode Decomposition

133
Fig. 7. The block diagram of the denoising algorithm
6. Results and discussions
To evaluate the performance of the algorithm, Test signal and Clinic noisy ECG signals are
processed.
6.1 Test signal
We choose time shifted sine signal which shapes similarly to ECG to test above method;
Gaussian White Noise is added as noise, which is zero mean and standard deviation change
with the SNR.
10lo
g
(var( ) / var( ))SNR signal noise , var means standard deviation. The SNR
of noisy test signals are 5. Figure8 shows the original clean signal; figure9 shows the noisy
signal; figure10 shows the denoised signal by the above algorithm, the SNR of which
achieve 14. Furthermore, the original characteristic waveform is preserved. Fig. 8. The clean time shifted sine signal

Fig. 10. The denoised time shifted sine signal
6.2 Clinical noisy ECG signal
The ECG signal as Figure.11 illustrates comes from clinical patient. Signal is sampled at 360
Hz; signal length is 1500; the ECG signal is corrupted by noise. Figure12 shows its phase
space diagram, which is a plot of the time derivative of the ECG signal against the ECG
signal itself. The derivative can accentuate the noisy and high frequency content in signal, so
it can better show dramatic improvement after denoising. The noisy ECG signal is processed
using the method mentioned above. For the generalized threshold function, m is selected as
2, which is Non-Negative Garrote threshold function. The noisy ECG signal is decomposed
into a series of IMFs by EEMD. The first seven IMFs are shown in figure13; the latter seven
IMFs are shown in figure 14.The First IMF is discarded owing to predominant noise. Obtain
the adaptive threshold value of each IMFs by formula (29). The values are 0.0422,
0.0297,
0.0210,
0.0148, 0.0104, 0.0074, 0.0052, 0.0037, 0.0026, 0.0018, 0.0013, 0.0009, 0.0006. Then
shrink the coefficients of each IMFs by the adaptive threshold values and Non-Negative

Adaptive Noise Removal of ECG Signal Based On Ensemble Empirical Mode Decomposition

135
Garrote threshold function. The first shrunken six IMFs are shown in figure15; the latter
shrunken seven IMFs are shown in figure16. Reconstruct the signal by the shrunken 13 IMFs
and obtain the denoised signal. The filtered ECG signal is illustrated as figure17. The phase
space diagram of filtered ECG signal is shown as figure 18. From visual inspection, the ECG
signal is much cleaner after being denoised; the original characteristic waveform, especially
the sharp Q, R, and S peaks is preserved, without distorting the P and T waves.The results
indicate that the method we have proposed significantly reduces noise and well preserves
the characteristics of ECG signal.