#  Identifying the Mislabeled Training Samples of ECG Signals using Machine
Learning

Yaoguang Li, Wei Cui and Cong Wang This work is supported by the National
Natural Science Foundation of China under Grant 11404113, and the Guangzhou
Key Laboratory of Brain Computer Interaction and Applications under Grant
201509010006. Corresponding author: Wei Cui. The authors are with the School
of Automation Science and Engineering, South China University of Technology,
Guangzhou 510641, China (e-mails<EMAIL_ADDRESS><EMAIL_ADDRESS>wangcong@scut.edu.cn).

###### Abstract

The classification accuracy of electrocardiogram signal is often affected by
diverse factors in which mislabeled training samples issue is one of the most
influential problems. In order to mitigate this negative effect, the method of
cross validation is introduced to identify the mislabeled samples. The method
utilizes the cooperative advantages of different classifiers to act as a
filter for the training samples. The filter removes the mislabeled training
samples and retains the correctly labeled ones with the help of 10-fold cross
validation. Consequently, a new training set is provided to the final
classifiers to acquire higher classification accuracies. Finally, we
numerically show the effectiveness of the proposed method with the MIT-BIH
arrhythmia database.

###### Index Terms:

ECG signal, mislabeled samples, cross validation, machine learning.

##  I Introduction

Cardiac disease has become a key issue that threatens human life safety. How
to detect heart disease as early as possible is the core of this problem,
because cardiac disease is very prone to sudden death, which makes the
detection of heart disease more urgent than treatment. Electrocardiogram (ECG)
signal contains a lot of useful information that can be used to diagnose
various cardiac diseases. Nowadays, ECG is the most effective tool for heart
disease detection and the use of machine learning algorithms for automatic
detection of ECG signal has become an increasingly significant topic in the
relevant areas [1, 2, 3, 4]. Moreover, body sensor network-based devices have
become widely accepted and novel ECG telemetry systems for cardiac health
monitoring applications have been proposed [5, 6, 7]. In order to improve the
classification accuracy of ECG signal, most works focus on two aspects, (1)
feature selection [8, 9], and (2) robustness of the machine learning
classifiers [9, 10]. It is clear that the use of ECG signal features that
maximize the distinction between different diseases can significantly improve
the classification accuracy such as temporal intervals [8, 9, 10],
morphological features [8], frequency domain features [11], high-order
statistics [12], wavelet transform coefficients [9, 10, 13]. There are also
some studies on discriminant function optimization of different classification
algorithms [9, 14]. Ref.[8] shows that linear models achieve good
classification accuracies, nevertheless, more interests have been attracted to
the nonlinear approaches in the past few years. Neural network [10, 11, 15,
16] and support vector machine (SVM) [12, 13, 17] are the most popular
algorithms. In addition, some optimization algorithms such as genetic
algorithm [18] and particle swarm algorithm [19, 20] are used to optimize the
parameters in the classification algorithms to improve the classification
accuracy.

Although these methods mentioned above have been proved to be effective in
previous works, they are based on an essential assumption that the samples
used to train the classifier are completely reliable, which can not always be
guaranteed in the real world. Many factors may make ECG data less reliable
such as medical expert diagnosis error, data encoding and communication
problems. The unreliable issues are distinguished into two types: feature (or
attribute) noise and class (label) noise. In particular, if the training set
has mislabeled samples, the classifier will be deteriorated heavily, and
remarkably reduce the actual classification accuracy. Therefore, the automatic
analysis of ECG signal with computer technology is still an auxiliary
equipment in the cardiac disease diagnosis. In Ref.[21], the mislabeled
samples issue is the most prominent factor that reduces the classification
accuracy of ECG signal. In Ref.[22], it is summarized that with mislabeled
samples, the prediction accuracy decreases quickly and the complexity of
classifiers increases rapidly. For instance, Ref.[22] and Ref.[23] demonstrate
that the model of decision tree (DT) would become more complicated with label
noises. In Fig. 1 and Fig. 2, we show the classification accuracies with
different label noises levels in Ref.[24] and this manuscript, respectively.
Obviously, even the four commonly used machine learning algorithms are
adopted, the classification accuracies of ECG signal are still deteriorated
badly due to the presence of mislabeled samples in the training set.

Figure 1: Classification accuracy with label noise in Ref.[24]
Figure 2: Classification accuracy with label noise in this paper

Ref.[25] takes advantage of the Gaussian mixture model to model normal ECG
heartbeat and develops a real-time abnormal ECG heartbeat detection and
removal system. Ref.[24] proposes an optimal subset search algorithm based on
the genetic algorithm to achieve the highest classification accuracy in the
test set. It is considered that the training samples outside of the optimal
subset are regarded as the mislabeled, which need to be removed from the
original training set. Even if some of the mislabeled samples are identified
correctly by this method, the identification rate is not high enough and the
classification result is still unsatisfactory. Besides, the highest label
noise level in the original training set considered in Ref.[24] is only 20%.

The purpose of this manuscript is to enhance the identification rate of the
mislabeled samples and improve the classification accuracy of ECG signal with
the cross validation method. Compared with the previous studies, the proposed
method can achieve higher identification rate of mislabeled samples. Moreover,
if the mislabeled samples existed in the training set is less than 20%, the
classification accuracy can be increased to the same level as there is no
mislabeled sample in the training set with the help of the proposed method. If
the proportion of mislabeled samples is 30%, the classification accuracy is
slightly lower than the case of no mislabeled sample in the training set. If
the proportion meets 40%, the classification accuracy is still much higher
than the circumstance without filtering.

The rest of the paper is organized as follows. Section 2 introduces the cross
validation method and machine learning algorithms briefly. In Section 3, We
summarize the basic procedure of ECG signal processing, including data
collection, data preprocessing and feature extraction etc. Section 4 presents
the experiment results with the MIT-BIH arrhythmia database [26]. Finally, We
conclude our results in Section 5.

##  II Cross validation and Machine learning algorithms

###  II-A Cross validation

Cross validation [27, 28] is a statistical analysis method used to verify the
performance of a classifier. Its basic idea is to divide the original data
into two parts. One part serves as the training set, and the other serves as
the validation set. The training set is used to train the classifier, and then
the classifier is tested on the validation set to evaluate its performance.
There are four main cross validation methods.

a. Hold-out method: The original data is randomly divided into two groups, one
as the training set, and the other as the validation set. The training set is
used to train the classifier, and then the classifier is tested on the
validation set to evaluate its performance of classification accuracy. Because
the way to divide the original data into groups by the hold-out method is
random, the final classification accuracy on the validation set has a direct
relationship to the result of grouping, which makes the performance evaluation
of the hold-out method not convincing enough.

b. Double cross validation (2-CV): 2-CV divides the data into two equal-sized
subsets and performs two rounds of classifier training. In the first round, a
subset serves as the training set and the other is used for the validation
set. In the second round, the training set and the validation set are
interchanged to train the classifier again. The two classification accuracies
represent the performance of the classifier. But the 2-CV is not commonly
used, because the number of samples in the training set is so small that the
training set can not represent the distribution of all the data.

c. K-fold cross validation (K-CV): The original data is divided into $k$
groups (usually equalized). Each group serves as the validation set for one
time with the remaining $k-1$ groups as the training set. Then we obtain $k$
models and the corresponding $k$ classification accuracies as the performance
of the classifier. The number $k$ is usually larger than three. The K-CV can
effectively avoid over-learning and under-learning and obtain persuasive
results.

d. Leave-one-out cross validation (LOO-CV): If there are $N$ samples in the
original data, then the LOO-CV is the same as N-CV, which means each sample
serves as the validation set and the rest $N-1$ samples are the training set.
Then we obtain $N$ models and the average of the classification accuracies of
the $N$ models serves as the performance of the classifier. Compared with the
K-CV and the 2-CV, in the LOO-CV, almost all the samples of the original data
are used to train the model. So the training set of the LOO-CV can represent
the distribution of all the original data and the results of the performance
evaluation are more reliable. Moreover, there is no random factor to make the
experiment unrepeatable in the LOO-CV. But the high computational cost makes
the LOO-CV unworkable when the data size is large.

In this paper, 10-CV is adopted in view of the maximization of the use of the
original training samples and the minimization of computational consumption.

###  II-B Machine learning algorithm

####  II-B1 Support vector machine (SVM)

SVM [29, 30, 31] maps the feature vector $x\in{R^{d}}$ to the high-dimensional
feature space $\phi(x)\in H$, and creates an optimal separation hyperplane
that maximizes the interval between support vectors and the hyperplane.
Different SVM classifiers are generated by different mapping rules. The
mapping function $\phi(\cdot)$ is determined by the kernel function
$K(x_{i},x_{j})$, which defines an inner product of $H$ space.

The optimization problem of the SVM interval can be written as follows:

