# Short-Sampled Blind Source Separation of Rotating Machinery Signals Based on Spectrum Correction

Authors: Xiangdong Huang; Xukang Jin; Haipeng Fu
Journal: Shock and Vibration
 (2016)

---

## Abstract

Nowadays, the existing blind source separation (BSS) algorithms in rotating machinery fault diagnosis can hardly meet the demand of fast response, high stability, and low complexity simultaneously. Therefore, this paper proposes a spectrum correction based BSS algorithm. Through the incorporation of FFT, spectrum correction, a screen procedure (consisting of frequency merging, candidate pattern selection, and single-source-component recognition), modifiedk-means based source number estimation, and mixing matrix estimation, the proposed BSS algorithm can accurately achieve harmonics sensing on field rotating machinery faults in case of short-sampled observations. Both numerical simulation and practical experiment verify the proposed BSS algorithm’s superiority in the recovery quality, stability to insufficient samples, and efficiency over the existing ICA-based methods. Besides rotating machinery fault diagnosis, the proposed BSS algorithm also possesses a vast potential in other harmonics-related application fields.


---

## Body

## 1. Introduction

As one of the most common classes of mechanical equipment, rotating machinery plays a significant role in industrial applications. Meanwhile, since it generally operates under harsh working conditions, it is likely to suffer from failures, which may cause the machinery to break down or decrease machinery service performance such as manufacturing quality and operation safety. Nowadays, rotating machineries in modern industry tend to be larger, more precise, and more automatic, which further increases the difficulty of the potential faults detection.Blind source separation (BSS), which can recover underlying sources from observations without the knowledge of the mixing system, is widely used in machinery fault diagnosis [1–5], speech recognition [6], wireless communication [7], and so on. Nowadays, BSS techniques applied in the machinery fault diagnosis mainly focus on two aspects: (
1
) removal of interferences and disturbances and (
2
) parameter modeling and feature detection for mechanical faults.On the one hand, as is known, rotating components (such as gears and bears) are the common and key components of modern machinery [8]. Affected by a lot of field factors (such as multiple motors that are fixed to the same structure or several fault events that happen simultaneously), the signal recorded from a sensor cannot solely reflect the operating state of a specific component. Furthermore, in industrial applications, these recorded signals are inevitably disrupted by the environment (ambient noise, other mechanical systems, etc.). Hence, BSS can act as an effective preprocessing procedure [9] to remove these interferences from other components or the disturbances arising from the environment. Wu et al. [10] proposed a BSS algorithm to remove the interferences of acoustic emission signals from a multiple cylinder diesel engine. In [11], an improved morphological component analysis (MCA) is proposed to diagnose compound faults of gearboxes. Cui et al. [4] put forward a null-space pursuit (NSP) BSS algorithm to diagnose compound faults of roller bearings.On the other hand, due to the effect of several rotor operations at some speeds, the signal recorded from a vibration sensor is mainly composed of multiple periodic harmonic components. For different categories of faults, the spectra of these recorded vibration signals exhibit distinct harmonics-related features. For example, a vibration signal caused by rotor misalignment is mainly characterized with the 2nd harmonic component [12]. The loosening of the bearing in the bearing block often generates components higher than 10th harmonic (even up to 20th harmonics). The fault of oil whirl [13] always gives rise to some subharmonics approximating half harmonic, and so forth. Hence, BSS is expected to accurately extract these harmonic features of individual sources. What is more, the model-based fault identification assumes that there exists a certain model to characterize a mechanical structure, in which the variation of model parameters can reflect the abnormal behaviors of the machinery system [14]. As a result, BSS can be utilized to identify the model parameters.Hence, a lot of studies of BSS problem have been made in the feature extracting and model identification fields. For example, sparse component analysis based [15] and independent component analysis (ICA) based [14] BSS methods were employed to estimate the vibration signals’ modal parameters. Following this, Žvokelj et al. [1] proposed the ensemble empirical mode decomposition based multiscale ICA (EEMD-MSICA) method and applied it into the bearing fault detection. Li et al. [16] proposed the supervised order tracking bounded component analysis (SOTBCA) based BSS algorithm for gear fault detection, which is suitable for dealing with the situation that the vibration signals do not satisfy the independent condition.To reduce the loss arising from fault accidents, it is urgently demanded in field operations that rotating machinery fault analysis should be as fast as possible. One possible solution is to implement the BSS in a short period of observations.However, these existing BSS methods can hardly work well in case of short-sampled observations. For example, the mainstream BSS method in rotating machinery fault diagnosis is the ICA [17]. A lot of ICA-based methods [18, 19], or improved ICA like second-order ICA [20], nonlinear adaptive ICA [21], and kernel ICA [22], are applied into the failure detection and analysis. As will be elaborated in this paper, ICA is likely to fall into nondeterministic solutions when provided only short-sampled observations. This arises from the fact that ICA is based on optimizing a kurtosis-related objective function. As a fourth-order cumulant statistic, the calculation of kurtosis needs to consume a large amount of samples. In fact, other statistics-based BSS methods, such as fourth-order-only blind identification (FOOBI) method [23] which is based on constructing high-order tensors, also exhibit poor performance in short-sampled situations.Hence, in this paper, we propose a novel blind source separation method which works well in both long observations and short observations. Due to the incorporation of spectrum correction and a phase coherence criterion, this BSS method can accurately extract harmonic features (frequency, amplitude, and phase) of individual sources. In case of short-sampled observations, which reduce the frequency resolution of fast Fourier transform (FFT) spectrum and thus deteriorate the picket-fence effect, the proposed BSS can also estimate harmonic parameters by means of spectrum correction. Moreover, a frequency screening procedure consisting of frequency merging, candidate pattern selection, and single-source-component recognition is able to exclude the interference between individual harmonics-related components. Therefore, unlike ICA or FOOBI method, the proposed BSS is competent in dealing with case of insufficient samples. In addition, the proposed BSS algorithm does not require thea priori source number. Both numerical simulation and practical experiment verify the proposed BSS algorithm’s superiority in efficiency and accuracy over the existing ICA-based methods.


## 2. Blind Source Separation Model

### 2.1. Temporal Model

ConsiderN underlying sources and M recording sensors. Suppose that the structure under investigation has a high rigidity, and the transmission delays in the mechanical structure are negligible compared to the sampling period [24]. In this case, the mixing system can be treated as an instantaneous one, which can be modeled as(1)
x
t
=
A
s
t
+
n
t
.In (1), s
(
t
)
=
[
s
1
(
t
)
,
s
2
(
t
)
,
…
,
s
N
(
t
)
]
T is the source vector, x
(
t
)
=
[
x
1
(
t
)
,
x
2
(
t
)
,
…
,
x
M
(
t
)
]
T is the observation vector, n
(
t
)
=
[
n
1
(
t
)
,
n
2
(
t
)
,
…
,
n
M
(
t
)
]
T is the additive noise vector, and A is the mixing matrix. The task of short-sampled BSS is to recover the sources s
1
(
t
)
~
s
N
(
t
) from the observations x
1
(
t
)
~
x
M
(
t
) without the knowledge of mixing matrix A in the small sample number situation.According to the relative relationship betweenN and M, the BSS problem can be divided into 2 conditions, the overdetermined or determined BSS (N
≤
M) and the underdetermined BSS (N
>
M). This paper focuses on the overdetermined condition.Since the vibration of some mechanical component stems from the rotation of the rotor,nth source can be formulated as a combination of individual harmonics: that is,(2)
s
n
t
=
∑
p
=
1
P
n
c
n
,
p
cos
⁡
2
π
f
n
,
p
t
+
θ
n
,
p
,where P
n is the number of components and c
n
,
p, f
n
,
p, and θ
n
,
p are the amplitude, frequency, and phase parameters of pth component of nth source, respectively.Based on this model, this paper aims to develop a BSS algorithm, which consumes a small amount of samples to estimate the mixing matrixA and recover all sources s
1
(
t
)
~
s
N
(
t
). Besides, it should be emphasized that, in industrial applications, the source number N is usually not known in advance. Therefore, this paper also addresses the problem of source number estimation.