$\displaystyle\left\\{\begin{aligned}
&\max\limits_{\alpha}~{}~{}\sum_{i=1}^{N}{\alpha_{i}}-\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}\cdot
y_{i}y_{j}\cdot K(x_{i},x_{j})\\\
&~{}\text{s.t.}~{}~{}~{}~{}\sum_{i=1}^{N}{\alpha_{i}y_{i}}=0\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}0\leqslant\alpha_{i}\leqslant
C~{}~{}~{}~{}~{}i=1,2,\dots,N,\end{aligned}\right.$ (1)

and the SVM decision function is

$\displaystyle f(x)=\text{sgn}\left(\sum_{i=1}^{N}y_{i}\alpha_{i}\cdot
K(x,x_{i})+b\right).$ (2)

####  II-B2 K-nearest neighbor (KNN)

In KNN classification [32, 33, 34, 35], the output is a class membership. An
object is classified by a majority vote of its $k$ nearest neighbors ($k$ is a
positive integer, typically small). If $k=1$, the object is simply assigned to
the class of its single nearest neighbor.

####  II-B3 Naive Bayes (NB)

The bayes classification method [36, 37] classifies a certain sample based on
its posterior probability calculated by priori probability and data.

In bayes classifier, if D is a certain sample with the feature
$x=(x_{1},x_{2},\dots,x_{n})$, its probability of belonging to label $y_{i}$
is: $P(Y=y_{i}|X_{1}=x_{1},X_{2}=x_{2},\dots,X_{n}=x_{n}),(i=1,2,\dots,m)$,
that is

$\displaystyle P(Y$ $\displaystyle=y_{j}|X=x)=\max\\{P(Y=y_{1}|X=x),$ (3)
$\displaystyle P(Y=y_{2}|X=x),\dots,P(Y=y_{m}|X=x)\\}.$

Based on the bayesian formula:

$\displaystyle P(Y=y_{j}|X=x)=\frac{P(X=x|Y=y_{j})P(Y=y_{j})}{P(X=x)},$ (4)

the key point to obtain the posterior probability is to calculate
$P(X=x|Y=y_{j})$.

It is difficult to calculate $P(X=x|Y=y_{j})$ directly. So NB is created with
the assumption that the features of the object are independent of each other,
that is $P(X=x|Y=y_{j})=P(x_{1}|y_{j})P(x_{2}|y_{j})\cdots P(x_{n}|y_{j})$.
Consequently, the posterior probability of each class can be calculated and
the class with the largest posterior probability is just the classification
label of the sample.

####  II-B4 Linear discriminant analysis (LDA)

The principle of the LDA [38, 39, 40, 41] is that the data points (vectors)
are projected to the space of lower dimension so that the projected vectors
can be easily classified. Thus, the points that belong to the same type are
much closer than those belong to different types in the projection space.
Suppose that the projection function is $y=w^{T}x$, we can calculate

a. the original center point (mean value) of the class $i$ is

$m_{i}=\frac{i}{n_{i}}\sum_{x\in D_{i}}x,$ (5)

where $D_{i}$ refers to the set of points that belong to class $i$.

b. The center of class $i$ after projection is

$\widetilde{m_{i}}=w^{T}m_{i}.$ (6)

c. The degree of dispersion (variance) between classes is

$\widetilde{S_{i}}=\sum_{y\in{Y_{i}}}{(y-\widetilde{m_{i}})^{2}}.$ (7)

d. The objective optimization function of LDA after projection to $w$ is

$\displaystyle
J(w)=\frac{|\widetilde{m_{1}}-\widetilde{m_{2}}|^{2}}{{\widetilde{S_{1}}}^{2}+\widetilde{S_{2}}^{2}}.$
(8)

####  II-B5 Decision tree (DT)

A DT is a flowchart-like structure in which each internal node represents a
“test” on a feature, each branch represents the outcome of the test, and each
leaf node represents a class label (decision taken after computing all
attributes). The paths from root to leaf represent classification rules.

C4.5 is an algorithm used to generate a DT [42, 43, 44]. C4.5 is an extension
of the earlier ID3 algorithm [45]. The DTs generated by C4.5 can be used for
classification. Both C4.5 and ID3 use the information entropy to build DTs
from a set of training data. C4.5 makes some improvements to ID3. Some of
these are handling both continuous and discrete attributes, handling training
data with missing attribute values, handling attributes with differing costs,
and pruning trees after creation.

##  III ECG signal processing

ECG is a technique for recording the changes in electrical activity of each
cardiac cycle of the heart, which provides diagnostic information about the
cardiac condition of a patient.

As shown in Fig. 3, in a normal ECG cycle, P wave is the first upward
deflection and represents the electric potential variation of the two atrial
depolarization. It has positive polarity and its duration is between 80 ms and
110 ms. The QRS complex wave, with the duration of 60 ms–100 ms, is the most
important part in the ECG signal automatic analysis, and represents the
electric potential variation of the two ventricular depolarization. Q wave is
the first negative waveform, R wave is always the first positive deflection
that follows the Q wave, and S wave is the downward deflection after the R
wave [46]. T wave represents the electric potential variation of the two
ventricular repolarization, and its direction is the same as the QRS main
wave. We can detect various cardiac abnormalities with the help of the
difference in amplitude and duration of different waves. The most
distinguishable features of ECG signal are P wave period and peak, QRS-complex
period, R-R interval, R peak, etc. Doctors observe the deviations of P, QRS
and T waves from the normal signal in terms of interval and amplitude to find
the abnormalities in heart.

Since the computer performs faster and more accurate than the human eyes in
observing these features. Thus, it is obvious that the use of computers for
ECG signal analysis is a better way in cardiac detection and monitoring.
Normal sinus rhythm is the cardiac rhythm that begins with the myocardial
contraction of the sinus node. A normal sinus rhythm shows the following five
characteristics [46].

a. normal P wave,

b. constant P-P and R-R interval,

c. constant P wave configuration in a given lead,

d. P-R interval and QRS interval within normal limit,

e. heart rate between 60 to 100 beats/min.

It is divided into 11 different arrhythmias types in the MIT-BIH arrhythmia
database based on the different situations of abnormality or disturbance [48]:
left bundle branch block beat, right bundle branch block beat, aberrated
atrial premature beat, premature ventricular contraction, fusion of
ventricular and normal beat, nodal premature beat, atrial premature beat,
premature or ectopic supraventricular beat, ventricular escape beat, nodal
escape beat ,and paced beat. There are 12 categories in total with the normal
beat.

Fig. 4 describes several common arrhythmia ECG waveforms. From top to bottom
are normal beat, atrial premature beat, premature ventricular contraction,
left bundle branch block beat and right bundle branch block beat,
respectively. Atrial premature beats are characterized by an abnormally shaped
P wave. Since the premature beat initiates outside the sinoatrial node, the
associated P wave appears different from those seen in normal sinus rhythm.
Typically, the atrial impulse propagates normally through the atrioventricular
node and into the cardiac ventricles, resulting in a normal, narrow QRS
complex. In general, premature ventricular contraction has early appearance of
QRS complex, with no previous P wave. The QRS duration is more than 0.12
seconds, with large wave deformity. The ST segment and T wave have opposite
direction to QRS main wave with a complete compensation intermittent. Bundle
branch block beats’ QRS duration is more than 0.12 seconds with possible
double R waves (R-r’). Sometimes, we can only find a notch gap between R and
r’ in the left bundle branch block beat in lead V5, V6. However, the QRS wave
of the right bundle branch block beat is like a ’M’ shape. Many ECG waveforms
of different arrhythmias have similarities, which leads to the possibility of
mislabeling especially with the presence of noise.

Figure 3: A real time ECG signal [47]
Figure 4: Several common arrhythmia ECG waveforms

###  III-A ECG signal preprocessing

The amplitude of a normal ECG signal ranges from $10\mu V$ to $4mV$, thus the
ECG signal acquisition and processing should be considered as weak signal
detection issues. Based on the properties of ECG signal, an important step is
to filter out the general noises in the signal and to retain the valuable
components of each waveform. Two main noises in the ECG signal are baseline
wander and power interference.

Baseline wander has the strongest negative impact on ECG signal. It is a low-
frequency signal ranging from 0.05Hz to several Hz. We filter out the baseline
wander effectively using the median filtering [49] in this paper. The formula
of the median filtering is

$\displaystyle
y(i)=\text{Med}\bigg{\\{}x(i-N),\cdots,x(i),\cdots,x(i+N)\bigg{\\}},$ (9)

where the $x(i)$ is the original value, $y(i)$ is the new one, and $N$ is
adjustable.

Power interference mainly refers to 50Hz and its high harmonic interference.
The human body has antenna effects because of its natural physiological
characteristics, nevertheless the ECG signal collection equipment usually has
long wires exposed with antenna effects, which makes power interference the
most common noise in human ECG signal [50]. Due to its good time-frequency
localization characteristics, wavelet transform has become a popular method in
signal denoising. It decomposes the noisy signal into multi-scale, and then
the wavelet coefficients that belong to the power interference are removed
according to the frequency ranges of different kinds of signals. Ultimately,
the remaining wavelet coefficients are utilized to reconstruct the signal.