### 2.2. Harmonics Based BSS Model

Combining (1) and (2), we can find that if an observation can be further linked to 3 harmonic-related parameters c
n
,
p, f
n
,
p, and θ
n
,
p, then the matrixA is expected to be estimated.Since a real signal contains two conjugate side spectra, we rewrites
n
(
t
) in (2) as(3)
s
n
t
=
s
~
n
t
+
s
~
n
∗
t
,where(4)
s
~
n
t
=
1
2
∑
p
=
1
P
n
c
n
,
p
exp
⁡
j
2
π
f
n
,
p
t
+
θ
n
,
p
.Further, if the harmonic frequencyf
n
,
p is far from direct component (DC), only a single side spectrum is enough to achieve BSS. In combination with (1), we have a frequency-domain model:(5)
X
~
f
=
A
S
~
f
.As is known, the ideal Fourier transform of a complex exponential signal is a dirac function. Hence, the spectrum ofnth source s
~
n
(
t
) in (4) is(6)
S
~
n
f
=
π
∑
p
=
1
P
n
c
n
,
p
δ
f
-
f
n
,
p
e
j
θ
n
,
p
.Denote the mixing matrix A as [
a
1
,
…
,
a
N
]. Substituting (6) into (5), we have(7)
X
~
1
f
⋮
X
~
m
f
⋮
X
~
M
f
=
π
a
1
,
…
,
a
N
∑
p
=
1
P
1
c
1
,
p
δ
f
-
f
1
,
p
e
j
θ
1
,
p
⋮
∑
p
=
1
P
n
c
n
,
p
δ
f
-
f
n
,
p
e
j
θ
n
,
p
⋮
∑
p
=
1
P
N
c
N
,
p
δ
f
-
f
N
,
p
e
j
θ
N
,
p
.To determine each column vector of the mixing matrixA, some particular frequency f
p
∗ which is only included in a single source and excluded by other sources, is considered; that is, f
p
∗ should satisfy(8)
f
p
∗
=
f
n
,
p
,
p
∈
1
,
…
,
P
n
f
p
∗
∉
f
n
¯
,
p
,
p
=
1
,
…
,
P
n
¯
,
n
¯
=
1
,
…
,
N
,
n
¯
≠
n
.Then, substituting (8) into (7) and combining with the sampling property of the dirac function “δ
(
·
)” in (7), we have(9)
X
~
1
f
p
∗
X
~
2
f
p
∗
⋮
X
~
M
f
p
∗
=
π
a
n
c
n
,
p
e
j
θ
n
,
p
.Then, it can be inferred from (9) that the frequency-domain vector X
~
(
f
p
∗
) corresponding to the component f
p
∗ is parallel to a
n. Hence, as long as sufficient single-source components f
p
∗ are collected, every column of the mixing matrix A can be sequentially determined.


### 2.3. Difficulty of Short-Sampled BSS

Note that (7) is an ideal Fourier model of the BSS system, in which the frequency f is a continuous variable. However, as is known, the ideal Fourier transform is unrealizable since it consumes infinite numbers of samples.In practice, the ideal Fourier transform is replaced by aL-point discrete Fourier transform (DFT) (“L” refers to the number of consumed samples), in which f in (7) only allows being one of L frequencies k
Δ
f, k
=
0,1
,
…
,
L
-
1 (Δ
f
=
f
s
/
L is the frequency resolution and f
s refers to the system sampling rate). Thus, the DFT spectrum of each observation will suffer from severe picket-fence effect.In addition, it is very likely that the frequencyf
n
,
p of nth source is not exactly the integer times of the DFT frequency resolution Δ
f
=
f
s
/
L, resulting in the fact that the dirac function δ
(
f
-
f
n
,
p
) in (7) cannot achieve an ideal sampling result. This deviation is also reflected in M observations’ DFT spectra X
~
m
(
k
Δ
f
)  (
m
=
1
,
…
,
M
), which exhibit the effect of the spectral leakage.Without loss of generality, denote the frequencyf
n
,
p of nth source as the summation of integer times and fractional time of Δ
f: that is,(10)
f
n
,
p
=
k
n
,
p
+
δ
n
,
p
Δ
f
,
k
n
,
p
∈
Z
+
,
δ
n
,
p
∈
-
0.5,0.5
.When the sample lengthL becomes smaller, the DFT frequency unit Δ
f
=
f
s
/
L gets larger and thus the DFT spectrum gets coarser. Limited by the picket-fence effect, in fact, the fractional item “δ
n
,
p
Δ
f” in (10) cannot be directly obtained from DFT bins and thus the frequency f
n
,
p has to be treated as the integer times of Δ
f (i.e., f
^
n
,
p
=
k
n
,
p
Δ
f), which corresponds to several peak DFT spectral bins of the observations. As a result, large deviation of frequency estimation inevitably occurs.Furthermore, as (7) shows, since an observation contains multiple components, severe interinterferences surely occur among distinct components when these frequency estimates f
^
n
,
p are inaccurate. As a result, the recovered spectrum of S
^
n
(
f
) is bound to be greatly different with the ideal spectrum, thereby increasing the BSS difficulty in the case of short-sampled observations.To overcome this difficulty, we introduce spectrum correction to solve this problem.


## 2.1. Temporal Model

ConsiderN underlying sources and M recording sensors. Suppose that the structure under investigation has a high rigidity, and the transmission delays in the mechanical structure are negligible compared to the sampling period [24]. In this case, the mixing system can be treated as an instantaneous one, which can be modeled as(1)
x
t
=
A
s
t
+
n
t
.In (1), s
(
t
)
=
[
s
1
(
t
)
,
s
2
(
t
)
,
…
,
s
N
(
t
)
]
T is the source vector, x
(
t
)
=
[
x
1
(
t
)
,
x
2
(
t
)
,
…
,
x
M
(
t
)
]
T is the observation vector, n
(
t
)
=
[
n
1
(
t
)
,
n
2
(
t
)
,
…
,
n
M
(
t
)
]
T is the additive noise vector, and A is the mixing matrix. The task of short-sampled BSS is to recover the sources s
1
(
t
)
~
s
N
(
t
) from the observations x
1
(
t
)
~
x
M
(
t
) without the knowledge of mixing matrix A in the small sample number situation.According to the relative relationship betweenN and M, the BSS problem can be divided into 2 conditions, the overdetermined or determined BSS (N
≤
M) and the underdetermined BSS (N
>
M). This paper focuses on the overdetermined condition.Since the vibration of some mechanical component stems from the rotation of the rotor,nth source can be formulated as a combination of individual harmonics: that is,(2)
s
n
t
=
∑
p
=
1
P
n
c
n
,
p
cos
⁡
2
π
f
n
,
p
t
+
θ
n
,
p
,where P
n is the number of components and c
n
,
p, f
n
,
p, and θ
n
,
p are the amplitude, frequency, and phase parameters of pth component of nth source, respectively.Based on this model, this paper aims to develop a BSS algorithm, which consumes a small amount of samples to estimate the mixing matrixA and recover all sources s
1
(
t
)
~
s
N
(
t
). Besides, it should be emphasized that, in industrial applications, the source number N is usually not known in advance. Therefore, this paper also addresses the problem of source number estimation.


## 2.2. Harmonics Based BSS Model

Combining (1) and (2), we can find that if an observation can be further linked to 3 harmonic-related parameters c
n
,
p, f
n
,
p, and θ
n
,
p, then the matrixA is expected to be estimated.Since a real signal contains two conjugate side spectra, we rewrites
n
(
t
) in (2) as(3)
s
n
t
=
s
~
n
t
+
s
~
n
∗
t
,where(4)
s
~
n
t
=
1
2
∑
p
=
1
P
n
c
n
,
p
exp
⁡
j
2
π
f
n
,
p
t
+
θ
n
,
p
.Further, if the harmonic frequencyf
n
,
p is far from direct component (DC), only a single side spectrum is enough to achieve BSS. In combination with (1), we have a frequency-domain model:(5)
X
~
f
=
A
S
~
f
.As is known, the ideal Fourier transform of a complex exponential signal is a dirac function. Hence, the spectrum ofnth source s
~
n
(
t
) in (4) is(6)
S
~
n
f
=
π
∑
p
=
1
P
n
c
n
,
p
δ
f
-
f
n
,
p
e
j
θ
n
,
p
.Denote the mixing matrix A as [
a
1
,
…
,
a
N
]. Substituting (6) into (5), we have(7)
X
~
1
f
⋮
X
~
m
f
⋮
X
~
M
f
=
π
a
1
,
…
,
a
N
∑
p
=
1
P
1
c
1
,
p
δ
f
-
f
1
,
p
e
j
θ
1
,
p
⋮
∑
p
=
1
P
n
c
n
,
p
δ
f
-
f
n
,
p
e
j
θ
n
,
p
⋮
∑
p
=
1
P
N
c
N
,
p
δ
f
-
f
N
,
p
e
j
θ
N
,
p
.To determine each column vector of the mixing matrixA, some particular frequency f
p
∗ which is only included in a single source and excluded by other sources, is considered; that is, f
p
∗ should satisfy(8)
f
p
∗
=
f
n
,
p
,
p
∈
1
,
…
,
P
n
f
p
∗
∉
f
n
¯
,
p
,
p
=
1
,
…
,
P
n
¯
,
n
¯
=
1
,
…
,
N
,
n
¯
≠
n
.Then, substituting (8) into (7) and combining with the sampling property of the dirac function “δ
(
·
)” in (7), we have(9)
X
~
1
f
p
∗
X
~
2
f
p
∗
⋮
X
~
M
f
p
∗
=
π
a
n
c
n
,
p
e
j
θ
n
,
p
.Then, it can be inferred from (9) that the frequency-domain vector X
~
(
f
p
∗
) corresponding to the component f
p
∗ is parallel to a
n. Hence, as long as sufficient single-source components f
p
∗ are collected, every column of the mixing matrix A can be sequentially determined.


## 2.3. Difficulty of Short-Sampled BSS

Note that (7) is an ideal Fourier model of the BSS system, in which the frequency f is a continuous variable. However, as is known, the ideal Fourier transform is unrealizable since it consumes infinite numbers of samples.In practice, the ideal Fourier transform is replaced by aL-point discrete Fourier transform (DFT) (“L” refers to the number of consumed samples), in which f in (7) only allows being one of L frequencies k
Δ
f, k
=
0,1
,
…
,
L
-
1 (Δ
f
=
f
s
/
L is the frequency resolution and f
s refers to the system sampling rate). Thus, the DFT spectrum of each observation will suffer from severe picket-fence effect.In addition, it is very likely that the frequencyf
n
,
p of nth source is not exactly the integer times of the DFT frequency resolution Δ
f
=
f
s
/
L, resulting in the fact that the dirac function δ
(
f
-
f
n
,
p
) in (7) cannot achieve an ideal sampling result. This deviation is also reflected in M observations’ DFT spectra X
~
m
(
k
Δ
f
)  (
m
=
1
,
…
,
M
), which exhibit the effect of the spectral leakage.Without loss of generality, denote the frequencyf
n
,
p of nth source as the summation of integer times and fractional time of Δ
f: that is,(10)
f
n
,
p
=
k
n
,
p
+
δ
n
,
p
Δ
f
,
k
n
,
p
∈
Z
+
,
δ
n
,
p
∈
-
0.5,0.5
.When the sample lengthL becomes smaller, the DFT frequency unit Δ
f
=
f
s
/
L gets larger and thus the DFT spectrum gets coarser. Limited by the picket-fence effect, in fact, the fractional item “δ
n
,
p
Δ
f” in (10) cannot be directly obtained from DFT bins and thus the frequency f
n
,
p has to be treated as the integer times of Δ
f (i.e., f
^
n
,
p
=
k
n
,
p
Δ
f), which corresponds to several peak DFT spectral bins of the observations. As a result, large deviation of frequency estimation inevitably occurs.Furthermore, as (7) shows, since an observation contains multiple components, severe interinterferences surely occur among distinct components when these frequency estimates f
^
n
,
p are inaccurate. As a result, the recovered spectrum of S
^
n
(
f
) is bound to be greatly different with the ideal spectrum, thereby increasing the BSS difficulty in the case of short-sampled observations.To overcome this difficulty, we introduce spectrum correction to solve this problem.


## 3. Spectrum Correction Based BSS

### 3.1. Spectrum Correction