###  III-B ECG signal feature extraction

As described in the introduction part, there are different representation
types of ECG signal such as temporal intervals, frequency domain features,
high-order statistics, and so forth. In order to compare our results with the
previous work [24], we adopt the same features (1) ECG morphology features,
and (2) three ECG temporal features, i.e., the QRS complex duration, the RR
interval (the time span between two consecutive R points representing the
distance between the QRS peaks of the present and previous beats), and the RR
interval averaged over the ten last beats. The morphology features are
extracted from the segmented ECG cycles. We define the middle sample of two
R-peak samples as M, and a ECG cycle ranges from the previous M of the current
R-peak to the next M of the current R-peak. The morphology features are
acquired by normalizing the duration of the segmented ECG cycles to the same
periodic length according to the procedure reported in [51]. That is, one of
the ECG segments $y_{i}$=$[y_{i}(1),y_{i}(2),\cdots,y_{i}(n^{*})]$ can be
converted into a segment $x_{i}$=$[x_{i}(1),x_{i}(2),\cdots,x_{i}(n)]$ that
holds the same signal morphology, but in different data length (i.e.,
$n^{*}\not=n$) using the following equation,

$\displaystyle
x_{i}(j)=y_{i}(j^{*})+(y_{i}(j^{*}+1)-y_{i}(j^{*}))(r_{j}-j^{*}),$ (10)

where $r_{j}=(j-1)(n^{*}-1)/(n-1)+1$, and $j^{*}$ is the integral part of
$r_{j}$. In this paper, we set $n=300$. Therefore, the various lengths of the
ECG segments will be compressed or extended into a set of ECG segments with
the same periodic length. Consequently, the total number of morphology and
temporal features is equal to 303.

###  III-C Feature normalization and dimensionality reduction

Temporal features and morphology features have different dimensions, which
damages the classification results seriously. Data normalization means
adjusting values measured on different scales to a common scale. In this
paper, we use the min-max normalization rule [52] to map the original data to
0–1 by linear transformation

$\displaystyle x^{*}=\frac{x-min}{max-min},$ (11)

where $max$ and $min$ are the maximum and minimum values of the feature
dimensions of $x$, respectively.

Finally, due to the high-dimension properties of the features, we use the
principal component analysis (PCA) technique [53, 54, 15] to project the
features into a lower dimensional feature space.

##  IV Experiment

###  IV-A Data description

The method proposed for the mislabeled samples identification is tested
experimentally on real ECG signals that obtained from the well-known MIT-BIH
arrhythmia database. In particular, in order to facilitate the comparison with
Ref.[24], the considered beats referred to the following six classes: normal
sinus rhythm (N), atrial premature beat (A), ventricular premature beat (V),
right bundle branch block (RB), paced beat (P), and left bundle branch block
(LB). According to Ref. [24], the beats were selected from the recordings of
20 patients, which corresponded to the following files: 100, 102, 104, 105,
106, 107, 118, 119, 200, 201, 202, 203, 205, 208, 209, 212, 213, 214, 215, and
217. With 36328 heart beats in total, there are 24150 N class, 338 A class,
2900 V class, 3689 RB class, 3450 P class and 1801 LB class. All beats are
divided into the training set and the test set by the class distribution.
There are 1500 N class, 100 A class, 1000 V class, 1000 RB class, 1000 P
class, 500 LB class, totally 5100 beats in the training set. The rest of the
beats are assigned to the test set.

After the data preprocessing, feature extraction, normalization and reduction,
we establish our experiment in the next subsection.

###  IV-B Experiment setup

We assume that the original ECG recording labels of the MIT-BIH arrhythmia
database are totally correct. We add some label noises with different levels
(5%, 10%, 20%, 30%, 40%) to the training set, namely, changing the label of
some training samples artificially. However, the test set remains unchanged.
That means all the operations are on the training set, and the test set is
only used to test the effectiveness of the operations on the training set. A
noise-free training set is used for comparison. The number of the changed
labels of each class is based on its proportion to the overall training set.
For example, when label noise is 5%, the number of changed labels of class N
equals 1500*5%=75.

As shown in Fig. 5 and Fig. 6, with noise-free training set, we utilize the
cooperative advantages of different classifiers to act as a filter for the
training samples. The filter removes the mislabeled training samples and
retains the correctly labeled ones with a 10-fold cross validation.
Consequently, a new training set is fed to the final classifiers to acquire
higher classification accuracy [21]. Then we create classifiers based on three
machine learning algorithms, the naive bayes, the $k$ nearest neighbor and the
linear discriminant analysis. These classifiers are used to verify the test
set and calculate the classification accuracy. Afterwards, with the help of
the cross validation, we dig out the mislabeled samples in the training set.
The detailed procedures are the following.

All the training samples are randomly divided into 10 folds. One fold serves
as the validation set and the remaining 9 folds are the sub-training set for
each time. Then, the sub-training set is fed to five machine learning
algorithms (SVM, C4.5, NB, KNN, LDA) to classify the validation set. The
results of the classification are compared with the original label of the
validation set samples. According to the comparison results, one can testify
whether a certain sample in the validation set is a mislabeled sample or not.
The above process is repeated 10 times. Each time the selected validation set
is different so that all the samples in the training set can be verified. When
all the mislabeled samples in the training set have been identified and
removed, we get a new training set, retrain classifiers, and finally obtain
satisfactory classification accuracies on the test set.

In the above procedure, different machine learning classifiers may give
different validation results. So, we design three criterions to solve this
problem

a. Standard (1): If all five classifiers determine a single sample mislabeled,
we regard it as truly mislabeled.

b. Standard (2): If four or more classifiers determine a single sample
mislabeled, we regard it as truly mislabeled.

c. Standard (3): If three or more classifiers determine a sample single
mislabeled, we regard it as truly mislabeled.

Due to the contingency of a single experiment and its possible errors, in
order to enhance the reliability of the experiment and improve the
persuasiveness of this method, the experimental process should be repeated
several times. The identification circumstances of mislabeled samples are
shown detailedly in Tables 1-3.

Figure 5: Illustration of the proposed automatic training sample validation
framework
Figure 6: Flow chart of the proposed automatic training sample validation
framework

Table 1
Detection performance of mislabeled samples with standard (1)
Noise level ANM INM AINM $P_{D}$ $P_{FA}$ 5% 255 267 227 89% 16 % 10% 510 463
424 83% 7.65% 20% 1020 700 675 66% 2.45% 30% 1530 829 796 52% 2.16% 40% 2040
640 586 29% 2.65%

Table 2
Detection performance of mislabeled samples with standard (2)
Noise level ANM INM AINM $P_{D}$ $P_{FA}$ 5% 255 334 239 93.73% 37% 10% 510
572 474 93% 19.22% 20% 1020 912 825 81% 8.53 % 30% 1530 1344 1216 79% 8.37 %
40% 2040 1435 1207 59% 11.18%

Table 3
Detection performance of mislabeled samples with standard (3)
Noise level ANM INM AINM $P_{D}$ $P_{FA}$ 5% 255 873 251 98.43% 242% 10% 510
1079 504 99% 112.75% 20% 1020 1346 891 87% 44.61% 30% 1530 2050 1401 92%
42.42% 40% 2040 2567 1640 80% 45.44%

In Tables 1-3, ANM refers to the number of actual mislabeled samples, INM
refers to the number of identified mislabeled samples, AINM refers to the
number of the mislabeled samples that are corrected identified. $P_{D}$ refers
to identification accuracy of mislabeled samples,

$\displaystyle P_{D}=\frac{AINM}{ANM}.$ (12)

$P_{FA}$ refers to identification error rate,

$\displaystyle P_{FA}=\frac{INM-AINM}{ANM}.$ (13)

We conclude that at the same noise level, $P_{D}$ obtains its lowest value in
standard (1) and the highest value in standard (3). Standard (1) is a cautious
standard, once one of the five classifiers determines a certain sample as
correctly labeled, we regard it as the truly correctly labeled. Standard (3)
is the most relaxed and gives more mislabeled samples. In contrast, $P_{FA}$
obtains its highest value in standard (1) and the lowest value in standard
(3). Obviously, there is a trade-off between $P_{D}$ and $P_{FA}$, and the two
objectives optimization is of particularly interest. In combination with these
two aspects, standard (2) is the best standard. In Ref.[24], when the label
noise are 5%, 10%, 20%, $P_{D}$ are 78.46%, 78.40%, 72.40%, respectively, and
$P_{FA}$ are 31.05%, 15.65%, 4.58%, respectively. In this paper, with standard
(2), $P_{D}$ are 93.73%, 93%, 81%, respectively, and $P_{FA}$ are 37%, 19.22%,
8.53%, respectively. In summary, $P_{D}$ with the proposed method is
significantly higher than Ref.[24], with the cost of slightly higher $P_{FA}$.