In this paper, we apply the ratio-based spectrum correction method addressed in [25] to mth (m
=
1
,
…
,
M) observation to overcome the short-sampled difficulty. The spectrum correction consists of the following steps:(1)
Implement Hanning-windowed DFT on theL-length observation and acquire its DFT spectrum X
m
(
k
)  (k
=
0,1
,
…
,
L
-
1).
(2)
Collect all the peak indices ofX
m
(
k
). For the peak index k
m
,
p
m, p
m
=
1
,
…
,
P
m (P
m is the peak number of X
m
(
k
)), calculate the amplitude ratio v
m
,
p
m between X
m
(
k
m
,
p
m
) and its subpeak neighbor: that is,(11)
v
m
,
p
m
=
X
m
k
m
,
p
m
max
⁡
X
m
k
m
,
p
m
-
1
,
X
m
k
m
,
p
m
+
1
.
Further, a variableu
m
,
p
m can be obtained as(12)
u
m
,
p
m
=
2
-
v
m
,
p
m
1
+
v
m
,
p
m
.
(3)
Adjustu
m
,
p
m to estimate the fractional number as(13)
δ
^
m
,
p
m
=
u
m
,
p
m
,
if
X
m
k
m
,
p
m
+
1
>
X
m
k
m
,
p
m
-
1
-
u
m
,
p
m
,
else.
Then, the accurate frequency estimate is(14)
f
^
m
,
p
m
=
k
m
,
p
m
+
δ
^
m
,
p
m
f
s
L
.
(4)
Acquire the corrected amplitude estimated
^
m
,
p
m and phase estimate ϕ
^
m
,
p
m as(15)
d
^
m
,
p
m
=
2
π
δ
^
m
,
p
m
1
-
δ
^
m
,
p
m
2
X
m
k
m
,
p
m
sin
⁡
π
δ
^
m
,
p
m
,
ϕ
^
m
,
p
m
=
angle
X
m
k
m
,
p
m
-
π
δ
^
m
,
p
m
L
-
1
L
,
where “a
n
g
l
e
(
·
)” is the acquiring angle operator.After spectrum correction, 3 harmonic parameter sets{
f
^
m
,
p
m
}, {
d
^
m
,
p
m
}, and {
ϕ
^
m
,
p
m
} of mth observation (m
=
1
,
…
,
M) can be acquired.Further, as (8) and (9) demonstrate, for an estimated frequency f
^
m
,
p
m, only when it is included by a single source can it be utilized to estimate a column of the mixing matrix A. Hence, it is necessary to screen those single-source related frequencies from {
f
^
m
,
p
m
}, m
=
1
,
…
,
M.


### 3.2. Screening Single-Source Components

The proposed scheme of screening single-source components consists of 3 stages: frequency merging, candidate pattern selection, and single-source-component recognition.


#### 3.2.1. Frequency Merging

Its noteworthy that, affected by noise and interferences, even for the same single-source component, its frequency estimates of all the observations obtained by spectrum correction still exhibit tiny differences. Hence, a frequency merging procedure should be implemented.If we put all these frequency estimates together and sort them in an ascending order, the aforementioned frequency estimates of tiny differences tend to converge into a cluster. Assuming altogether thatP clusters are formed, without loss of generality, denote pth (p
=
1
,
…
,
P) cluster as {
f
^
p
,
q
,
q
=
1
,
…
,
Γ
p
} (Γ
p refers to pth cluster’s size). Then, Γ
p elements of this cluster can be merged by their average:(16)
f
¯
p
=
1
Γ
p
∑
q
=
1
Γ
p
f
^
p
,
q
.


#### 3.2.2. Candidate Pattern Selection

Theoretically, in terms of the BSS model (1), as long as the mixing matrix A does not contain zero elements, any source component should be included in all the observations. In other words, those corrected frequencies not contained by all the observations can be treated as fake components and should be removed.In practice, given a small thresholdε
>
0, for a merged frequency f
¯
p, if for each observation index m (m
=
1
,
…
,
M) there exists only one peak subscript p
¯
m satisfying(17)
f
^
m
,
p
¯
m
-
f
¯
p
<
ε
,
f
¯
p can be regarded as an effective component. Accordingly, in combination with (9), a pattern vector z
p
¯ relevant to this component’s M corrected parameter pairs (amplitude and phase) can be selected as a candidate vector to estimate a column of the matrix A; that is,(18)
z
p
¯
=
d
^
1
,
p
¯
1
e
j
ϕ
^
1
,
p
¯
1
⋮
d
^
m
,
p
¯
m
e
j
ϕ
^
m
,
p
¯
m
⋮
d
^
M
,
p
¯
M
e
j
ϕ
^
M
,
p
¯
M
,
p
¯
=
1
,
…
,
P
¯
.After candidate pattern selection, the number of merged frequencies is reduced fromP to P
¯.


#### 3.2.3. Single-Source-Component Criterion

In rotational machinery fault analysis, it is likely that multiple sources contain some common harmonic components (i.e., the overlapping frequencies). Obviously, these frequencies are not in accordance with (8) and (9) and should not be adopted to estimate the mixing matrix A. Hence, these overlapping components are invalid and should be removed from the candidate frequencies {
f
¯
p
}.Assume that amongP
¯ candidate frequencies {
f
¯
p
}, only P
∗ frequencies {
f
p
∗
,
p
∗
=
1
,
…
,
P
∗
} are single-source components. Since f
p
∗ only belongs to a single source, in combination with (9), its corresponding single-source-component vector z
p
∗ (i.e., the item X
~
(
f
p
∗
) in (9)) should be parallel to a column of the mixing matrix A.Furthermore, since the matrixA is real-valued, from (8) and (9), one can find that M phases of all the entries of X
~
(
f
p
∗
) originate from the same phase θ
n
,
p of a single source’s component f
p
∗ (i.e., f
n
,
p in (8)) and thus should be equal to each other.Thus, a single-source-component vectorz
p
∗ should exhibit two special properties:(1)
Its amplitude vector is parallel to a column of the mixing matrixA.
(2)
Its phase vector possesses a property of coherence, in which any two phase entries ofz
p
∗ should approximately point to the same direction. In other words, the following inequality of single-source-component criterion should be satisfied:(19)
1
C
M
2
∑
r
,
l
ϕ
^
r
,
p
∗
-
ϕ
^
l
,
p
∗
<
ξ
,
where1
≤
r, l
≤
M, r
≠
l, and ξ is a small positive value.


### 3.3. DB-Index Based Source Number Estimation andk-Means Clustering

If the source number “N” is known, one can directly employ a clustering algorithm (such as k-means clustering) on single-source-component vectors {
z
p
∗
,
p
∗
=
1
,
…
,
P
∗
} to estimate all the N columns of the mixing matrix A. However, in industrial applications, the source number “N” is usually unknown in advance. Therefore, this section combines DB-index [26] with k-means clustering to estimate N and A.Clearly, if the number of clusters is specified asI, then, the conventional k-means clustering algorithm can classify {
z
p
∗
,
p
∗
=
1
,
…
,
P
∗
} into I clusters C
i (i
=
1
,
…
,
I), whose entries can be denoted as(20)
C
i
=
z
˙
i
,
r
i
,
r
i
=
1,2
,
…
,
R
i
.The relationship between these clusters and the entire set of single-source-component vectors can be expressed as(21)
⋃
i
=
1
I
C
i
=
z
p
∗
,
p
∗
=
1
,
…
,
P
∗
,
∑
i
=
1
I
R
i
=
P
∗
.Davies Bouldin index (DB-index) is used to evaluate the appropriateness of data partitions [26] of a clustering algorithm. The definition of the DB-index is formulated as(22)
D
I
=
1
I
∑
i
=
1
I
max
j
⁡
G
i
+
G
j
M
i
,
j
,where G
i
,
G
j represents the dispersion measurement of two distinct groups C
i
,
C
j (assuming their cluster centers are c
i
,
c
j) and M
i
,
j refers to the similarity between these two groups. They are calculated with the following two formulas:(23)
G
i
=
1
R
i
∑
r
i
=
1
R
i
z
˙
i
,
r
i
-
c
i
,
M
i
,
j
=
c
i
-
c
j
.Apparently, on the one hand, the largerM
i
,
j is, the less the similarity between ith and jth clusters is, that is, the better the partition discrimination is. On the other hand, the smaller the dispersion degree G
i is, the higher the concentration degree of the group C
i is. As a result, the smaller the DB-index is, the more appropriate the data partition is. Therefore, the source number estimation can be realized by searching out the minimum DB-index of the k-means algorithm: that is,(24)
N
=
arg
⁡
min
I
⁡
D
I
.Once the source numberN is determined, the magnitude parts of cluster centers c
1
,
…
,
c
N of groups C
1
,
…
,
C
N generated by k-means algorithm can be directly treated as the columns of the mixing matrix estimate A
^.


### 3.4. Summary of the Proposed BSS Recovery Algorithm

Having obtained the overdetermined mixing matrix estimateA
^, the sources can be recovered by(25)
s
^
t
=
A
^
-
1
x
t
,where A
^
-
1 refers to the pseudoinverse of A
^.To summarize, the proposed BSS algorithm is listed as follows.Step 1.
Implement the procedure of spectrum correction addressed in Section3.1 on x
m
(
t
)
,
m
=
1
,
…
,
M, to acquire the corrected frequency set f
^
m
,
p
m, amplitude set d
^
m
,
p
m, and phase set ϕ
^
m
,
p
m.Step 2.
Merge the corrected frequenciesf
^
m
,
p
m using (16). Further, use (17) and (18) to acquire the candidate vectors {
z
p
-
,
p
-
=
1
,
…
,
P
-
}. Then, in terms of the screening criterion (19), pick out single-source-component vectors {
z
p
∗
,
p
∗
=
1
,
…
,
P
∗
} from these P
- candidate patterns.Step 3.
Implement the modifiedk-means clustering on single-source-component vectors {
z
p
∗
,
p
∗
=
1
,
…
,
P
∗
} to obtain the final estimate of the source number N and mixing matrix A
^.Step 4.
Calculate the pseudoinverseA
^
-
1 of A
^ and recover the source s
^
(
t
) by (25).


## 3.1. Spectrum Correction

In this paper, we apply the ratio-based spectrum correction method addressed in [25] to mth (m
=
1
,
…
,
M) observation to overcome the short-sampled difficulty. The spectrum correction consists of the following steps:(1)
Implement Hanning-windowed DFT on theL-length observation and acquire its DFT spectrum X
m
(
k
)  (k
=
0,1
,
…
,
L
-
1).
(2)
Collect all the peak indices ofX
m
(
k
). For the peak index k
m
,
p
m, p
m
=
1
,
…
,
P
m (P
m is the peak number of X
m
(
k
)), calculate the amplitude ratio v
m
,
p
m between X
m
(
k
m
,
p
m
) and its subpeak neighbor: that is,(11)
v
m
,
p
m
=
X
m
k
m
,
p
m
max
⁡
X
m
k
m
,
p
m
-
1
,
X
m
k
m
,
p
m
+
1
.
Further, a variableu
m
,
p
m can be obtained as(12)
u
m
,
p
m
=
2
-
v
m
,
p
m
1
+
v
m
,
p
m
.
(3)
Adjustu
m
,
p
m to estimate the fractional number as(13)
δ
^
m
,
p
m
=
u
m
,
p
m
,
if
X
m
k
m
,
p
m
+
1
>
X
m
k
m
,
p
m
-
1
-
u
m
,
p
m
,
else.
Then, the accurate frequency estimate is(14)
f
^
m
,
p
m
=
k
m
,
p
m
+
δ
^
m
,
p
m
f
s
L
.
(4)
Acquire the corrected amplitude estimated
^
m
,
p
m and phase estimate ϕ
^
m
,
p
m as(15)
d
^
m
,
p
m
=
2
π
δ
^
m
,
p
m
1
-
δ
^
m
,
p
m
2
X
m
k
m
,
p
m
sin
⁡
π
δ
^
m
,
p
m
,
ϕ
^
m
,
p
m
=
angle
X
m
k
m
,
p
m
-
π
δ
^
m
,
p
m
L
-
1
L
,
where “a
n
g
l
e
(
·
)” is the acquiring angle operator.After spectrum correction, 3 harmonic parameter sets{
f
^
m
,
p
m
}, {
d
^
m
,
p
m
}, and {
ϕ
^
m
,
p
m
} of mth observation (m
=
1
,
…
,
M) can be acquired.Further, as (8) and (9) demonstrate, for an estimated frequency f
^
m
,
p
m, only when it is included by a single source can it be utilized to estimate a column of the mixing matrix A. Hence, it is necessary to screen those single-source related frequencies from {
f
^
m
,
p
m
}, m
=
1
,
…
,
M.


## 3.2. Screening Single-Source Components

The proposed scheme of screening single-source components consists of 3 stages: frequency merging, candidate pattern selection, and single-source-component recognition.


### 3.2.1. Frequency Merging

Its noteworthy that, affected by noise and interferences, even for the same single-source component, its frequency estimates of all the observations obtained by spectrum correction still exhibit tiny differences. Hence, a frequency merging procedure should be implemented.If we put all these frequency estimates together and sort them in an ascending order, the aforementioned frequency estimates of tiny differences tend to converge into a cluster. Assuming altogether thatP clusters are formed, without loss of generality, denote pth (p
=
1
,
…
,
P) cluster as {
f
^
p
,
q
,
q
=
1
,
…
,
Γ
p
} (Γ
p refers to pth cluster’s size). Then, Γ
p elements of this cluster can be merged by their average:(16)
f
¯
p
=
1
Γ
p
∑
q
=
1
Γ
p
f
^
p
,
q
.


### 3.2.2. Candidate Pattern Selection

Theoretically, in terms of the BSS model (1), as long as the mixing matrix A does not contain zero elements, any source component should be included in all the observations. In other words, those corrected frequencies not contained by all the observations can be treated as fake components and should be removed.In practice, given a small thresholdε
>
0, for a merged frequency f
¯
p, if for each observation index m (m
=
1
,
…
,
M) there exists only one peak subscript p
¯
m satisfying(17)
f
^
m
,
p
¯
m
-
f
¯
p
<
ε
,
f
¯
p can be regarded as an effective component. Accordingly, in combination with (9), a pattern vector z
p
¯ relevant to this component’s M corrected parameter pairs (amplitude and phase) can be selected as a candidate vector to estimate a column of the matrix A; that is,(18)
z
p
¯
=
d
^
1
,
p
¯
1
e
j
ϕ
^
1
,
p
¯
1
⋮
d
^
m
,
p
¯
m
e
j
ϕ
^
m
,
p
¯
m
⋮
d
^
M
,
p
¯
M
e
j
ϕ
^
M
,
p
¯
M
,
p
¯
=
1
,
…
,
P
¯
.After candidate pattern selection, the number of merged frequencies is reduced fromP to P
¯.


### 3.2.3. Single-Source-Component Criterion