Besides, with different standards in the procedure, we can identify the
corresponding mislabeled samples and obtain new training sets. Meanwhile, new
classifiers are updated to improve the classification accuracies. The details
are shown in Tables 4-6:

Table 4
Classification accuracy obtained by naive bayes classifier on the ECG beats
for all the considered training set scenarios
Naive bayes Noise level NF IF S1 S2 S3 0 73.30% – – – – 5% 71.96% 73.13%
74.67% 75.09% 75.48% 10% 69.33% 73.20% 74.16% 75.53% 75.23% 20% 62.01% 73.04%
71.60% 74.15% 73.87% 30% 51.85% 73.06% 64.95% 72.34% 70.71% 40% 38.77% 72.74%
49.68% 57.87% 54.79%

Table 5
Classification accuracy obtained by KNN classifier on the ECG beats for all
the considered training set scenarios
KNN Noise level NF IF S1 S2 S3 0 97.50% – – – – 5% 96.60% 97.29% 97.00% 97.25%
95.22% 10% 94.97% 97.20% 96.91% 96.80% 94.38% 20% 87.84% 97.12% 93.97% 95.53%
94.73% 30% 77.48% 96.96% 86.69% 92.22% 88.81% 40% 64.16% 96.86% 71.20% 78.15%
71.12%

Table 6
Classification accuracy obtained by LDA classifier on the ECG beats for all
the considered training set scenarios
LDA Noise level NF IF S1 S2 S3 0 74.50% – – – – 5% 73.88% 74.27% 74.28% 74.38%
73.93% 10% 73.41% 74.35% 74.12% 74.23% 73.35% 20% 71.27% 74.62% 73.27% 73.82%
72.49% 30% 63.78% 74.35% 68.87% 70.95% 67.82% 40% 48.89% 74.40% 55.16% 58.27%
50.12%

In Tables 4–6 and Figs. LABEL:NB,KNN,LDA , NF refers to no filtering. IF means
ideal filtering (artificially removes all the mislabeled samples from the
training set). S1, S2, S3 represent standard (1), standard (2) and standard
(3), respectively. Obviously, the classification accuracies of NB, KNN and LDA
reduce remarkably with the increase of the label noise level in the training
set without filtering. In the case of ideal filtering , the classification
accuracies are improved to almost the same level as the noise level equals 0.
With different standards, the proposed method gets different filtering
results. Specifically, the classification accuracies on the test set have been
significantly improved. In some individual cases, such as the label noise
equals 10%, the accuracy of the proposed method is even higher than that of
the ideal filtering, because our method can not only filters out the
artificially added label noise, but also removes some undiscovered noise.
These noises may come from feature extraction or other processes after ECG
denoising. From the experimental results in Tables 4–6, we note that standard
(2) has the best classification accuracy on the test set, which is consistent
with what we mentioned in Sec.IV(B). As shown in Figs. 6–8, if the mislabeled
samples existed in the training set is lower than 20%, the classification
accuracy can be increased to the same level as there is no mislabeled sample
in the training set with the help of the proposed method. If the proportion is
30%, the classification accuracy is slightly lower than that with no
mislabeled sample. If the proportion meets 40%, the classification accuracy is
still much higher than the circumstances without filtering.

Figure 7: NB classification accuracy
Figure 8: KNN classification accuracy
Figure 9: LDA classification accuracy

##  V Conclusion

In this paper, we investigate how to identify and eliminate the mislabeled
training samples that are widely existed in ECG analysis. We use the cross
validation method to improve the identification rate and the classification
accuracy of ECG signal. Simulation results demonstrate the effectiveness of
our proposed method. The mislabeled samples in the training set directly
damage classifiers and their classification accuracy seriously. We use five
different machine learning classifiers to classify the validation set, and to
improve the reliability of the mislabeled samples identification. Especially,
if the label noise level is not higher than 20%, the classification accuracy
can be improved to the same level as there is no mislabeled samples in the
training set. It is noted that, the mislabeled noises are randomly added
rather from the similarity of different arrhythmias in the actual ECG
waveform. There are also some points needing for improvement in our further
works. For example, the computational load is relatively expensive due to the
large number of classifiers, and also, our method can not be competent for too
high noise levels. Further researches are necessary to address these problems
and to make this method more practical and effective.