In rotational machinery fault analysis, it is likely that multiple sources contain some common harmonic components (i.e., the overlapping frequencies). Obviously, these frequencies are not in accordance with (8) and (9) and should not be adopted to estimate the mixing matrix A. Hence, these overlapping components are invalid and should be removed from the candidate frequencies {
f
¯
p
}.Assume that amongP
¯ candidate frequencies {
f
¯
p
}, only P
∗ frequencies {
f
p
∗
,
p
∗
=
1
,
…
,
P
∗
} are single-source components. Since f
p
∗ only belongs to a single source, in combination with (9), its corresponding single-source-component vector z
p
∗ (i.e., the item X
~
(
f
p
∗
) in (9)) should be parallel to a column of the mixing matrix A.Furthermore, since the matrixA is real-valued, from (8) and (9), one can find that M phases of all the entries of X
~
(
f
p
∗
) originate from the same phase θ
n
,
p of a single source’s component f
p
∗ (i.e., f
n
,
p in (8)) and thus should be equal to each other.Thus, a single-source-component vectorz
p
∗ should exhibit two special properties:(1)
Its amplitude vector is parallel to a column of the mixing matrixA.
(2)
Its phase vector possesses a property of coherence, in which any two phase entries ofz
p
∗ should approximately point to the same direction. In other words, the following inequality of single-source-component criterion should be satisfied:(19)
1
C
M
2
∑
r
,
l
ϕ
^
r
,
p
∗
-
ϕ
^
l
,
p
∗
<
ξ
,
where1
≤
r, l
≤
M, r
≠
l, and ξ is a small positive value.


## 3.2.1. Frequency Merging

Its noteworthy that, affected by noise and interferences, even for the same single-source component, its frequency estimates of all the observations obtained by spectrum correction still exhibit tiny differences. Hence, a frequency merging procedure should be implemented.If we put all these frequency estimates together and sort them in an ascending order, the aforementioned frequency estimates of tiny differences tend to converge into a cluster. Assuming altogether thatP clusters are formed, without loss of generality, denote pth (p
=
1
,
…
,
P) cluster as {
f
^
p
,
q
,
q
=
1
,
…
,
Γ
p
} (Γ
p refers to pth cluster’s size). Then, Γ
p elements of this cluster can be merged by their average:(16)
f
¯
p
=
1
Γ
p
∑
q
=
1
Γ
p
f
^
p
,
q
.


## 3.2.2. Candidate Pattern Selection

Theoretically, in terms of the BSS model (1), as long as the mixing matrix A does not contain zero elements, any source component should be included in all the observations. In other words, those corrected frequencies not contained by all the observations can be treated as fake components and should be removed.In practice, given a small thresholdε
>
0, for a merged frequency f
¯
p, if for each observation index m (m
=
1
,
…
,
M) there exists only one peak subscript p
¯
m satisfying(17)
f
^
m
,
p
¯
m
-
f
¯
p
<
ε
,
f
¯
p can be regarded as an effective component. Accordingly, in combination with (9), a pattern vector z
p
¯ relevant to this component’s M corrected parameter pairs (amplitude and phase) can be selected as a candidate vector to estimate a column of the matrix A; that is,(18)
z
p
¯
=
d
^
1
,
p
¯
1
e
j
ϕ
^
1
,
p
¯
1
⋮
d
^
m
,
p
¯
m
e
j
ϕ
^
m
,
p
¯
m
⋮
d
^
M
,
p
¯
M
e
j
ϕ
^
M
,
p
¯
M
,
p
¯
=
1
,
…
,
P
¯
.After candidate pattern selection, the number of merged frequencies is reduced fromP to P
¯.


## 3.2.3. Single-Source-Component Criterion

In rotational machinery fault analysis, it is likely that multiple sources contain some common harmonic components (i.e., the overlapping frequencies). Obviously, these frequencies are not in accordance with (8) and (9) and should not be adopted to estimate the mixing matrix A. Hence, these overlapping components are invalid and should be removed from the candidate frequencies {
f
¯
p
}.Assume that amongP
¯ candidate frequencies {
f
¯
p
}, only P
∗ frequencies {
f
p
∗
,
p
∗
=
1
,
…
,
P
∗
} are single-source components. Since f
p
∗ only belongs to a single source, in combination with (9), its corresponding single-source-component vector z
p
∗ (i.e., the item X
~
(
f
p
∗
) in (9)) should be parallel to a column of the mixing matrix A.Furthermore, since the matrixA is real-valued, from (8) and (9), one can find that M phases of all the entries of X
~
(
f
p
∗
) originate from the same phase θ
n
,
p of a single source’s component f
p
∗ (i.e., f
n
,
p in (8)) and thus should be equal to each other.Thus, a single-source-component vectorz
p
∗ should exhibit two special properties:(1)
Its amplitude vector is parallel to a column of the mixing matrixA.
(2)
Its phase vector possesses a property of coherence, in which any two phase entries ofz
p
∗ should approximately point to the same direction. In other words, the following inequality of single-source-component criterion should be satisfied:(19)
1
C
M
2
∑
r
,
l
ϕ
^
r
,
p
∗
-
ϕ
^
l
,
p
∗
<
ξ
,
where1
≤
r, l
≤
M, r
≠
l, and ξ is a small positive value.


## 3.3. DB-Index Based Source Number Estimation andk-Means Clustering

If the source number “N” is known, one can directly employ a clustering algorithm (such as k-means clustering) on single-source-component vectors {
z
p
∗
,
p
∗
=
1
,
…
,
P
∗
} to estimate all the N columns of the mixing matrix A. However, in industrial applications, the source number “N” is usually unknown in advance. Therefore, this section combines DB-index [26] with k-means clustering to estimate N and A.Clearly, if the number of clusters is specified asI, then, the conventional k-means clustering algorithm can classify {
z
p
∗
,
p
∗
=
1
,
…
,
P
∗
} into I clusters C
i (i
=
1
,
…
,
I), whose entries can be denoted as(20)
C
i
=
z
˙
i
,
r
i
,
r
i
=
1,2
,
…
,
R
i
.The relationship between these clusters and the entire set of single-source-component vectors can be expressed as(21)
⋃
i
=
1
I
C
i
=
z
p
∗
,
p
∗
=
1
,
…
,
P
∗
,
∑
i
=
1
I
R
i
=
P
∗
.Davies Bouldin index (DB-index) is used to evaluate the appropriateness of data partitions [26] of a clustering algorithm. The definition of the DB-index is formulated as(22)
D
I
=
1
I
∑
i
=
1
I
max
j
⁡
G
i
+
G
j
M
i
,
j
,where G
i
,
G
j represents the dispersion measurement of two distinct groups C
i
,
C
j (assuming their cluster centers are c
i
,
c
j) and M
i
,
j refers to the similarity between these two groups. They are calculated with the following two formulas:(23)
G
i
=
1
R
i
∑
r
i
=
1
R
i
z
˙
i
,
r
i
-
c
i
,
M
i
,
j
=
c
i
-
c
j
.Apparently, on the one hand, the largerM
i
,
j is, the less the similarity between ith and jth clusters is, that is, the better the partition discrimination is. On the other hand, the smaller the dispersion degree G
i is, the higher the concentration degree of the group C
i is. As a result, the smaller the DB-index is, the more appropriate the data partition is. Therefore, the source number estimation can be realized by searching out the minimum DB-index of the k-means algorithm: that is,(24)
N
=
arg
⁡
min
I
⁡
D
I
.Once the source numberN is determined, the magnitude parts of cluster centers c
1
,
…
,
c
N of groups C
1
,
…
,
C
N generated by k-means algorithm can be directly treated as the columns of the mixing matrix estimate A
^.