## References

  * [1] A. Batra and V. Jawa, “Classification of arrhythmia using conjunction of machine learning algorithms and ECG diagnostic criteria,” Int. J. Biol. Biomed., vol. 1, pp. 1-7, 2016.
  * [2] M. M. Al Rahhal, Y. Bazi, H. AlHichri, N. Alajlan, F. Melgani, and R. R. Yager, “Deep learning approach for active classification of electrocardiogram signals,” Inf. Sci. (Ny)., vol. 345, pp. 340-354, 2016.
  * [3] M. Huanhuan and Z. Yue, “Classification of electrocardiogram signals with deep belief networks,” Computational Science and Engineering (CSE), 2014 IEEE 17th International Conference on. pp. 7-12, 2014.
  * [4] T. P. Exarchos, C. Papaloukas, D. I. Fotiadis, and L. K. Michalis, “An association rule mining-based methodology for automated detection of ischemic ECG beats,” IEEE Trans. Biomed. Eng., vol. 53, no. 8, pp. 1531-1540, 2006.
  * [5] S. Y. Lee, J. H. Hong, C. H. Hsieh, M. C. Liang, S. Y. C. Chien, and K. H. Lin, “Low-power wireless ECG acquisition and classification system for body sensor networks,” IEEE Journal of Biomedical and Health Informatics, vol. 19, pp.236-246, 2015.
  * [6] W. Chen, S. B. Oetomo, D. Tetteroo, F. Versteegh, T. Mamagkaki, M. S. Pereira, L. Janssen, and A. van Meurs, “Mimo pillow–an intelligent cushion designed with maternal heart beat vibrations for comforting newborn infants,” IEEE Journal of Biomedical and Health Informatics, vol. 19, pp.979-985, 2015.
  * [7] L. Wang, G. Z. Yang, J. Huang, J. Zhang, L. Yu, Z. Nie, and D. R. S. Cumming, “A wireless biomedical signal interface system-on-chip for body sensor networks”, IEEE Transactions on Biomedical Circuits and Systems, vol. 4, pp. 112-117, 2010.
  * [8] P. de Chazal, M. O Dwyer, and R. Reilly, “Automatic classification of heartbeats using ECG morphology and heartbeat interval features,” IEEE Trans. Biomed. Eng., vol. 51, no. 7, pp. 1196-1206, 2004.
  * [9] O. T. Inan, L. Giovangrandi, and G. T. A. Kovacs, “Robust neural-networkbased classification of premature ventricular contractions using wavelet transform and timing interval features,” IEEE Trans. Biomed. Eng., vol. 53, no. 12, pp. 2507-2515, 2006.
  * [10] T. Ince, S. Kiranyaz, and M. Gabbouj, “A generic and robust system for automated patient-specific classification of ecg signals,” IEEE Trans. Biomed. Eng., vol. 56, no. 5, pp. 1415-1426, 2009.
  * [11] C. H. Lin, “Frequency-domain features for ECG beat discrimination using grey relational analysis based classifier,” Comput. Math. Appl., vol. 55, no. 4, pp. 680-690, 2008.
  * [12] S. Osowski, L. T. Hoai, and T. Markiewicz, “Support vector machine-based expert system for reliable heartbeat recognition,” IEEE Trans. Biomed. Eng., vol. 51, no. 4, pp. 582-589, 2004.
  * [13] A. Daamouche, L. Hamami, N. Alajlan and F. Melgani. “A wavelet optimization approach for ECG signal classification,” Biomedical Signal Processing and Control., vol. 7, no. 4, pp. 342-349, 2012.
  * [14] J. Fayn, “A classification tree approach for cardiac ischemia detection using spatiotemporal information from three standard ECG leads,” IEEE Trans. Biomed. Eng., vol. 58, no. 1, pp. 95-102, 2011.
  * [15] R. J. Martis, U. R. Acharya and L. C. Min. “ECG beat classification using PCA, LDA, ICA and Discrete Wavelet Transform,” Biomedical Signal Processing and Control., vol. 8, no. 5, pp. 437-448, 2013.
  * [16] W. Jiang and S. Kong, “Block-based neural networks for personalized ECG signal classification,” IEEE Trans. Biomed. Eng., vol. 18, no. 6, pp. 1750-1761, 2007.
  * [17] A. Kampouraki, G. Manis, and C. Nikou, “Heartbeat time series classification with support vector machines,” IEEE Trans. Inf. Technol. Biomed., vol. 13, no. 4, pp. 512-518, 2009.
  * [18] J. A. Nasiri, M. Naghibzadeh, H. S. Yazdi, and B. Naghibzadeh, “ECG arrhythmia classification with support vector machines and genetic algorithm,” in Proc. Comput. Model. Simul., pp. 187-192, 2009.
  * [19] F. Melgani and Y. Bazi, “Classification of electrocardiogram signals with support vector machines and particle swarm optimization,” IEEE Trans. Inf. Technol. Biomed., vol. 12, no. 5, pp. 667-677, 2008.
  * [20] M. Korürek and B. Dogan, “ECG beat classification using particle swarm optimization and radial basis function neural network,” Expert Syst. Appl., vol. 37, no. 12, pp. 7563-7569, 2010.
  * [21] C. E. Brodley and M. A. Friedl, “Identifying mislabeled training data,” J. Artif. Intell. Res., vol. 11, pp. 131-167, 1999.
  * [22] B. Frénay and M. Verleysen, “Classification in the presence of label noise: a survey,” IEEE Trans. Neural. Netw. Learn. Syst., vol. 25, no. 5, pp. 845-869, 2014.
  * [23] J. R. Quinlan, “Induction of decision trees,” Mach. Learn., vol. 1, no. 1, pp. 81-106, 1986.
  * [24] E. Pasolli and F. Melgni, “Genetic algorithm-based method for mitigating label noise issue in ECG signal classification,” Biomedical Signal Processing and Control., vol. 19, pp. 130-136, 2015.
  * [25] W. Louis, M. Komeili and D. Hatzinakos, “Real-time heartbeat outlier removal in electrocardiogram(ECG) biometrie system,” IEEE Canadian Conference on Electrical Computer Engineering., Vancouver, BC, May. 2016, pp. 1-4.
  * [26] R. Mark and G. Moody, MIT-BIH Arrhythmia Database.(1997).[Online]. Available: http://ecg.mit.edu/dbinfo.html.
  * [27] S. Arlot and A. Celisse, “A survey of cross-validation procedures for model selection,” Statistics Surveys., vol. 4, pp. 40-79, 2010.
  * [28] Y. Zhang and Y. Yang, “Cross-validation for selecting a model selection procedure,” J. Econometrics., vol. 187, no. 1, pp. 95-112, 2015.
  * [29] S. Theodoridis and K. Koutroumbas, “Pattern Recognition,” 4th ed. New York: Academic, 2008.
  * [30] N. Cristianini and J. Shawe-Taylor, “An Introduction to Support Vector Machines,” Cambridge Univ. Press, 2000
  * [31] J. A. K. Suykens, “Support vector machines: a nonlinear modelling and control perspective,” Euro. J. Control, Special Issue on Fundamental Issues in Control, vol. 7, no. 2-3, pp. 311-327, Aug. 2001.
  * [32] D. Coomans and D. L. Massart, “Alternative k-nearest neighbour rules in supervised pattern recognition: Part 1. k-nearest neighbour classification by using alternative voting rules,” Analytica Chimica Acta, vol. 136, pp. 15 C 27, 1982.
  * [33] T.M. Cover and P.E. Hart, Nearest Neighbor Pattern Classification, IEEE Trans.Information Theory, vol. 13,pp. 21-27,Jan.1967.
  * [34] M. L. Zhang and Z. H. Zhou, “ML-KNN: A lazy learning approach to multilabel learning,” Pattern Recognit., vol. 40, no. 7, pp. 2038-2048, 2007.
  * [35] J. Keller, M. Gray, and J. Givens, “A fuzzy k-nearest neighbor algorithm,” IEEE Trans. Syst. Man, Cybern., vol. SMC-15, pp. 580-585, 1985.
  * [36] S. Viaene, A. Derrig and G. Dedene, “A case study of applying boosting naive bayes to claim fraud diagnosis,” IEEE Trans. Knowl. Data Eng., vol. 16, no. 5, pp. 612-620, 2004.
  * [37] S. B. Kim, K. S. Han, H. C. Rim, and S. H. Myaeng, “Some effective techniques for naive bayes text classification,” IEEE Trans. Knowl. Data Eng., vol. 18, no. 11, pp. 1457-1466, 2006.
  * [38] R.A. Fisher, “The Use of Multiple Measures in Taxonomic Problems,” Ann. Eugenics, vol. 7, pp. 179-188, 1936.
  * [39] G. J. McLachlan, “Discriminant Analysis and Statistical Pattern Recognition.” Hoboken, NJ, USA: Wiley, 2004.
  * [40] J. Ye, R. Janardan, and Q. Li, “Two-dimensional linear discriminant analysis,” in Advances in Neural Information Processing Systems (NIPS). Cambridge, MA: MIT Press, 2004, pp. 1569-1576.
  * [41] R. Haeb-Umbach and H. Ney, “Linear discriminant analysis for improved large vocabulary spech recognition,” in Proc. ICASSP, vol. 1, pp. 13-16, 1992.
  * [42] J.R. Quinlan, “C4.5: Programs for Machine Learning.” San Mateo, Calif.: Morgan Kaufmann, 1993.
  * [43] J.R. Quinlan, “Improved Use of Continuous Attributes in C4.5,” J. Artificial Intelligence Research, vol. 4, pp. 77-90, 1996.
  * [44] S. Ruggieri, “Efficient c4.5,” IEEE Trans. Knowledge and Data Eng., vol. 14, no. 2, pp. 438-444, 2002.
  * [45] N. R. Pal and S. Chakraborty, “Fuzzy rule extraction from ID3-type decision trees for real data,” IEEE Trans. Syst, Man, Cybern. B., vol. 31, pp. 745-754, 2001.
  * [46] T. Tabassum, M. Islam, “An Approach of Cardiac Disease Prediction by Analyzing ECG Signal,” International Conference on Electrical Engineering and Information Communication Technology., Dhaka, Bangladesh, September. 2016, pp. 1-5.
  * [47] http://www.todayifoundout.com/index.php/2011/10/how-to-read-an-ekg-electrocardiograph/
  * [48] X. Qin, Y. Yan, J. Fan and L. Wang, “Real-time electrocardiogram analysis based on multi-kernel extreme learning machine,” Journal of Integration Technology., vol. 4, no. 5, pp. 36-45, 2015.
  * [49] G. Wischermann, “Median filtering of video signals-A powerful alternative,” SMPTE J., pp. 541-546, 1991.
  * [50] J. A. van Alste and T. S. Schilder, “Removal of base-line wander and power-line interference from the ECG by an efficient FIR filter with a reduced number of taps,” IEEE Trans. Biomed. Eng., vol. BME-32, pp. 1052-1060, 1985.
  * [51] J. J. Wei, C. J. Chang, N. K. Shou and G. J. Jan, “ECG data compression using truncated singular value decomposition,” IEEE Trans. Biomed. Eng., vol. 5, no. 4, pp. 290-299, 2001.
  * [52] Y. K. Jain and S. K. Bhandare, “Min max normalization based data perturbation method for privacy protection,” in International Journal of Computer and Communication Technology (IJCCT)., vol. 2, pp. 45-50, 2011.
  * [53] S. Wold, K. Esbensen, and P. Geladi, Principal component analysis, Chemometr. Intell. Lab. Syst., vol. 2, no. 1 C3, pp. 37 C52, Aug. 1987.
  * [54] I.T. Jolliffe, “Principal Component Analysis,” second ed. Springer-Verlag, 2002.


# ACOUSTIC Room COMPENSATION USING Local PCA-BASED ROOM AVERAGE Power Response
ESTIMATION

###### Abstract

Acoustic room compensation techniques, which allow a sound reproduction system
to counteract undesired alteration to the sound scene due to excessive room
resonances, have been extensively studied. Extensive efforts have been
reported to enlarge the region over which room equalization is effective and
to contrast variations of room transfer functions in space. A speaker-tuning
technology “Trueplay” allows users to compensate for undesired room effects
over an extended listening area based on a spatially averaged power response
of the room, which is conventionally measured using microphones on portable
devices when users move around the room. In this work, we propose a novel
system that leverages measured speaker echo path self-responses to predict the
room average power responses using a local PCA based approach. Experimental
results confirm the effectiveness of the proposed estimation method, which
further leads to a room compensation filter design that achieves a good sound
similarity compared to the reference system with the ground-truth room average
power response while outperforming other systems that do not use the proposed
estimator.

##  1 Introduction

When sound is reproduced by a sounding unit, e.g. a loudspeaker, the
perception of the sound that is delivered to the listener is modified by the
listening environment. Excessive reflections or resonances may result in
detrimental impacts on the original sound, e.g. altered music timbre,
incorrect sound localization. The sound reproduction system can also introduce
undesired artifacts, e.g. frequency band-extension, non-linearities and etc..