## 3.4. Summary of the Proposed BSS Recovery Algorithm

Having obtained the overdetermined mixing matrix estimateA
^, the sources can be recovered by(25)
s
^
t
=
A
^
-
1
x
t
,where A
^
-
1 refers to the pseudoinverse of A
^.To summarize, the proposed BSS algorithm is listed as follows.Step 1.
Implement the procedure of spectrum correction addressed in Section3.1 on x
m
(
t
)
,
m
=
1
,
…
,
M, to acquire the corrected frequency set f
^
m
,
p
m, amplitude set d
^
m
,
p
m, and phase set ϕ
^
m
,
p
m.Step 2.
Merge the corrected frequenciesf
^
m
,
p
m using (16). Further, use (17) and (18) to acquire the candidate vectors {
z
p
-
,
p
-
=
1
,
…
,
P
-
}. Then, in terms of the screening criterion (19), pick out single-source-component vectors {
z
p
∗
,
p
∗
=
1
,
…
,
P
∗
} from these P
- candidate patterns.Step 3.
Implement the modifiedk-means clustering on single-source-component vectors {
z
p
∗
,
p
∗
=
1
,
…
,
P
∗
} to obtain the final estimate of the source number N and mixing matrix A
^.Step 4.
Calculate the pseudoinverseA
^
-
1 of A
^ and recover the source s
^
(
t
) by (25).


## 4. Experiment

In this section, both numerical simulation of synthesis signals and practical mechanical diagnosis experiment are conducted to verify the performance of proposed BSS algorithm. As a comparison, the results of fast-ICA are also presented.


### 4.1. Numerical Simulation

Consider a3
×
2 mixing system expressed as(26)
A
=
0.59
0.89
0.97
0.47
0.51
0.53
.Two sources s
1
(
t
) and s
2
(
t
) are formulated as(27)
s
1
t
=
cos
⁡
2
π
100
t
+
π
9
+
1.4
cos
⁡
2
π
200
t
+
2
π
9
+
1.85
cos
⁡
2
π
400
t
+
π
3
,
s
2
t
=
1.7
cos
⁡
2
π
50
t
+
π
36
+
0.8
cos
⁡
2
π
80
t
+
5
π
18
+
1.2
cos
⁡
2
π
800
t
+
5
π
12
.The sampling rate was fixed as f
s
=
2000Hz and 4 cases of sample length (L
=
400
,
200
,
100
,
70) were taken into account. Since fast-ICA needs several iterative operations to optimize a kurtosis-related objective function, which starts from a random initialization on the demixing matrix, it is likely to fall into failure in case of insufficient samples. Hence, for each sample length case, 1000 trials were conducted. The times of successful trials of both BSS algorithms were recorded in Table 1. Moreover, among these successful trials, correlation coefficients between the recovered signals and the sources were statistically averaged and also listed in Table 1. Figures 1 and 2 present the recovered results of these two BSS algorithms in case of long observations (L
=
400), while Figures 3 and 4 present the short observation case (L
=
70).Table 1
The recovery correlation coefficients and successful times of numerical experiment.
Sample lengthL
400
200
100
70
ICA
s
1
(
t
)
Successful times
999
685
575
544
Correlation coefficient
0.9724
0.9961
0.9727
0.9093
s
2
(
t
)
Successful times
1000
456
645
620
Correlation coefficient
0.9987
0.9386
0.9027
0.8334
Proposed
s
1
(
t
)
Successful times
1000
1000
1000
1000
Correlation coefficient
1.0000
1.0000
1.0000
1.0000
s
2
(
t
)
Successful times
1000
1000
1000
1000
Correlation coefficient
1.0000
1.0000
1.0000
1.0000Figure 1
Numerical results of fast-ICA (L
=
400).Figure 2
Numerical results of proposed BSS algorithm (L
=
400).Figure 3
Numerical results of fast-ICA (L
=
70).Figure 4
Numerical results of proposed BSS algorithm (L
=
70).As Figures1 and 2 depict, both the fast-ICA and proposed algorithm can acquire high-quality recovered waveforms in case of long observations (L
=
400, limited by page layout, only half-duration waveforms are plotted). However, when the sample length reduces into L
=
70, one can observe that obvious distortions appear in the waveforms recovered by fast-ICA in Figure 3. In contrast, there exist no distortions in the recovered waveforms in Figure 4, reflecting that the proposed BSS algorithm outperforms fast-ICA in dealing with insufficient samples.Table1 shows that as the sample length decreases, the times of successful recovery of fast-ICA decline sharply, and the average correlation coefficient also tends to be slightly smaller, accordingly. In contrast, as Table 1 lists, all the trials of the proposed BSS algorithm for different sample lengths are successfully conducted and all correlation coefficients remain 1. This is because, unlike fast-ICA, the proposed BSS algorithm is based on spectrum correction related harmonics analysis rather than statistical analysis and thus it is insensitive to the sample length.


### 4.2. Mechanical Diagnosis Experiment

In this section, two practical fault signalss
1
(
t
), s
2
(
t
) collected from field rotating machineries are treated as sources. s
1
(
t
) is an imbalance fault signal with the rotating frequency 89.6853 Hz, and s
2
(
t
) is a misalignment fault signal with the rotating frequency 102.8811 Hz. The mixing system is the same as the matrix A in (26). Different sample lengths (L
=
400
,
200
,
100
,
70) were considered. In each case, 1000 trials were conducted. Figures 5–8 present the recovery results of both BSS algorithms. Table 2 lists their recovery performance indexes.Table 2
The recovery correlation coefficients and successful times of numerical experiment.
Sample lengthL
400
200
100
70
ICA
s
1
(
t
)
Successful times
908
922
912
874
Correlation coefficient
0.9734
0.9390
0.9778
0.9724
s
2
(
t
)
Successful times
912
902
900
734
Correlation coefficient
0.9270
0.9140
0.9193
0.9316
Proposed
s
1
(
t
)
Successful times
1000
1000
1000
1000
Correlation coefficient
0.9995
0.9973
0.9992
0.9870
s
2
(
t
)
Successful times
1000
1000
1000
1000
Correlation coefficient
0.9908
0.9803
0.9727
0.9881Figure 5
Practical results of fast-ICA (L
=
400).Figure 6
Practical results of the proposed BSS algorithm (L
=
400).Figure 7
Practical result of fast-ICA (L
=
70).Figure 8
Practical result of proposed BSS algorithm (L
=
70).From Figures5 and 6, one can see that, just like the recovery of synthesis signals in (27), both the fast-ICA and the proposed BSS algorithm can achieve excellent recovery results in the long-sample situation (L
=
400). Nevertheless, when it comes to the short-sample situation (L
=
70), the proposed BSS algorithm exhibits better performance than the fast-ICA does.From Table2, one can see that, as the sample length decreases from L
=
400 to 70, the proposed BSS algorithm’s superiority over fast-ICA becomes more obvious. In particular, due to the effect of field noise, the correlation coefficients resulting from the proposed algorithm do not remain 1 but approximate to 1. Hence, the proposed BSS algorithm outperforms the fast-ICA in rotating machinery fault diagnosis.