Room response compensation (RRC) has been widely studied [1, 2, 3, 4, 5, 6]
for improving sound reproduction quality, contrasting detrimental effects of
the room environment. In an RRC system, the room transfer function (RTF)
characterizing the path from the sound reproduction system to the listener is
equalized with a suitably designed equalizer through its inversion. From [7],
RRC approaches are categorized into single-point and multi-point room
equalizers. A single-point room equalizer derives the equalization filter
based on the measurement in a single location of the RTF [8], which is
effective only in a limited zone around the measured point. In real-world
scenarios, RTFs vary significantly with respect to the position in the room
[9] and time [10, 11]. To enlarge the equalization zone and to contrast the
room response variations, multi-point equalizers have been proposed [12]. A
multi-point room equalization system leverages multiple measurements of RTFs
at different locations in order to design the equalizer in a either fixed or
adaptive manner. However, the placement and measurements at numerous separate
locations that covers the extended room area could be a time-consuming
process.

Alternatively, RTFs can be estimated from a finite set of room impulse
response measurements at specific positions in the room through model-based
approaches. Several works have been reported in the literature, such as
common-acoustical pole and residue modelling [13], compressive sensing [14,
15] and the planewave decomposition method [16, 17]. However, depending on the
formulation of sparsity models that are used to characterize the entire
soundfield in a room, these approaches still require a relatively large number
of microphone measurements to achieve sufficient accuracy and microphone
positions in space are also required. The estimation of RTFs is often severely
under-determined, which is typically due to the large number of independent
parameters needed to represent the acoustic system in reverberant
environments. To this end, Fozunbal et al. [18] deployed Principal Component
Analysis (PCA) and proposed a supervised algorithm by forming a model based on
a training set that consists of pre-measured acoustic impulse responses from
several known locations. Koren et al. [19] further improved the PCA-based
estimation approach by exploiting spatial information of localized sources in
order to properly train local models that corresponds to different regions in
the room. PCA is also applied in [20] to facilitate the estimation of the
sound pressure at the eardrum using an inward-facing microphone at the in-ear
earpiece.

An efficient variant of the multi-point equalization system can be formulated
based on a spatially-averaged power response measurement. While a power
measurement necessarily eliminates phase information, such ”room curves” are
commonly used by speaker designers and reveal problems that are common to all
locations, such as boundary effects [21]. Sheen [22] proposed a speaker tuning
technology “Trueplay” [23] based on room power response measurements. The room
average power response is typically derived using microphones on portable
devices such as mobile phones. The user moves around the room with the
portable device while the speaker under test plays a cyclic stimulus signal.
The power response of each cycle is averaged to create the room average power
response.

In this paper we present a novel supervised approach to estimate the room
average power response by forming models based on a training set. The system
requires only an on-board microphone while eliminating the need for additional
microphone measurements and the knowledge of speaker position in space. The
proposed room compensation design uses the echo path, i.e. the acoustic
impulse response between the speaker and on-board microphone, to predict the
room average power response. Specifically, the echo path measurement is
projected onto a lower dimensional space using PCA and transformed to obtain
an estimate of the room average power response. It also exploits local
features of echo path self-responses (e.g. reverberation time) and divides the
global training dataset into local subsets, thereby deriving separate local
PCA models that corresponds to different room conditions. Results demonstrate
the improved performance of the local PCA approach over other conventional
estimation schemes and when this estimate is incorporated into the design of
“Trueplay” for room compensation, the proposed approach achieves good
similarity to the reference system with the ground-truth room average power
response.

##  2 Problem Formulation

Fig. 1: Loudspeaker room compensation acoustic setup.

The acoustic setup of “Trueplay” is depicted in Fig. 1. The room average power
response $r(\omega)$ can be measured at radial frequency $\omega$ using an
external microphone as described in Sec. 1. The on-board microphone measures
the echo path self-response $s(\omega)$ of the device. The goal is to derive
an equalization filter $G_{eq}(\omega)$ such that the room average power
response meets a target magnitude response $t(\omega)$:

$\vspace{-0.1cm}r(\omega)|G_{eq}(\omega)|^{2}=t(\omega).$ (1)

The target response is a subjective choice made by the loudspeaker designer.
Targets are commonly flat across frequency with some boosts in the bass
portion of the response. The magnitude response of the ideal filter is then
given by

$\vspace{-0.1cm}|G_{eq}(\omega)|=\sqrt{\frac{t(\omega)}{r(\omega)}}.\vspace{-0.1cm}$
(2)

The magnitude response is then turned into a minimum phase finite impulse
response (FIR) filter by taking it’s real cepstrum [24], providing a
realizable filter. The FIR can be used directly for equalization. To reduce
computational complexity, the Trueplay system instead computes a set of second
order infinite impulse response (IIR) filters that approximate the magnitude
and phase response of the ”ideal” FIR. Various techniques exist for the design
of arbitrary IIR filters [25][26].

The main challenge of this paper is to replace $r(\omega)$ in Eq. (2) with an
appropriate estimator of the room average power response $\hat{r}(\omega)$
based on the in-situ measurement of the echo path self-response, which will be
described in the following section.

##  3 Estimation of room average power response

In this section, we present an estimation scheme to map the measured echo path
self-response at the on-board microphone to the room average power response.
The first subsection introduces a linear least-sqaures (LS) regression method
that minimizes mean-squared error of the estimated room average power response
coefficients in the frequency domain. Subsequently, we propose an estimation
method that benefits from the numerical robustness and efficiency of PCA.
Finally, we motivate a local PCA based method that trains individual local
models corresponding to different room reverberant conditions.

###  3.1  Linear least-squares regression

Let
$\mathcal{M}=\\{\mathbf{s_{j}},\mathbf{r_{j}}\in\mathbb{R}^{\frac{N_{f}}{2}+1}|j=1,\ldots,J\\}$
be a set of magnitude spectrum coefficients of echo path self-responses
$s_{j}(\omega)$ and room average power responses $r_{j}(\omega)$ in the log-
energy form. $N_{f}$ is the FFT length. The optimal filter $\mathbf{g_{LS}}$
should minimize the difference between the estimated room average power
response to the measured counterpart via linear mapping in the frequency
domain, i.e., the following least-squares cost function

$E(\mathbf{g_{LS}})=\|\mathbf{D_{s}}\mathbf{g_{LS}}-\mathbf{d_{r}}\|^{2}_{2}+\mu\|\mathbf{g_{LS}}\|^{2}_{2},$
(3)

where $\mathbf{D_{s}}$ ($J(\frac{N_{f}}{2}+1)\times(\frac{N_{f}}{2}+1)$) are
stacked diagonal matrices containing magnitude coefficients of all measured
echo path self-responses $s_{j}(\omega)$ and $\mathbf{d_{r}}$
($J(\frac{N_{f}}{2}+1)\times 1$) is the stacked vectors containing all
$r_{j}(\omega)$. The Tikhonov regularization factor $\mu=0.001$ is considered
in this work to prevent over-amplification of $\mathbf{g_{LS}}$. The optimum
with respect to $\mathbf{g_{LS}}$ is given by

$\mathbf{\hat{g}_{LS}}=(\mathbf{D_{s}}^{H}\mathbf{D_{s}}+\mu\mathbf{I})^{-1}\mathbf{D_{s}}^{H}\mathbf{d_{r}}.$
(4)

Using the estimated filter $\mathbf{\hat{g}_{LS}}$ from the training stage,
room average power response coefficients $\mathbf{\hat{r}_{LS}}$ are estimated
from a measurement of the echo path self-response $\mathbf{s}$ during run-time
as follows:

$\mathbf{\hat{r}_{LS}}=\mathbf{s}\odot\mathbf{\hat{g}_{LS}},$ (5)

where $\odot$ denotes element-wise product.

###  3.2  Global PCA-based estimation

As can be seen from (3), a direct linear mapping of the frequency domain
vectors $\mathbf{s_{j}},\mathbf{r_{j}}$ is possible but would require a large
set of training data to ensure sufficient accuracy if the FFT length is large.
In this section, a PCA-based estimator of $\mathbf{r}$ is designed based on
features of measured echo path self-responses $\mathbf{s}$.

PCA is commonly used for dimensionality reduction by projecting each data
point onto only the first few principal components to obtain lower-dimensional
data that helps to avoid over-fitting [27]. By conducting PCA, we extract the
first $K_{s}$ and $K_{r}$ principal components
$\mathbf{u}_{s,k},\mathbf{u}_{r,k}\in\mathbb{R}^{\frac{N_{f}}{2}+1}$ of
$\mathbf{s_{j}}$ and $\mathbf{r_{j}}$, respectively. The room average power
response principal components matrix is defined as

$\mathbf{U}_{r}=[\mathbf{u}_{r,1},\mathbf{u}_{r,2},\ldots,\mathbf{u}_{r,K_{r}}].$
(6)

Let $\mathbf{\bar{r}}$ be the ensemble average of $\mathbf{r_{j}}$:
$\mathbf{\bar{r}}=\sum_{j=1}^{J}\mathbf{r_{j}}/J$. To obtain gain vectors
$\mathbf{g}_{r,j}$ that minimize the Euclidean distance between the
reconstructed frequency domain vectors

$\mathbf{\hat{r}}_{j}=\mathbf{\bar{r}}+\mathbf{U}_{r}\mathbf{g}_{r,j}$ (7)

and the raw $\mathbf{r}_{j}$, we utilize the orthonormality of the principal
components and get

$\mathbf{g}_{r,j}=\mathbf{U}_{r}^{H}(\mathbf{r}_{j}-\mathbf{\bar{r}}).$ (8)

Similarly we obtain the gains $\mathbf{g}_{s,j}$ for the echo path.

After converting frequency domain coefficients into the principal component
domain, the problem is to find a linear map
$\mathbf{A}\in\mathbb{R}^{K_{s}\times K_{r}}$ that projects the echo path
self-response gain vectors onto the room average power response gain vectors.
The following cost function is defined:

$\vspace{-0.2cm}C(\mathbf{A})=\sum_{j=1}^{J}\|\mathbf{g}_{r,j}-\mathbf{A}\mathbf{g}_{s,j}\|^{2}.\vspace{-0.1cm}$
(9)

This projection allows us to estimate the room average power response gain
vector given the echo path self-response gain vector. To minimize
$C(\mathbf{A})$, we have

$\vspace{-0.2cm}\mathbf{\hat{A}}=\underset{\mathbf{A}}{\mbox{argmin}}\
C(\mathbf{A})=\sum_{j=1}^{J}\mathbf{g}_{r,j}\mathbf{g}_{s,j}^{H}(\sum_{j=1}^{J}\mathbf{g}_{s,j}\mathbf{g}_{s,j}^{H})^{-1}.\vspace{-0.1cm}$
(10)

The above-demonstrated steps are the training stage for room average power
response estimation based on a training set
$\mathbf{s_{j}},\mathbf{r_{j}}\in\mathcal{M}$. After measuring the echo path
self-response $\mathbf{s}$ at runtime, the gain vector for $\mathbf{s}$ in the
the principal component domain can be calculated as follows:

$\mathbf{g}_{s}=\mathbf{U}_{s}^{H}(\mathbf{s}-\mathbf{\bar{s}}),$ (11)

where $\mathbf{\bar{s}}$ is the ensemble average of $\mathbf{s}_{j}$ from the
training stage. Then an estimate of $\mathbf{g}_{r}$ can be obtained by
$\mathbf{\hat{g}}_{r}=\mathbf{\hat{A}}\mathbf{g}_{s}$. Finally an estimate
$\mathbf{\hat{r}_{pca}}$ for $\mathbf{r}$ in the frequency domain with the
ensemble average $\mathbf{\bar{r}}$ and $\mathbf{U}_{r}$ from the training
stage

$\mathbf{\hat{r}_{pca}}=\mathbf{\bar{r}}+\mathbf{U}_{r}\mathbf{\hat{g}}_{r}.$
(12)

###  3.3  Local PCA-based estimation

The PCA-based estimator of the room average power response in Sec. 3.2 is
derived based on the global training dataset that consists of a mixture of
measured echo path self-response/room average power response sets under a
large variety of room conditions. In this section, we propose to exploit local
features of echo path self-responses to “cluster” the global training dataset
into small subsets and to train local PCA models that corresponds to different
room conditions accordingly, which facilitates PCA to better capture the
structure of acoustic responses and their variation patterns.

The following features that divide the global dataset into sub-groups are
investigated in the training stage:

  * •

Reverberation time RT30 [28] of the measured echo path self-response to
quantify the reverberant condition at the position where the speaker is
placed;

  * •

Low-frequency roll-off gain of the echo path frequency response. In this work,
we choose to monitor the magnitude roll-off between 120 Hz and 60 Hz.

At the inference stage, the selection criteria that indicates which local PCA
model to be used for deriving an estimated room average power response
$\mathbf{\hat{r}_{lpca}}$ when a new echo path self-response becomes available
is also established accordingly. The effectiveness of the proposed local PCA
estimation approach is evaluated in Sec. 4.

##  4 System Evaluation

In this section, we evaluate the effectiveness of incorporating the proposed
room average power response estimator into the ”Trueplay” room compensation
filter design. First the acoustic setup and the database of measurements are
introduced. Second, we present results that demonstrate the estimation
accuracy of the proposed PCA-based approaches. Finally, results of a
perceptual evaluation of the proposed room compensation design utilizing the
multi stimulus hidden reference and anchor (MUSHRA) framework are presented.

###  4.1 Self-response/Room average power response database

Fig. 2: Echo path self-responses and room average power responses from the
collected database. Bold lines denote ensemble average.

Fig. 3: The distribution of RT30 and low-frequency roll-off gains of echo path
self-responses for Sonos Move database.

For data collection of echo path self-response/room average power response
sets, a Sonos Move speaker (with built-in microphones) and one pre-calibrated
external microphone were used. Room average power responses were measured as
described in Sec. 1 and echo path self-responses were derived by playing sine
sweep signals from the speaker. The Sonos Move speaker was placed in different
rooms with a wide diversity of room types and at various locations in rooms
(e.g. center, corner, near walls and etc.). Overall, as shown in Fig. 2, the
database includes a total of 1207 echo path self-response/room average power
response sets at a sampling rate of 16 kHz and FFT length $N_{f}=2048$. The
distribution histograms of the two investigated local features of echo path
self-responses in Sec. 3.3, i.e. RT30 and low-frequency roll-off gains, are
demonstrated in Fig. 3. From Fig. 3, it can be seen that both features follow
a normal distribution well. Therefore, we fitted them with two normal
distribution models and used the derived mean and standard deviation values to
define the local subgroup division schemes for each of the features shown in
Table 1.

Table 1: Local subgroup division schemes | RT30 | Low-freq roll-off
---|---|---
Group 1 | [0,87) ms, 265 sets |  ($-\infty$,10) dB,211 sets
Group 2 | [87,151] ms, 679 sets | [10,18] dB, 758 sets
Group 3 |  (151,+$\infty$) ms, 263 sets |  (18,$\infty$) dB, 208 sets

When evaluating the performance we used a repetitive cross-validation
approach, i.e., we used data from $N_{tr}$ echo path self-response/Room
average power response sets for training then evaluated the performance for
$N_{val}$ unseen sets and repeated this process for 50 times. Note that the
training data size for the local PCA-based method was $N_{tr}/3$ for each of
the subgroup.

###  4.2 Estimation of room average power response

To evaluate the performance of the room average power response estimators, we
used the following measure to quantify the dB level error between the
estimated room average power response $\hat{r}(\omega)_{j}$ and the ground-
truth ${r}(\omega)_{j}$ for the $j$th sample across frequencies:

$\epsilon(\omega)_{j}=|10\log(|{r}(\omega)_{j}|^{2})-10\log(|\hat{r}(\omega)_{j}|^{2})|.$
(13)

Fig. 4: 95 percentile estimation error bounds of room average power response
estimators using a repetitive cross-validation.

Fig. 5: Average estimation error (bold lines) for three PCA-based estimators.
Shaded areas indicate estimation error standard deviation.

In particular, we focused on the 95th percentile error in dB level of the
repetitive cross-validation process as we aimed to check if the majority of
estimation processes deliver sufficient accuracy. In this evaluation, we chose
$N_{tr}=300$ and $N_{val}=100$. In Fig. 4, estimation performance of PCA-based
systems (including both the “Global PCA” method and local PCA methods “Local
PCA RT” and “Local PCA Roff” that represents systems based on local features
of RT30 and Low-frequency roll-off gain, respectively) is compared with the
linear LS regression estimator (“LS”) in Sec. 3.1 and the ensemble average
estimator ”Average”. For all PCA-based estimation methods, we selected a fixed
value of $K_{r}=32$ as we found that it consistently reconstructed raw room
average power response coefficients with sufficient accuracy using only the
first 32 orders of principal components. Note that the system of equations in
Eq. (9) is under-determined if $K_{s}<K_{r}$, hence, it is important to choose
$K_{s}\geq K_{r}$. Therefore, we choose $K_{s}=240$ for the global PCA method
and $K_{s}=80$ when training individual local PCA models.

From Fig. 4, it shows that a noticeably improved estimation performance can be
achieved by all PCA–based approaches compared to the linear LS regression
method using the same training dataset, which verifies advantages of the PCA
due to its efficiency and dimensionality reduction. Both of the local PCA
methods also outperform the global PCA method. Especially, ”Local PCA RT”
achieves an 95th percentile error bound that is within 5 dB across all
frequencies. It clearly shows that the division of the global training dataset
into sub-groups based on RT30 facilitates trained PCA models to capture the
variability of acoustic impulse responses due to local reverberation
conditions more accurately, which leads to superior estimation performance in
comparison to the system using low-frequency roll-off gains as local features.
Fig. 5 presents error regions across frequencies for all three PCA-based
estimation schemes, which demonstrates the same performance pattern as in Fig.
4. For the remaining part of the paper, the “local PCA RT” system is referred
to as the local PCA-based method.

###  4.3  Acoustic room compensation