## 4.1. Numerical Simulation

Consider a3
×
2 mixing system expressed as(26)
A
=
0.59
0.89
0.97
0.47
0.51
0.53
.Two sources s
1
(
t
) and s
2
(
t
) are formulated as(27)
s
1
t
=
cos
⁡
2
π
100
t
+
π
9
+
1.4
cos
⁡
2
π
200
t
+
2
π
9
+
1.85
cos
⁡
2
π
400
t
+
π
3
,
s
2
t
=
1.7
cos
⁡
2
π
50
t
+
π
36
+
0.8
cos
⁡
2
π
80
t
+
5
π
18
+
1.2
cos
⁡
2
π
800
t
+
5
π
12
.The sampling rate was fixed as f
s
=
2000Hz and 4 cases of sample length (L
=
400
,
200
,
100
,
70) were taken into account. Since fast-ICA needs several iterative operations to optimize a kurtosis-related objective function, which starts from a random initialization on the demixing matrix, it is likely to fall into failure in case of insufficient samples. Hence, for each sample length case, 1000 trials were conducted. The times of successful trials of both BSS algorithms were recorded in Table 1. Moreover, among these successful trials, correlation coefficients between the recovered signals and the sources were statistically averaged and also listed in Table 1. Figures 1 and 2 present the recovered results of these two BSS algorithms in case of long observations (L
=
400), while Figures 3 and 4 present the short observation case (L
=
70).Table 1
The recovery correlation coefficients and successful times of numerical experiment.
Sample lengthL
400
200
100
70
ICA
s
1
(
t
)
Successful times
999
685
575
544
Correlation coefficient
0.9724
0.9961
0.9727
0.9093
s
2
(
t
)
Successful times
1000
456
645
620
Correlation coefficient
0.9987
0.9386
0.9027
0.8334
Proposed
s
1
(
t
)
Successful times
1000
1000
1000
1000
Correlation coefficient
1.0000
1.0000
1.0000
1.0000
s
2
(
t
)
Successful times
1000
1000
1000
1000
Correlation coefficient
1.0000
1.0000
1.0000
1.0000Figure 1
Numerical results of fast-ICA (L
=
400).Figure 2
Numerical results of proposed BSS algorithm (L
=
400).Figure 3
Numerical results of fast-ICA (L
=
70).Figure 4
Numerical results of proposed BSS algorithm (L
=
70).As Figures1 and 2 depict, both the fast-ICA and proposed algorithm can acquire high-quality recovered waveforms in case of long observations (L
=
400, limited by page layout, only half-duration waveforms are plotted). However, when the sample length reduces into L
=
70, one can observe that obvious distortions appear in the waveforms recovered by fast-ICA in Figure 3. In contrast, there exist no distortions in the recovered waveforms in Figure 4, reflecting that the proposed BSS algorithm outperforms fast-ICA in dealing with insufficient samples.Table1 shows that as the sample length decreases, the times of successful recovery of fast-ICA decline sharply, and the average correlation coefficient also tends to be slightly smaller, accordingly. In contrast, as Table 1 lists, all the trials of the proposed BSS algorithm for different sample lengths are successfully conducted and all correlation coefficients remain 1. This is because, unlike fast-ICA, the proposed BSS algorithm is based on spectrum correction related harmonics analysis rather than statistical analysis and thus it is insensitive to the sample length.


## 4.2. Mechanical Diagnosis Experiment

In this section, two practical fault signalss
1
(
t
), s
2
(
t
) collected from field rotating machineries are treated as sources. s
1
(
t
) is an imbalance fault signal with the rotating frequency 89.6853 Hz, and s
2
(
t
) is a misalignment fault signal with the rotating frequency 102.8811 Hz. The mixing system is the same as the matrix A in (26). Different sample lengths (L
=
400
,
200
,
100
,
70) were considered. In each case, 1000 trials were conducted. Figures 5–8 present the recovery results of both BSS algorithms. Table 2 lists their recovery performance indexes.Table 2
The recovery correlation coefficients and successful times of numerical experiment.
Sample lengthL
400
200
100
70
ICA
s
1
(
t
)
Successful times
908
922
912
874
Correlation coefficient
0.9734
0.9390
0.9778
0.9724
s
2
(
t
)
Successful times
912
902
900
734
Correlation coefficient
0.9270
0.9140
0.9193
0.9316
Proposed
s
1
(
t
)
Successful times
1000
1000
1000
1000
Correlation coefficient
0.9995
0.9973
0.9992
0.9870
s
2
(
t
)
Successful times
1000
1000
1000
1000
Correlation coefficient
0.9908
0.9803
0.9727
0.9881Figure 5
Practical results of fast-ICA (L
=
400).Figure 6
Practical results of the proposed BSS algorithm (L
=
400).Figure 7
Practical result of fast-ICA (L
=
70).Figure 8
Practical result of proposed BSS algorithm (L
=
70).From Figures5 and 6, one can see that, just like the recovery of synthesis signals in (27), both the fast-ICA and the proposed BSS algorithm can achieve excellent recovery results in the long-sample situation (L
=
400). Nevertheless, when it comes to the short-sample situation (L
=
70), the proposed BSS algorithm exhibits better performance than the fast-ICA does.From Table2, one can see that, as the sample length decreases from L
=
400 to 70, the proposed BSS algorithm’s superiority over fast-ICA becomes more obvious. In particular, due to the effect of field noise, the correlation coefficients resulting from the proposed algorithm do not remain 1 but approximate to 1. Hence, the proposed BSS algorithm outperforms the fast-ICA in rotating machinery fault diagnosis.


## 5. Conclusion

This paper proposes a novel blind source separation algorithm based on spectrum correction. Both numerical simulation and practical experiment verify the proposed BSS algorithm’s excellent performance. In general, this algorithm possesses the following 4 merits:(1)
Compared to classical fast-ICA algorithm, the proposed algorithm can achieve a higher-quality source recovery even in case of short-sample observations. This meets the demand of fast response of the rotating machinery fault analysis.
(2)
The spectrum correction involved in the proposed algorithm does well in harmonics information extraction and thus is especially suitable for rotating machinery fault analysis. As is known, most of these faults arise from the rotor malfunction, which generates a lot of rotating-frequency related harmonics.
(3)
The proposed BSS algorithm can accurately determine the underlying source number by means of the modifiedk-means clustering, which is in accordance with practical situation of rotating machinery operations.
(4)
Unlike fast-ICA, the proposed BSS algorithm does not involve random initialization and iterative operations and thus possesses a higher stability and lower complexity, which enhances the reliability and efficiency of the rotating machinery fault analysis.In fact, besides rotating machinery fault analysis, harmonics analysis is also frequently encountered in a lot of fields such as power harmonics analysis, channel estimation in communication, radar, and sonar. Hence, the proposed BSS algorithm possesses a vast potential in a wide range of applications.



---

*Source: 9564938-2016-09-20.xml*