To investigate how well the estimate of the room average power response
$\hat{r}(\omega)$ leads to a design that matches the reference “Trueplay”
system with the ground-truth room average power response, we conducted a
MUSHRA-style headphone listening test [29] with 11 listeners, including 6
professionally/self-reported trained listeners. For each listener, we randomly
selected 4 self-response/room average power response sets from the validation
dataset to simulate signals. Three stimulus signals were considered (including
vocal music, jazz music and movie clips), which were then presented over a
pre-equalized Beyerdynamic DT-880 headphone at 75 dBSPL.

In this evaluation, the following processing conditions are included: 1)
Reference system with ground-truth room average power responses; 2) Stock
system without room compensation processing; 3) “LS300”, denotes a filter
design system using linear LS regression estimator with $N_{tr}=300$;
4)“LSMove”, denotes the system using linear LS regression estimator with all
available training data from the Sonos Move database, i.e. $N_{tr}\approx
1000$; 5)“GPCA”, denotes the system using global PCA-based estimator with
$N_{tr}=300$; 6)“LPCA300”, denotes the system using local PCA-based estimator
with $N_{tr}=300$; 7)“LPCA150”, denotes the system using local PCA-based
estimator with $N_{tr}=150$. Subjects were asked to rate the perceived sound
similarity in comparison to provided reference samples.

Fig. 6 presents results of the MUSHRA test from all 11 subjects, which shows
that most subjects were able to identify the hidden reference signals.
However, several naïve listeners reported having difficulty distinguishing
multiple samples that sounded highly similar. For trained listeners, results
in Fig. 7 show that the stock condition was rated least similar to the
reference. The remaining six conditions were analyzed using a linear mixed
effects model [30]. Apart from hidden references, highest similarity ratings
were obtained by the proposed ”LPCA300” system ($p<0.001$ compared to “GPCA”,
“LS300” and even “LSMove” that uses a significantly larger training dataset).
In addition, “LPCA150” produces a close performance to ”LPCA300” with
$p>0.001$, which suggests that there is a potential to further reduce the size
of the training dataset while maintaining comparable perceptual improvements
for the PCA-based method. Overall, listening test results validates the
proposed local PCA-based room average power response estimation approach.

Fig. 6: Sound similarity ratings of the MUSHRA study evaluating room
compensation filters from all 11 listeners.

Fig. 7: Sound similarity ratings of the MUSHRA study from 6 trained listeners
(*** indicates p $<$ 0.001).

##  5 Conclusion

In this work, a room compensation design by leveraging an estimate of room
average power response is proposed. We present a PCA-based approach to predict
room average power responses using measured echo path self-response at the on-
board microphone. By exploiting reverberation time of echo path self-responses
and deriving individual local PCA models that capture the local variability of
impulse responses, the proposed method obtains an 95th percentile estimation
error within $5$ dB, which facilitates the proposed room compensation filter
to deliver a good sound similarity compared to the reference system. A MUSHRA-
style listening test was conducted and results confirm the effectiveness of
the approach.

## References

  * [1] Matti Karjalainen, Tuomoas Paatero, John N Mourjopoulos, and Panagiotis D Hatziantoniou,  “About room response equalization and dereverberation,”  in IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, 2005. IEEE, 2005, pp. 183–186.
  * [2] Matti Karjalainen, Aki Makivirta, Poju Antsalo, and Vesa Valimaki,  “Low-frequency modal equalization of loudspeaker-room responses,”  in Audio Engineering Society Convention 111. Audio Engineering Society, 2001.
  * [3] Matti Karjalainen, Paulo AA Esquef, Poju Antsalo, Aki Mäkivirta, and Vesa Välimäki,  “Frequency-zooming arma modeling of resonant and reverberant systems,”  Journal of the Audio Engineering Society, vol. 50, no. 12, pp. 1012–1029, 2002.
  * [4] John N Mourjopoulos,  “Comments on’analysis of traditional and reverberation-reducing methods of room equalization’,”  Journal of the Audio Engineering Society, vol. 51, no. 12, pp. 1186–1188, 2003.
  * [5] Wenyu Jin, W. Bastiaan Kleijn, and David Virette,  “Multizone soundfield reproduction using orthogonal basis expansion,”  in 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, 2013, pp. 311–315.
  * [6] Wenyu Jin and W. Bastiaan Kleijn,  “Theory and design of multizone soundfield reproduction using sparse methods,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 23, no. 12, pp. 2343–2355, 2015.
  * [7] Stefania Cecchi, Alberto Carini, and Sascha Spors,  “Room response equalization – a review,”  Applied Sciences, vol. 8, no. 1, 2018.
  * [8] Stephen T Neely and Jont B Allen,  “Invertibility of a room impulse response,”  The Journal of the Acoustical Society of America, vol. 66, no. 1, pp. 165–169, 1979.
  * [9] J Mourjopoulos,  “On the variation and invertibility of room impulse response functions,”  Journal of Sound and Vibration, vol. 102, no. 2, pp. 217–228, 1985\.
  * [10] Panagiotis D Hatziantoniou and John N Mourjopoulos,  “Errors in real-time room acoustics dereverberation,”  Journal of the Audio Engineering Society, vol. 52, no. 9, pp. 883–899, 2004.
  * [11] Wenyu Jin,  “Adaptive reverberation cancelation for multizone soundfield reproduction using sparse methods,”  in 2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2016, pp. 509–513.
  * [12] Sunil Bharitkar and Chris Kyriakakis,  Immersive audio signal processing,  Springer, 2006.
  * [13] Y. Haneda, Y. Kaneda, and N. Kitawaki,  “Common-acoustical-pole and residue model and its application to spatial interpolation and extrapolation of a room transfer function,”  IEEE Transactions on Speech and Audio Processing, vol. 7, no. 6, pp. 709–717, 1999.
  * [14] Rémi Mignot, Laurent Daudet, and François Ollivier,  “Room reverberation reconstruction: Interpolation of the early part using compressed sensing,”  IEEE Transactions on Audio, Speech, and Language Processing, vol. 21, no. 11, pp. 2301–2312, 2013.
  * [15] Samuel A. Verburg and Efren Fernandez-Grande,  “Reconstruction of the sound field in a room using compressive sensing,”  The Journal of the Acoustical Society of America, vol. 143, no. 6, pp. 3770–3779, 2018.
  * [16] Wenyu Jin and W. Bastiaan Kleijn,  “Multizone soundfield reproduction in reverberant rooms using compressed sensing techniques,”  in 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2014, pp. 4728–4732.
  * [17] Niccolò Antonello, Enzo De Sena, Marc Moonen, Patrick A. Naylor, and Toon van Waterschoot,  “Room impulse response interpolation using a sparse spatio-temporal representation of the sound field,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 25, no. 10, pp. 1929–1941, 2017.
  * [18] Majid Fozunbal, Ton Kalker, and Ronald W. Schafer,  “Multi-channel echo control by model learning,”  in Proc. Intern. Workshop Acoustic Echo and Noise Control, 2008\.
  * [19] Tomer Koren, Ronen Talmon, and Israel Cohen,  “Supervised system identification based on local pca models,”  in 2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2012, pp. 541–544.
  * [20] Wenyu Jin, Tim Schoof, and Henning Schepker,  “Individualized hear-through for acoustic transparency using pca-based sound pressure estimation at the eardrum,”  arXiv preprint arXiv:2110.04385, 2021.
  * [21] Floyd E. Toole,  Sound Reproduction: The Acoustics and Psychoacoustics of Loudspeakers and Rooms,  Routledge, USA, 3rd edition, 2017.
  * [22] “Trueplay spectral correction,” https://tech-blog.sonos.com/posts/trueplay-spectral-correction,  Accessed: 2022-04-06.
  * [23] Sonos Inc.,  “Tune your sonos speakers with trueplay,”  .
  * [24] S.-C. Pei and H.-S. Lin,  “Minimum-phase fir filter design using real cepstrum,”  IEEE Transactions on Circuits and Systems II: Express Briefs, vol. 53, no. 10, pp. 1113–1117, 2006.
  * [25] Peter Dodds,  “A flexible numerical optimization approach to the design of biquad filter cascades,”  Journal of the Audio Engineering Society, October 2020.
  * [26] Richard Lee,  “Simple arbitary iirs,”  Journal of the Audio Engineering Society, October 2008.
  * [27] G. W. Stewart,  “On the early history of the singular value decomposition,”  SIAM Review, vol. 35, no. 4, pp. 551–566, 1993.
  * [28] M. R. Schroeder,  “New method of measuring reverberation time,”  The Journal of the Acoustical Society of America, vol. 37, no. 3, pp. 409–412, 1965.
  * [29] M. Schoeffler et al.,  “webmushra — a comprehensive framework for web-based listening tests,”  Journal of Open Research Software, vol. 6, no. 1, pp. 8, 2018.
  * [30] D.A. Magezi,  “Linear mixed-effects models for within-participant psychology experiments: an introductory tutorial and free, graphical user interface (lmmgui),”  Frontiers in Psychology, 2015.