# Wavelet Packet Division Multiplexing (WPDM)-Aided Industrial WSNs ††thanks:
This material is based upon work supported by Science Foundation Ireland (SFI)
and is co-funded under the European Regional Development Fund under Grant
Numbers 13/RC/2077 and 13/RC/2077-P2.
©2023 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

_Indrakshi Dey and Nicola Marchetti_

###### Abstract

Industrial Internet-of-Things (IIoT) involve multiple groups of sensors, each
group sending its observations on a particular phenomenon to a central
computing platform over a multiple access channel (MAC). The central platform
incorporates a decision fusion center (DFC) that arrives at global decisions
regarding each set of phenomena by combining the received local sensor
decisions. Owing to the diverse nature of the sensors and heterogeneous nature
of the information they report, it becomes extremely challenging for the DFC
to denoise the signals and arrive at multiple reliable global decisions
regarding multiple phenomena. The industrial environment represents a specific
indoor scenario devoid of windows and filled with different noisy electrical
and measuring units. In that case, the MAC is modelled as a large-scale
shadowed and slowly-faded channel corrupted with a combination of Gaussian and
impulsive noise. The primary contribution of this paper is to propose a
flexible, robust and highly noise-resilient multi-signal transmission
framework based on Wavelet packet division multiplexing (WPDM). The local
sensor observations from each group of sensors are waveform coded onto wavelet
packet basis functions before reporting them over the MAC. We assume a multi-
antenna DFC where the waveform-coded sensor observations can be separated by a
bank of linear filters or a correlator receiver, owing to the orthogonality of
the received waveforms. At the DFC we formulate and compare fusion rules for
fusing received multiple sensor decisions, to arrive at reliable conclusions
regarding multiple phenomena. Simulation results show that WPDM-aided wireless
sensor network (WSN) for IIoT environments offer higher immunity to noise by
more than 10 times over performance without WPDM in terms of probability of
false detection.

##  1 Introduction

Industrial Internet-of-Things (IoT)/Industry 4.0 conceptually will
interconnect everything within an industry, including employees, machines and
products. They will interact with each other to ensure seamless production
without any human intervention. The idea is to deploy different kinds of
sensors and related devices that will collect information regarding customer
requirements, conditions of the operating machines, and the surrounding
conditions in which they operate. Through groups of heterogeneous sensors,
machines will communicate among themselves through a virtual social network
and will communicate with human beings (managers, owners, customers) through
an actual digital network [1, 2, 3].

Deploying a large-scale wireless sensor network (WSN) that can enable the
deployment of industrial IoT is extremely challenging owing to different
factors. All sensors have to transmit their decisions simultaneously over a
multiple access channel (MAC) to a decision fusion center (DFC) owing to
limited bandwidth and as a result, the DFC receives a superposition of
heterogeneous sensor decisions [4]. Such superposition results in interfering
sensor signals corrupted with noise. Besides, the wireless MAC suffers from
random time-varying fading and shadowing [5]. Additionally, industrial
environment is different from traditional indoor environments like homes,
offices and buildings. Owing to the large dimensions and presence of noisy
instruments and machines, industrial environments suffer from impulsive noise
in addition to background and electronic Gaussian noise [6].

In order to combat the above-mentioned challenges concomitant with large-scale
industrial sensing networks, we propose a robust wireless sensor communication
system that combines wavelet packet division multiplexing (WPDM) on the
transmit side with distributed multi-antenna based decision fusion (DF) on the
receive side. Implementing WPDM has been recommended in [7, 8] for providing
substantial immunity to both impulsive and Gaussian noise. Furthermore,
deploying multiple antennas at the DFC can improve fusion performance over
deep faded and shadowed MAC [9, 10]. Therefore, combining WPDM and multi-
antenna DF can provide considerable reduction in vulnerability to noise,
fading and shadowing in a large-scale industrial sensing network.

The primary contribution of this paper is to propose i) waveform coding of
local sensor observations from each group of sensors onto wavelet packet basis
functions before transmitting them over wireless MAC and ii) fusing the
decisions at multi-antenna DFC using linear filter-bank or correlator receiver
for separating the received waveform-coded sensor observations, with the aim
of incorporating resilience to noise, fading and shadowing. Another major
contribution of the paper is to evaluate the performance of our proposed
system against a realistic industrial environment, which suffers from a
combination of Gaussian distributed additive background noise and impulsive
noise. We formulate impulsive noise as either Middleton Class A distributed or
Bernoulli-Gaussian distributed [11, 12]. Local sensor decisions from each
sensor group within an industrial sensing environment are coded with wavelet
packet basis functions, and each group is waveform-coded at each level of
binary tree structure thereby formulating a WPDM system. The WPDM multiplexed
signals are transmitted over the MAC to a DFC equipped with multiple antennas.
At the DFC, we formulate linear multi-rate filter-based optimum and sub-
optimum fusion rules generalized to the WPDM set-up. Performance of our
formulated system is simulated and compared for three different wavelet
scaling functions [13], Haar, Shannon and Piecewise linear spline, in presence
of varying degrees of impulse noise.

The rest of the paper is organized as follows. Section II introduces the
system design, Section III formulates the filters and the fusion rules,
Section IV analyzes performance of our proposed system, while concluding
remarks are provided in Section V.

##  2 System Design

The infographic representation for a large-scale industrial sensing network
consisting of heterogeneous sensor groups is provided in Fig. 1. We assume
that there are different groups of sensors; $z\triangleq\\{1,\dotso,Z\\}$. For
example, one group of sensors could be mobility sensors reporting if there is
a moving target or not, another group could sense if a certain harmful
chemical is present or not, another group could sense if carcinogenic material
like asbestos is present or not etc. Each of the $Z$ sensor groups transmit
their observations to a DFC equipped with $N$ antennas. The observations are
waveform coded onto wavelet packet basis functions, thereby implementing a
wavelet packet division multiplexing (WPDM) structure.

###  2.1 System Model

We design a binary tree structure with $L$ number of levels such that
$Z=2^{L}$ or $L=\log_{2}Z$. If $Z$ is odd, then $L=\lceil\log_{2}Z\rceil$,
where $\lceil\rceil$ denotes the ceiling function. Let us also assume that the
$z$th sensor group consists of $M$ transmit sensors
$(m\triangleq\\{1,\dotso,M\\})$. Each of the $M$ sensors in each group
transmits its independent observation over a $T_{l}$ symbol duration using a
time division multiplexed system, where $l$ is the level of the binary tree
structure of the WPDM system, $l\geq 0;l\triangleq\\{1,\dotso,L\\}$. Here
$T_{l}=2^{l}T_{0}$ where $T_{0}$ is the time interval between symbols. It is
worth-mentioning here that the tree can be pruned and developed depending on
the total number of sensor groups present.

Figure 1: Representation of WPDM-aided large-scale industrial sensing network
with sensors belonging to the $z$th group sending their observations to a
multi-antenna aided DFC.

An example scenario with $Z=4$ sensor groups is presented in the inset of Fig.
1. In Fig. 1, $\phi_{lz,m}$ is the family of scaling functions derived from
the family of wavelets such that these functions form the terminals (leaves)
of the tree. These functions can also be referred to as the wavelet packet
basis function or simply a wavelet packet. For a given tree structure, the
scaling function, $\phi_{lz,m}$ can be constituted using,

$\displaystyle\mathbf{\Phi}_{l}^{z}=[\phi_{lz,1},\phi_{lz,2},\dotso,\phi_{lz,M}]=\sum_{q=1}^{Q}\mathbf{f}_{l}^{z}[q]\leavevmode\nobreak\
\phi_{00}(qT_{0})$ (1)

where
$\mathbf{f}_{l}^{z}[q]=[f_{l1}^{z}[q],f_{l2}^{z}[q],\dotso,f_{lM}^{z}[q]]$ is
the equivalent sequence filter built from the combination of the up-sampling
and down-sampling finite impulse response (FIR) filters [14] and $\phi_{00}$
is the root of the tree structure. It is worth-mentioning here that
$\mathbf{\Phi}_{l}^{z}$ functions are self and mutually orthogonal at integral
multiples of $T_{l}$ and have a finite duration.

Each of the $M$ sensors in the $z$th sensor group sends its local observation
based on a binary hypothesis test, i.e. $\mathcal{H}_{l1}^{z}$, if the
intended target is present, and $\mathcal{H}_{l0}^{z}$ if it is absent or at a
level below the recommended threshold. The individual sensor observations are
mapped onto a binary phase shift keying (BPSK)-modulated symbol, such that,
$x^{z}_{lm}\triangleq\\{+1,-1\\}$. Here the $m$th sensor observation to a
$z$th group is mapped onto the $l$th level of the binary tree such that,

$\displaystyle
s^{z}_{lm}=x^{z}_{lm}\sum_{q=1}^{Q}\sum_{k=1}^{K}f^{z}_{lm}[k]\leavevmode\nobreak\
\phi_{00}((Zq+k)T_{0})$ (2)

where $1\leq K\leq Q/2$; $Q$ is the (even) total number of filter
coefficients, and $K$ is the total number of zeros of the FIR filter,
$F^{z}_{lm}(\zeta)$ at $\zeta=-1$ where $F^{z}_{lm}(\cdot)$ is the
$Z$-transform of $f^{z}_{lm}[\cdot]$. Therefore,
$s^{z}_{lm}=x^{z}_{lm}\sigma^{z}_{lm}[k,q]$ where
$\sigma^{z}_{lm}[k,q]=\sum_{q=1}^{Q}\sum_{k=1}^{K}f^{z}_{lm}[k]\leavevmode\nobreak\
\phi_{00}((Zq+k)T_{0})$ which makes up the encoded set of sensor decisions on
the $l$th level as,

$\displaystyle\mathbf{s}^{z}_{l}\triangleq[s^{z}_{l1},s^{z}_{l2},\dotso
s^{z}_{lM}]^{t}=[x^{z}_{l1}\sigma^{z}_{l1}[k,q],\dotso,x^{z}_{lM}\sigma^{z}_{lM}[k,q]]^{t}.$
(3)

Owing to the orthogonality between the wavelet packet basis functions
$\mathbf{\Phi}^{z}_{l}$, different groups of sensors
($z\triangleq\\{1,\dotso,Z\\}$) can be time division multiplexed, while arrays
can be implemented at the DFC for fusing sensor observations and obtaining a
final decision (refer to Fig. 1). At the $(l,z)$th terminal (the $z$th group
of sensors at the $l$th level), $\mathbf{s}^{z}_{l}$ is time-multiplexed with
other $(Z-1)$ group of sensors to be communicated over the wireless channel.
The bandwidth of the channel is however assumed to be sufficiently large as
not to cause any significant distortion of the waveform of the transmitted
sensor observations.

###  2.2 Signal and Channel Models

The signal received at the DFC can be written in the discrete-time matrix form
over the $l$th level as,

$\displaystyle\mathbf{Y}^{z}_{l}=\sqrt{\rho_{l}}\mathbf{G}^{z}_{l}\mathbf{s}^{z}_{l}+\mathbf{e}^{z}_{l}\triangleq[y^{z}_{l1},y^{z}_{l2},\dotso
y^{z}_{lM}]^{t}$ (4)

where $\mathbf{Y}^{z}_{l}$ is the received signal vector from the $z$th group
of sensors that has been mapped on the $l$th level, $\mathbf{G}^{z}_{l}$ is
the channel matrix between the $z$th group of sensors and the DFC,
$\mathbf{e}^{z}_{l}\sim\mathcal{N}_{\mathbb{C}}\big{(}\mathbf{0}_{N\times
T_{l}},\leavevmode\nobreak\ \Sigma^{2}_{e}\mathbf{I}_{N\times T_{l}}\big{)}$
is the noise vector and $\mathcal{N}_{\mathbb{C}}(\mu,\varphi)$ denotes
complex normal distribution with mean $\mu$ and co-variance $\varphi$. Here
$\mathbf{e}^{z}_{l}$ accounts for heterogeneous noise encountered in
industrial environment including both impulsive and Gaussian noise. Therefore,
$\Sigma^{2}_{e}=\Sigma_{w}^{2}+\sum_{\gamma=1}^{5}\Sigma_{I}^{2}/\gamma A$ for
Middleton Class A noise and
$\Sigma^{2}_{e}=\Sigma_{w}^{2}+M\Sigma_{I}^{2}/\varrho$ for Bernoulli-Gaussian
distributed noise, where $\Sigma_{w}^{2}$ is the variance of the Gaussian
noise and $\Sigma_{I}^{2}$ is the variance of the impulsive noise,
$A=\eta\tau/\tau_{0}$, if the impulsive noise follows the Middleton Class A
distribution, $\eta$ is the average number of impulses per second, $\tau_{0}$
is the interval between two impulses, and $\tau$ is the average duration of
each impulse. If the impulsive noise is Bernoulli-Gaussian distributed, it
will follow that distribution with a probability of $\varrho$. Each element of
the channel matrix $\mathbf{G}^{z}_{l}$ is given by
$g^{z}_{n,lm}=\sqrt{\lambda_{lm}}h^{z}_{n,lm}$ with the geometric attenuation
and the shadow fading $\lambda_{m}$ and fast fading coefficients
$h^{z}_{n,lm}$. Consequently,
$\mathbf{G}^{z}_{l}=\mathbf{H}^{z}_{l}\sqrt{\mathbf{D}_{l}}\in\mathbb{C}^{N\times
M}$ denotes the matrix of the generic channel coefficients,
$\mathbf{H}^{z}_{l}\in\mathbb{C}^{N\times M}$ denotes the matrix of the fast
fading coefficients, and $\mathbf{D}_{l}$ is a diagonal matrix with
$d_{lm,lm}=\lambda_{lm}$.

At the DFC, by applying a multirate filter bank with filters having impulse
responses $\mathbf{f}^{z}_{l}[-q]$, a process called wavelet packet
reconstruction, we can arrive at the linearized equivalent received signal
model as,

$\displaystyle\hat{\mathbf{Y}}^{z}_{l}=\sqrt{\rho_{l}}\hat{\mathbf{G}}^{z}_{l}\mathbf{r}^{z}_{l}\leavevmode\nobreak\
\hat{\mathbf{\Theta^{z}_{l}}}[k,q]+\mathbf{e}^{z}_{l}$ (5)

where $\hat{\mathbf{G}}^{z}_{l}$ is the estimated channel matrix. The DFC
estimates the channel state information (CSI), where half of the coherence
interval of the channel is used for training to estimate the channel and
establish the frequency and timing synchronization. A recommended channel
estimator is the minimum mean-squared error (MMSE) estimator. In this paper,
we do not separately estimate the channel and calculate the channel estimation
error. We just assume that the channel estimation error is included within
$\Sigma^{2}_{w}$, the variance of the Gaussian distributed noise. In (5),
$\mathbf{r}^{z}_{l}=[r^{z}_{l1},r^{z}_{l2},\dotso,r^{z}_{lM}]^{t}$ are the
recovered sensor observations in order to obtain the equivalent sequence at
the root level $l$, as, $r^{z}_{l0}[k]=x^{z}_{lm}\sum_{q}f^{z}_{lm}[k-Zq]$. It
is possible to recover the original transmitted $m$th sensor observation as,
$x^{z}_{lm}=\sum_{k}f^{z}_{lm}[k-Zq]r^{z}_{l0}[k]$. In (5),
$\hat{\mathbf{\Theta^{z}_{l}}}$ is the vector of the auto-correlation
functions of the wavelet packets at the $l$th level, expressed as,
$\hat{\mathbf{\Theta^{z}_{l}}}\triangleq[R_{\phi}(T_{0}-qT_{0}),R_{\phi}(2T_{0}-qT_{0}),\dotso,R_{\phi}(KT_{0}-qT_{0})]$,
where $R_{\phi}(\cdot)$ is the auto-correlation function of
$\phi_{00}((Zq+k)T_{0})$.

##  3 Filters and Fusion

At the DFC, we employ linear filter-based optimum and sub-optimum decision
fusion rules, by recovering decision from each sensor from each group and then
fusing them to obtain the final decision for each group of sensors. For
example, the received set of vectors at the DFC,
$\hat{Y}^{1}_{1},\hat{Y}^{2}_{1},\dotso,\hat{Y}^{Z-1}_{\lceil\log_{2}Z\rceil},\hat{Y}^{Z}_{\lceil\log_{2}Z\rceil}$
can be analyzed at the DFC to obtain a reliable set of decisions;
$\mathcal{H}^{1}_{i}$, $\mathcal{H}^{2}_{i}$, $\dotso$, $\mathcal{H}^{Z}_{i}$
where $i=1$ or 0 depending on the tested scenario.

###  3.1 Design of Fusion Rules

We start by deriving the optimum fusion rule based on the linear multi-rate
filter bank as,

$\displaystyle\ln[p(\mathbf{r}^{z}_{l}|\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q],\mathcal{H}^{z}_{i,l})\big{]}$
$\displaystyle\quad\approx\sum_{m=1}^{M}\ln[\sum_{x_{m}}\psi(r^{z}_{lm}|x^{z}_{lm})P(x^{z}_{lm}|\mathcal{H}^{z}_{i,l})\big{]}$
(6)

for parallel access sensor channels [15] where $\psi(\cdot|\cdot)$ is the
conditional distribution of individual recovered sensor decision with respect
to the transmitted one. In case of multiple access sensor channel, the log-
likelihood ratio (LLR) [10] of $\mathbf{r}^{z}_{l}$ can be approximated as,

$\displaystyle\Lambda^{z}_{l}\triangleq\ln[\frac{p(\mathbf{r}^{z}_{l}|\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q],\mathcal{H}^{z}_{1,l})}{p(\mathbf{r}^{z}_{l}|\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q],\mathcal{H}^{z}_{0,l})}\Bigg{]}$
(7)

where
$\hat{\sigma^{z}_{lm}[k,q]}=\sum_{q}\sum_{k}\hat{f}^{z}_{lm}[k]R_{\phi}(kT_{0}-qT_{0}+\Delta)$,
$R_{\phi}$ is the auto-correlation function of $\phi_{00}((Zq+k)T_{0})$,
$P(\cdot)$ and $P(\cdot|\cdot)$ are the individual and conditional probability
mass functions respectively and $p(\cdot)$ and $p(\cdot|\cdot)$ are the
individual and conditional probability density functions respectively,
$\hat{f}^{z}_{lm}[k]$ is the set of modified filter coefficients and $\Delta$
is the random variable that represents the timing discrepancy between sensor
transmissions and reception at the DFC in absence of any synchronization
between them. Since $\Delta$ varies slowly over a group of symbols, a constant
timing error is visible. Equation (7) can then be expressed as,

$\displaystyle\Lambda^{z}_{l}\approx\sum_{m=1}^{M}\ln[\frac{\psi(r^{z}_{lm}|x^{z}_{lm}=1)P^{z}_{D,lm}+\psi(r^{z}_{lm}|x^{z}_{lm}=-1)(1-P^{z}_{D,lm})}{\psi(r^{z}_{lm}|x^{z}_{lm}=1)P^{z}_{F,lm}+\psi(r^{z}_{lm}|x^{z}_{lm}=-1)(1-P^{z}_{F,lm})}\Bigg{]}$
(8)

where $P^{z}_{D,lm}$ and $P^{z}_{F,lm}$ are the probabilities of detection and
false alarm of the $m$th sensor on the $l$th layer, with the $m$th sensor
belonging to the $z$th group of sensors.

Using linear filter bank processing, we can therefore formulate a set of sub-
optimum rules such that,

$\displaystyle{\mathbf{r}}^{z}_{l}=(\mathbf{A}^{z}_{l})^{\dagger}\hat{\mathbf{Y}}^{z}_{l}$
(9)

where $\mathbf{A}^{z}_{l}$ can be realized using either Matched Filter (MF) or
a zero-forcing (ZF) detector and $(\cdot)^{\dagger}$ denotes the conjugate
transpose. In essence,

$\displaystyle{\mathbf{r}}^{z}_{l,\text{MF}}$
$\displaystyle\triangleq(\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q])^{\dagger}\hat{\mathbf{Y}}^{z}_{l}$
$\displaystyle{\mathbf{r}}^{z}_{l,\text{ZF}}$
$\displaystyle\triangleq(\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q])^{\dagger}\big{(}(\mathbf{\mathcal{D}}^{z}_{l})^{-1}\big{)}^{\dagger}\hat{\mathbf{Y}}^{z}_{l}$
(10)

where
$\mathbf{\mathcal{D}}^{z}_{l}=\frac{1}{N}(\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q])^{\dagger}(\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q])$
is a diagonal matrix for $N>>M$.

###  3.2 Performance Measures

If we assume that
$P(\hat{\mathbf{Y}}^{z}_{l}|\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q],\mathcal{H}^{z}_{i,l})$
is Gaussian mixture distributed, then we can express
$\mathbf{r}^{z}_{l}|\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q],\mathcal{H}^{z}_{i,l}$
as Gaussian mixture distributed such that,
$\mathbf{r}^{z}_{l}|\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q],\mathbf{x}^{z}_{l}\sim\mathcal{N}\big{(}\mathbb{E}\\{\mathbf{r}^{z}_{l}|\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q],\mathbf{x}^{z}_{l}\\},\mathbb{V}\\{\mathbf{r}^{z}_{l}|\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q],\mathbf{x}^{z}_{l}\\}\big{)}$
where $\mathcal{N}(\mathbb{E}\\{\cdot\\},\mathbb{V}\\{\cdot\\})$ represents
normal distribution with mean $\mathbb{E}\\{\cdot\\}$ and covariance
$\mathbb{V}\\{\cdot\\}$. Therefore,

$\displaystyle{\mathbf{r}}^{z}_{l,\text{MF}}|\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q],\mathbf{x}^{z}_{l}$
$\displaystyle\sim\mathcal{N}\big{(}N\mathbf{\mathcal{D}}^{z}_{l}\sqrt{\rho_{l}}\mathbf{x}^{z}_{l},\Sigma^{2}_{e}||N\mathbf{\mathcal{D}}^{z}_{l}||\big{)}$
$\displaystyle{\mathbf{r}}^{z}_{l,\text{ZF}}|\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q],\mathbf{x}^{z}_{l}$
$\displaystyle\sim\mathcal{N}\big{(}N\sqrt{\rho_{l}}\mathbf{x}^{z}_{l},\Sigma^{2}_{e}||N(\mathbf{\mathcal{D}}^{z}_{l})^{-1}||\big{)}$
(11)

where $\mathbf{x}^{z}_{l}=[x^{z}_{l1},x^{z}_{l2},\dotso,x^{z}_{lM}]$ is the
transmit vector of $M$ sensor observations multiplexed on the $l$th level of
the wavelet packet binary tree. Now since,
$p(\mathbf{r}^{z}_{l}|\hat{\mathbf{G}}^{z}_{l}\hat{\sigma}^{z}_{lm}[k,q],\mathbf{x}^{z}_{l})\approx\prod_{m=1}^{M}\psi(r^{z}_{lm}|x^{z}_{lm})$,
we can express (3.2) in terms of individual sensor observations assuming they
are independent of each other as,

$\displaystyle\psi(r^{z}_{lm,MF}|x^{z}_{lm})$
$\displaystyle\triangleq\mathcal{N}\big{(}Nd^{z}_{lm}\sqrt{\rho_{l}}x^{z}_{lm},\Sigma^{2}_{e}Nd^{z}_{lm}\big{)}$
$\displaystyle\psi(r^{z}_{lm,ZF}|x^{z}_{lm})$
$\displaystyle\triangleq\mathcal{N}\big{(}N\sqrt{\rho_{l}}x^{z}_{lm},N\Sigma^{2}_{e}/d^{z}_{lm}\big{)}$
(12)

where $d^{z}_{lm}$ is the $m$th element of $\mathbf{\mathcal{D}}^{z}_{l}$.
From (3.2), it will be possible to estimate the transmit observations which
are BPSK modulated.

The probability of error in estimating the transmitted sensor observations can
be calculated by finding the probability of error in the final obtained
decision on whether a target or phenomenon is present ($P^{z}_{F1,l}$) or
absent ($P^{z}_{F0,l}$) at the $l$th level for the $z$th sensor group.
Therefore,

$\displaystyle
P^{z}_{F0,l}=\lim_{N\to\infty}\mathcal{Q}\Bigg{(}\frac{{\mathbf{r}}^{z}_{l}-\sqrt{N\mathbf{\mathcal{D}}^{z}_{l}\rho_{l}}\mathbf{x}^{z}_{l}/\Sigma_{e}}{\sqrt{1/2(M(1-P^{z}_{F,l})+\Sigma^{2}_{e})}}\Bigg{)}$
(13)

where $\mathcal{Q}(\cdot)$ is used to denote the complementary cumulative
distribution function (CCDF) and $P^{z}_{F,l}$ is the probability of wrongly
detecting a ‘0’ or absence of a phenomenon or a target.

##  4 Performance Analysis

Before simulating the performance of our designed WPT aided industrial WSN, we
need to decide on the design criteria of three components; noise model for
$\mathbf{e}^{z}_{l}$, the filter coefficients $f^{z}_{lm}[q]$ and the scaling
functions for the wavelets $\phi_{lz,m}$.

####  4.0.1 Noise Model

The combination of Gaussian and impulsive noise follows the pdf,
$p(\mathbf{e})=(1-\wp)\mathcal{N}(0,\Sigma_{g}^{2})+\wp\mathcal{N}(0,\sqrt{\kappa}\Sigma_{I}^{2}/A)$,
where $\wp$ denotes the probability of having an impulsive noise with zero
mean and variance $\sqrt{\kappa}\Sigma_{I}^{2}/A$, while the probability of
having a Gaussian distributed background noise is $(1-\wp)$ with zero mean and
variance $\Sigma_{g}^{2}$. We also define a factor
$\Gamma=\Sigma_{g}^{2}/\Sigma_{I}^{2}=0.25$ implying that the impulsive noise
is 25 times stronger than the background Gaussian noise. The variance of the
impulse noise is given by $\sqrt{\kappa}\Sigma_{I}^{2}/A$ where $\kappa$ is
the Poisson distributed sequence whose pdf is characterized by the impulse
index $A$. For our case, we choose $A=0.1$. For Bernoulli-Gaussian distributed
noise, the combined noise pdf can be modified to
$p(\mathbf{e})=(1-\wp)\mathcal{N}(0,\Sigma_{g}^{2})+\wp\mathcal{B}(\varrho,\varsigma\varrho(1-\varrho))$
where $\mathcal{B}$ represents Bernoulli distribution with mean $\varrho$,
variance $\varrho(1-\varrho)$ and $\varsigma$ is the frequency of occurrence
of the impulse noise. For our case, we have $\varrho=0.3$.

Figure 2: Comparative ROC ($P_{D_{0}}$ v/s $P_{F_{0}}$) of WPDM and DF aided
industrial sensing network with different wavelet scaling functions and fusion
rules in presence of Middleton Class A impulse noise with parameters
$Z=4,M=8,N=64,\wp=0.3$ at a fixed SNR of 10dB. Figure 3: Comparative ROC
($P_{D_{0}}$ v/s $P_{F_{0}}$) of WPDM and DF aided industrial sensing network
with different wavelet scaling functions and fusion rules in presence of
Middleton Class A impulse noise with parameters $Z=4,M=8,N=64,\wp=0.5$ at a
fixed SNR of 10dB.

####  4.0.2 Filter Coefficients

For our synthetic system, we choose an FIR filter of length $Q=14$ and $K=2$.
The filter bandwidth $B$ is chosen to be equal to $B=\sqrt{2}$ such that
$K_{0}=2K-1-\lceil 2\log_{2}B\rceil\geq 1$ where $K_{0}$ is the number of
continuous derivatives of the auto-correlation functions, $R_{\phi}(\cdot)$.
Using Nyquist sampling rate, the sampling frequency $f_{s}$, can be calculated
as $B=f_{s}/2$. Using these information, we calculated the filter
coefficients, $h[q]=\operatorname{sinc}(q-D)$ and $g[q]=(-1)^{q}h[2k+1-q]$
where $D$ is the total delay of the FIR filters given by $D=(Q-1)/4B$, $h[q]$
and $g[q]$ are the filter coefficients for the two branches of the binary tree
structure at each level of the tree (refer to Fig. 1). The filter coefficients
$h[q]$ and $g[q]$ for $q=0,1,\dotso,13$ satisfy the orthonormality and
regularity constraints. The equivalent filter coefficients $f^{z}_{lm}[q]$ can
then be built from $h^{z}_{lm}[q]$ and $g^{z}_{lm}[q]$, from which the filter
sequence $\mathbf{f}^{z}_{l}[q]$ can be built for $M$ sensors belonging to the
$z$th group of sensors.

####  4.0.3 Scaling Functions

Three different scaling functions are constructed using standard wavelets from
[13]. Performance of our designed transmission system is simulated using the
three scaling functions separately, comparing them to find out the set of
wavelet scaling functions that offers the best fusion performance under
different conditions. The wavelet scaling functions selected are,

$\displaystyle\text{Haar :}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \phi_{00}(x)=1\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ 0\leq x<1$
$\displaystyle\text{Shannon :}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \phi_{00}(x)=\operatorname{sinc}(x)$
$\displaystyle\text{Piecewise linear spline :}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\phi_{00}(x)=1-|x|$ (14)

where $\phi_{00}(x)$ are the root scaling functions which are fed back to (1)
to generate the family of scaling functions to obtain,

$\displaystyle\text{Haar :}\leavevmode\nobreak\
\mathbf{\Phi}^{z}_{l}=\sum_{q=1}^{Q}\mathbf{f}^{z}_{l}[q];\leavevmode\nobreak\
\text{Shannon :}\leavevmode\nobreak\
\mathbf{\Phi}^{z}_{l}=\sum_{q=1}^{Q}\mathbf{f}^{z}_{l}[q]\leavevmode\nobreak\
\operatorname{sinc}(qT_{0})$ $\displaystyle\text{Piecewise linear spline
:}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\mathbf{\Phi}^{z}_{l}=\sum_{q=1}^{Q}\mathbf{f}^{z}_{l}[q](1-|qT_{0}|)$ (15)

Using the scaling functions from (4.0.3), we simulate the performance of a WPT
aided industrial WSN. Each of the $Z$ sensor groups are randomly deployed and
uniformly distributed in a circular annulus around the DFC with radii
$\varphi_{\text{max}}=1000$ m and $\varphi_{\text{min}}=100$ m. Within each
sensor cluster, $M$ sensors are randomly deployed within a rectangular area
(size of 50m $\times$ 100m) following a half-normal distribution [12]. The
channel coefficients are log-normal shadowed and Rayleigh faded such that,
$h^{z}_{n,lm}\sim\mathcal{N}_{\mathbb{C}}(0,\text{diag}(B^{z}_{lm}))$ where
$\lambda^{z}_{lm}=\upsilon_{lm}\big{(}\frac{\varphi_{\text{min}}}{\varphi_{lm}}\big{)}^{\eta}$,
$B^{z}_{lm}=\big{(}\beta^{z}_{lm}(0),\dotso,\beta^{z}_{lm}(qT_{0})\big{)}^{t}$,
$10\log_{10}(\upsilon_{lm})\sim\mathcal{N}\big{(}\mu_{\lambda}dB,\sigma^{2}_{\lambda}dB\big{)}$,
$\eta$ is the pathloss exponent, $\varphi_{lm}$ is the distance of the mth
sensor from the DFC at the $l$th level of the binary wavelet packet tree,
$\rho_{l}=1/\sqrt{N}$ and $\\{P^{z}_{D,l},P^{z}_{F,l}\\}=\\{0.5,0.05\\}$.

Figure 4: Comparative ROC ($P_{D_{0}}$ v/s $P_{F_{0}}$) of WPDM and DF aided
industrial sensing network with different wavelet scaling functions and fusion
rules in presence of Middleton Class A impulse noise with parameters
$Z=4,M=8,N=64,\wp=0.7$ at a fixed SNR of 10dB.

Figure 5: Comparative probability of erroneous/false detection as a function
of SNR (dB) for ZF-based fusion applied to a WPDM-aided ($Z=4,M=8,N=64$)
industrial sensor network with different wavelet scaling functions when the
impulse noise ($\wp=0.3$) is either Middleton Class A or Bernoulli-Gaussian
distributed.

Figure 6: Comparative probability of erroneous/false detection as a function
of SNR (dB) for ZF-based fusion applied to a WPDM-aided ($Z=4,M=8,N=64$)
industrial sensor network with different wavelet scaling functions when the
impulse noise ($\wp=0.5$) is either Middleton Class A or Bernoulli-Gaussian
distributed.

Figure 7: Comparative probability of erroneous/false detection as a function
of SNR (dB) for ZF-based fusion applied to a WPDM-aided ($Z=4,M=8,N=64$)
industrial sensor network with different wavelet scaling functions when the
impulse noise ($\wp=0.7$) is either Middleton Class A or Bernoulli-Gaussian
distributed.

We generate two sets of figures, (Fig. 2, Fig. 3, Fig. 4) demonstrating
receiver operating characteristics (ROC) and (Fig. 7, Fig. 7, Fig. 7)
depicting the probability of false detection against signal-to-noise ratio
(SNR) in dB. Both the set of figures are generated for the configurations,
$Z=4,M=8,N=64$, thereby representing a virtual large-scale distributed MIMO
set-up. In this case, we are considering 4 groups of sensors, each group
containing 8 sensors, while the DFC is equipped with 64 antennas. Therefore,
the number of levels of the binary tree, $L=\log_{2}Z=\log_{2}4=2$. Also the
$T_{0}=10^{-3}$ for the interval between BPSK modulated symbols that are
transmitted as the sensor observations; hence
$T_{l}=\\{2^{1}T_{0},2^{2}T_{0}\\}=\\{2\times 10^{-3},4\times 10^{-3}\\}$.
Also the pathloss exponent for the MAC channel is 2, with
$(\mu_{\lambda},\sigma_{\lambda})=(4,2)$ dB representing indoor smart
industrial environment [16].

All the figures represent the MAC scenario suffering from different
combinations of Gaussian and impulsive noise by varying $\wp=(0.3,0.5,0.7)$.
The performance of the system with $Z=4,M=8,N=64$ where the sensor
observations are sent over the MAC without being encoded into wavelet packets
with maximal ratio combining (MRC) for fusing the observations at the DFC, is
plotted as the benchmark in Fig. 2, Fig. 3, Fig. 4. The benchmark curves are
labelled as ‘w/o WPDM+DF’. For ROCs in Fig. 2, Fig. 3 and Fig. 4, three sets
of scaling functions are compared, Haar, Shannon and Piecewise linear spline
when two different fusion rules are used, MF and ZF. In all cases, ZF offers
better fusion performance as compared to MF. The SNR is fixed at 10 dB. Fig.
2, $\wp=0.3$ depicts a condition where the impulse noise is low with higher
background Gaussian noise. In this case, all the three wavelet functions
perform very close to each other and in such a condition, the Haar scaling
function can be preferred over the others owing to its low computational
complexity. In Fig. 3, $\wp=0.5$, presents a condition with equal amount of
Gaussian and impulse noise. Piecewise linear spline function offers the best
performance. The reason can be attributed to the fact that spline wavelets are
formed by linear combination of B-splines [13] and inherit the properties of
the basis functions that form the wavelets. Therefore, they are suitable for
isolating noise from the original signal depending on the application
scenario. In Fig. 4, impulse noise is higher than the Gaussian noise
($\wp=0.7$). In this case, the Shannon wavelet performs best, as with complex-
valued wavelets based on $\operatorname{sinc}$, it is highly successful in
denoising heterogeneous sensor observation signals.

For the plots of probability of erroneous/false detection, we compare
performances of three wavelet scaling functions with the case without WPDM.
For the DFC side, we just use ZF-based decision fusion. WPDM is capable of
denoising the sensor signals and offers improvement in performance by more
than 10 times in the worst-case high impulse noise scenario ($\wp=0.7$, Fig.
7) for a particular value of SNR. Transmission of sensor observations aided by
WPDM offers a large improvement in performance in a noisy environment. This is
because the transmitted waveforms overlap in time, thereby dispersing the
impulse noise energy over several symbols at each binary tree terminal. A
moderate noise burst that is capable of resulting in erroneous detection is
distributed over several waveforms without causing any false detection.

##  5 Conclusion

In order to design a robust large-scale sensing network for the emerging smart
Industry/Industry 4.0 environments, we conceive the novel design of combining
WPDM of sensor decisions at the transmit side with DF of heterogeneous sensor
decisions on the receive side. Evaluating performance of our proposed system
in presence of Gaussian and impulsive noise, Rayleigh block fading and indoor
shadowing, we demonstrated that WPDM-aided industrial WSN outperforms
conventional large-scale multi-antenna WSN arrangements substantially, in
presence of impulsive noise particularly. However, in this paper we assume
that each sensor group, $z$, consists of an equal number of sensors, $M$. In
future, we will relax this assumption to include different number of sensors
belonging to each sensor group. We also plan to conduct a detailed study of
the wavelet scaling functions and how to optimize their design depending on
the environment at hand.

## References

  * [1] What is Industry 4.0—the Industrial Internet of Things (IIoT), Available from: https://www.epicor.com/en-in/resource-center/articles/what-is-industry-4-0/.
  * [2] M. Wollschlaeger, T. Sauter, and J. Jasperneite, “The future of industrial communication: automation networks in the era of the internet of things and industry 4.0", _IEEE Industrial electronics magazine_ , vol. 11, no. 1, pp. 17–27, 2017.
  * [3] B. Afzal, M. Umair, G. A. Shah and E. Ahmed, “Enabling IoT platforms for social IoT applications: Vision, feature mapping, and challenges." _Future Gener. Comput. Syst._ , vol. 92, 2019, pp. 718–731.
  * [4] D. Ciuonzo, G. Romano, and P. Salvo Rossi, “Channel-aware decision fusion in distributed MIMO wireless sensor networks: Decode-and-fuse vs. decode-then-fuse,” _IEEE Trans. Wireless Commun._ , vol. 11, no. 8, pp. 2976–2985, Aug. 2012.
  * [5] A. Chawla, A. Patel, A. K. Jagannatham, and P. K. Varshney, “Distributed detection in massive MIMO wireless sensor networks under perfect and imperfect CSI," _IEEE Trans. Signal Processing_ , vol. 67, no. 15, pp. 4055–4068, 2019.
  * [6] M. Cheffena, “Propagation channel characteristics of industrial wireless sensor networks", _IEEE Antennas Propag. Mag._ , vol. 58, 2016, pp. 66–73.
  * [7] A. R. Lindsey, “Wavelet packet modulation for orthogonally multiplexed communication," _IEEE Trans. Signal Processing_ , vol. 45, pp. 1336–1339, May 1997.
  * [8] K. M. Wong, J. Wu, T. N. Davidson, and Q. Jin, “Wavelet packet division multiplexing and wavelet packet design under timing error effects," _IEEE Trans. Signal Processing_ , vol. 45, pp. 2877–2890, Dec. 1997.
  * [9] I. Dey, P. Salvo Rossi, M. M. Butt, and N. Marchetti, “Virtual MIMO wireless sensor networks: Propagation measurements and fusion perfor- mance," _IEEE Trans. Antennas and Prop._ , vol. 67, no. 8, pp. 5555–5568, 2019\.
  * [10] D. Ciuonzo, P. Salvo Rossi, and S. Dey, “Massive MIMO channel-aware decision fusion," _IEEE Trans. Signal Process._ , vol. 63, no. 3, pp. 604– 619, 2015.
  * [11] D. Middleton, “Statistical-physical models of electromagnetic interference," _IEEE Trans. Electromag. Compat._ , vol. EMC-19, pp. 106–127, Aug. 1977.
  * [12] A. Papoulis, _Probability, Random Variables, and Stochastic Processes_ , 3rd ed: McGraw-Hill, 1991.
  * [13] I. Daubechies, _Ten Lectures on Wavelets._ Philadelphia, PA: SIAM, 1992.
  * [14] A. N. Akansu and R. A. Haddad, _Multiresolution Signal Decomposition. Transforms, Subbands, and Wavelets._ Boston, MA: Academic, 1992.
  * [15] A. Lei and R. Schober, “Coherent Max-Log decision fusion in wireless sensor networks," _IEEE Trans. Commun._ , vol. 58, no. 5, pp. 1327–1332, May 2010.
  * [16] I. Dey, P. Salvo Rossi, M. M. Butt, and N. Marchetti, “Virtual MIMO wireless sensor networks: Propagation measurements and fusion performance," _IEEE Trans. Antennas and Prop._ , vol. 67, no. 8, pp. 5555–5568, 2019\.


11institutetext:  Algorithm Engineering, Institute of Theoretical Informatics,
KIT-Department of Informatics, Karlsruhe Institute of Technology (KIT),
Germany

#  SATViz: Real-Time Visualization of Clausal Proofs

Tim Holzenkamp  Kevin Kuryshev  Thomas Oltmann  Lucas Wäldele  Johann Zuber
Tobias Heuer Project Supervision Markus Iser Project Definition and
Supervision

###### Abstract

Visual layouts of graphs representing SAT instances can highlight the
community structure of SAT instances. The community structure of SAT instances
has been associated with both instance hardness and known clause quality
heuristics. Our tool SATViz visualizes CNF formulas using the variable
interaction graph and a force-directed layout algorithm. With SATViz, clause
proofs can be animated to continuously highlight variables that occur in a
moving window of recently learned clauses. If needed, SATViz can also create
new layouts of the variable interaction graph with the adjusted edge weights.
In this paper, we describe the structure and feature set of SATViz. We also
present some interesting visualizations created with SATViz.

##  0.1 Introduction

Visual representations of algorithms can improve our understanding of
algorithmic concepts and the nature of the underlying problems. Visualizations
of the Conflict-Driven Clause Learning algorithm (CDCL) are of high interest,
as this is currently the most successful algorithm for solving the
propositional satisfiability problem (SAT problem). Despite the complexity of
the SAT problem, implementations of CDCL in so-called SAT solvers are
successfully used in industrial practice, e.g., in software verification [4],
hardware verification [12], product configuration [21], or planning [18].

For satisfiable instances, SAT solvers output a variable assignment which
serves as a certificate of satisfiability. Such a certificate can be checked
in polynomial time and hence SAT $\in$ NP. However, modern SAT solvers also
output certificates of unsatisfiability – a proof given in the DRAT format
[8]. Such a proof can be checked efficiently by a procedure with runtime
polynomial in the proof length. The output of verifiable proofs increases the
trust in the correctness of SAT algorithms. Trust is crucial in scenarios
where SAT solvers, e.g., verify the absence of bugs in safety-critical
software.

Famous proofs of unsatisfiability that solved previously open mathematical
problems, e.g., the Boolean Pythagorean triples problem, have been generated
by SAT solvers [9, 10]. The explainability of such automatically generated
proofs is of great interest. This is mainly due to their sheer size, e.g., the
total size of the proof of the Boolean Pythagorean triples problem is 200 TB.

The source of such generated proofs are the clauses learned by a CDCL SAT
solver. In CDCL, learned clauses, i.e., small proofs, also help to accelerate
search-space navigation. The structure of the resulting reasoning is hard to
grasp as millions of clauses are involved. The average instance of the SAT
competition 2021 benchmark contains $2.3$ millions of clauses. When running
Kissat [2] on these instances we observe an average of 8676 (with a median of
4040) learned clauses per second.

To this end, visual representations illustrating learned clauses during the
execution of SAT solvers can bridge the gap to gain a better understanding of
the structure of automated reasoning and resolution-space navigation. In this
paper, we present a tool which creates real-time visualizations of proof-
sequences as they are created by a given SAT solver or proof file.

##  0.2 Background

Instances of the SAT problem are given in _Conjunctive Normal Form_ (CNF) and
are represented by a _set of clauses_. Each clause is a is a _set of literals_
and represents the disjunction of those literals. Each literal is a negated or
non-negated Boolean variable.

CDCL combines search (cf. DPLL algorithm) and resolution (cf. saturation
algorithm) [3]. CDCL solvers maintain and extend a current partial assignment
which is created from repeated decision and subsequent Boolean constraint
propagation. On conflict, CDCL infers a new clause by resolution of the
conflict-reason clauses which are determined by analyzing the implication
graph of the current partial assignment [7]. If the empty clauses can be
inferred, then the set of learned clauses forms a testable proof of the
unsatisfiability of an instance. Clause learning is also an important means of
narrowing down the search space of satisfiable instances, e.g., to escape from
unsatisfiable regions of the search space early on. As the relevance of a
learned clause for the overall search process can not be known in advance,
clause learning and forgetting is controlled by heuristics [15].

A SAT formula in conjunctive normal form (CNF) can be represented as a
hypergraph. Hypergraphs are generalizations of regular graphs in which each
_hyperedge_ can connect more than two nodes. In this model, each variable
corresponds to a node and each clause is modeled as a hyperedge. While there
is little work on visualizing hypergraphs, we transform the hypergraph to a
regular graph. The variant in which each hyperedge is replaced by a clique is
known as the _Variable Interaction Graph_ (VIG) in the SAT domain (cf. Section
0.4.2).

##  0.3 Related Work

| DPvis | 3Dvis | iSAT | SATGraf | SATViz
---|---|---|---|---|---
Interactive | ✓ | ✗ | ✓ | ✓ | ✓
Multiple Reduction Types | ✗ | ✗ | ✓ | ✓ | ✓
3D Layout | ✗ | ✓ | ✗ | ✗ | ✗
Real-time Animations | ✓ | ✗ | ✗ | ✓ | ✓
Graph Contractions | ✗ | ✗ | ✗ | ✗ | ✓
IPASIR Interface | ✗ | ✗ | ✗ | ✗ | ✓
DRAT Interface | ✗ | ✗ | ✗ | ✗ | ✓
Table 1: Feature Coverage of CNF Formula Visualizers

In the area of SAT solving, the tools DPvis and 3Dvis by Sinz have been used
to visualize the structure of SAT instances based on the force-directed layout
of their graph representations [20]. DPvis can also visualize formula
evolution by simplifying according to the branching and propagation steps in
runs of the integrated SAT solver Minisat [6]. Visualizations of 3Dvis appear
in Knuth’s “The Art of Computer Programming” Vol. 4.6 on “Satisfiability” (the
author’s favorite pages) [13]. Recent SAT competitions use the intriguing
visualizations of 3Dvis as their
logos.111https://satcompetition.github.io/2022/logo2022-large.png

iSAT is a tool for instrumentation and interactive control of SAT solvers.
iSAT facilitates the analysis of intermediate solver states. For external
visualization, iSAT can export several types of graph representations of
formulas to files [16]. The tool SATGraf by Newsham et al. visualizes
community structure of SAT formulas. SATGraf can also display statistics and
solving progress of specially instrumented SAT solvers running on small SAT
instances [14].

Table 1 shows a comparison of the features covered by SATViz and related
tools. Except for the more performance-critical 3Dvis, all tools considered
offer some form of interactive control to reinitialize the visualization with
the current solver state. While DPvis and 3Dvis are limited to the use of the
VIG, the other tools offer additional options and some form of extensibility.
3Dvis is the only tool that can generate three-dimensional layouts of SAT
instances.

Real-time animations are supported by DPvis, SATGraf and SATViz, but each of
them highlights a different aspect. DPvis focuses on search progress via unit
resolution, SATGraf focuses on the survival of the pre-computed communities,
and SATViz highlights recently learned clauses. Both SATGraf and SATViz can
relayout an instance which has been modified by learned and forgotten clauses.

However, SATViz is the only tool that facilitates the handling of large
industrial SAT instances in real-time. This is possible since SATViz offers
graph contractions for reducing the number of edges while preserving the graph
structure. Like this, even large SAT instances can be visualized efficiently
and also become less cluttered.

SATViz instruments SAT solvers with standard methods of the IPASIR [1] or can
moreover replay proofs given in the DRAT format [8]. In contrast to other
tools, SATViz can visualize the course of any SAT solver which implements the
IPASIR or exports DRAT proofs.

##  0.4 The Tool SATViz

With our tool SATViz, we can transform CNF formulas to graphs and visualize
them using a graph drawing algorithm. A SAT solver running on the local or a
remote machine can then connect and send learned clauses to the application,
which are then highlighted in the visual graph representation. Thereby, the
color of a node represents how frequently the corresponding variable appeared
in learned clauses (similar to a heat-map). Through this visualization we can
better analyze the working of CDCL algorithms and better understand the
structure of different problem instances. Code and instructions for using
SATViz can be found on GitHub.222https://github.com/satviz/satviz

###  0.4.1 Architecture

Clause Consumer Clause Producer $\oplus$ IPASIR Solver Clausal Proof CNF
Formula TCP/IP Graph
Transformer

Figure 1: The Architecture of SATViz

Figure 1 illustrates the core components of SATViz’ architecture. It can be
divided into two central components: the _clause consumer_ and the _clause
producer_. The clause consumer is responsible for visualizing SAT instances
and learned clauses, while the clause producer sends a stream of learned
clauses to the clause consumer over the network. This decoupling adds some
flexibility to the architecture, as learned clauses can be sent from different
machines from any software implementing the clause producer protocol.

The clause producer of SATViz can read clauses from two types of sources: a
file containing a clausal proof or a shared library of an SAT solver
implementing the IPASIR interface [1]. In the former, SATViz reads a clausal
proof, e.g., in DRAT file format [8], and streams it to the clause consumer.
In the latter, the clause producer starts a SAT solver initialized with a
given CNF formula and streams learned clauses to the clause consumer at the
time they are exposed over the IPASIR interface (intercepted in an
`ipasir_set_delete` callback function).

The clause consumer buffers the incoming stream of learned clauses and sends
chunks of incoming clauses with a fixed frame rate to the _graph transformer_
component. The graph transformer is responsible for visualizing learned
clauses in the graph representation. It adds new edges or updates their
weight, and changes the colors of the corresponding nodes contained in the
learned clauses. This transformation is configurable and extensible. The graph
transformer configuration includes parameters to control the heat-map which
highlights learned clauses (cf. Section 0.4.4), different graph models to
represent CNF formulas (cf. Section 0.4.2), as well as options for reducing
the size of larger instances to improve the readability of the visualization
and to speedup the performance of the animation (cf. Section 0.4.3).

Given a graph, the layout algorithm uses force-directed placement to create
the visual representation of the graph. We use the implementation provided by
the Open Graph Drawing Framework (OGDF) [5]. The heat-map algorithm
continuously updates node colors according to the incoming learned clauses.
The placement algorithm creates and uses an initial layout based on the
original formula. But on demand, it can recalculate a new layout based on the
modified graph induced by the clauses which have been learned and forgotten so
far. The user interface also allows to pause, stop, rewind and replay the
solver run and to step in at any time into the proof. For larger graph the
magnification feature using the mouse wheel is very practical.

SATViz also supports to record the visualization including user interactions
and exports this to a video file. SATViz can also operate in a headless mode
and silently create videos from given proofs or SAT solver runs in the
background. Some sample videos can be found in the dedicated YouTube playlist
of
Kuryshev.333https://www.youtube.com/playlist?list=PLhFsx92qubLxMsFQ4lUlCS6vFmJvTbYJS

###  0.4.2 Hypergraph Reduction

A CNF formula induces a hypergraph such that the variables correspond to nodes
and each clause forms a hyperedge spanning its variables. Hypergraphs are hard
to visualize directly, which is why reductions to graph representations are
used in practice [20]. The two most common representations are the _clique
expansion_ and the _bipartite graph representation_ [11, 19]. The former adds
a clique between all nodes contained in a hyperedge, while the latter
additionally models the hyperedges as nodes and connects each with their
corresponding nodes. However, the clique expansion unnecessarily increases the
size of our graph model in presence of large hyperedges, while the bipartite
graph representation adds much more nodes than necessary (the number of
clauses is usually considerably larger than the number of variables).

We therefore present a third reduction to which we refer to as the _ring
reduction_. Here, we sort the node IDs of each hyperedge in increasing order
and connect consecutive nodes with an edge. Additionally, we add an edge
between the two nodes with the lowest and highest ID in each hyperedge. This
representation combines the advantages of both standard representations as the
number of nodes equals to the nodes in clique expansion, while the number of
edges equals to the edges in the bipartite graph representation.

In SATViz, we support the ring reduction and clique expansion as graph models.
Additionally, each hyperedge contributes an edge weight which is a
configurable function inversely proportional to its own size.

###  0.4.3 Graph Contraction

Figure 2: Visualization of instance AProVE09-15 (94663 variables) highlighting
clauses learned by Kissat with no contraction (left) and 10 iterations of max-
weight contractions (right).

The size of SAT instances have a large impact on visualization performance. In
initial experiments, the time to layout SAT instances with more than 50k
variables reached a for usability critical time limit of one minute. Moreover,
visualizations of larger SAT instances can become cluttered and their
structure concealed as too many nodes compete for space. For these reasons,
SATViz features a pre-processing step in which it runs a graph contraction
algorithm which is based on the well-known label propagation algorithm [17].
The algorithm assigns labels to the nodes. Initially, each node has its own
label. Then, the algorithm iterates over the nodes in random order – also
called a _round_ – and whenever a node is visited, it is assigned the label
appearing most frequently in its neighborhood. This is repeated for a fixed
number of rounds or until none of the nodes changed its label in a round.
Subsequently, we collapse nodes with the same label into a single node and
connect them to nodes corresponding to adjacent labels. We repeat the
clustering and contraction procedure recursively until the size of the graph
becomes manageable for our graph drawing algorithm.

Figure 2 shows two visualizations of the SAT instance
AProVE09-15444https://gbd.iti.kit.edu/file/8a244542d09e20e9e8813ce7089c4135
which has 94663 variables. The left visualization represents the input
instance, while the right visualization was obtained with our graph
contraction algorithm. As it can be seen, the graph contractions facilitate
the analysis of large SAT instances and can unveil hidden structures.

###  0.4.4 Highlighting Areas of Interest

Figure 3: Visualization of instance AProVE09-15 highlighting clauses learned
by Kissat – with a heat-map frame width of 1000 clauses – after 1000, 1690,
and 3090 conflicts (from left to right).

SATViz uses configurable heat-map colors to highlight recently learned
clauses. We use a _heat value_ to determine the color for each node. This
value is determined by either (i) normalizing each variables’ occurrence count
for the last $k$ learned clauses or (ii) assigning the maximum heat value to
the most recently used variables and then reduce the value to zero again
within the next $k$ time steps (unless it appears again in a learned clause).
The parameter $k$ is configurable and assists in switching between coarse- and
fine-grained proof analysis. The color of node contractions is determined by
the average heat value of the variables which it represents. SATViz also has a
parameter for adjusting the speed of the animation.

Figure 3 shows three snapshots of a run of Kissat [2] on the instance
AProVE09-15 within the first 5000 learned clauses. The heat-map highlights
variables occurring in the $k=1000$ most recently learned clauses. The
displayed sequence of resolution steps appears in the animation like a
windshield wiper.

###  0.4.5 Evolution of a Proof via Relayouting

Figure 4: Visualization of instance AProVE09-15 solved by Kissat. The right
graph shows the instance after a relayout using adapted weights after 6000
learned or forgotten clauses. The left graph shows the magnified center of the
circular region of the right graph.

SATViz offers to pause the animation and relayout the instance with the
weights adjusted by the learned and deleted clauses. Figure 4 shows the
structural changes in the AProVE09-15 instance after 6000 clauses learned by
Kissat. We observe several variable connections become loose and degenerate to
simple structures, which means that many original clauses become redundant and
get deleted within the first 6000 conflicts. However, some variable
connections get stronger, as the graph forms a densely connected community at
the center of the circular structure of the relayouted instance in Figure 4.

###  0.4.6 Observing the Evolution of Formula Structure

Figure 5 shows the evolution of the SAT instance
Newton.5.1.i.smt2-cvc4555https://gbd.iti.kit.edu/file/92b98ab8055b143e0283a215a70cf001
as it was solved by Kissat. We observe that the original structure decays and
a new structure emerges in which several nodes form circular shapes and
tentacles. However, the resulting tentacles and circles are arranged around a
densely connected core of nodes. This is due to several variables becoming
weakly connected or even disconnected from the densely connected instance
core. This densely connected core is where further searching and learning
takes place, while nothing happens to the nodes that are now loosely
connected. As learning progresses, more and more nodes detach from this core,
so that it thins out and becomes flatter, while the inner circles around this
core become more densely populated.

Figure 5: Visualization of instance Newton.5.1.i.smt2-cvc4 solved by Kissat.
The first row shows the instance in its original layout (left) and after 5000
learned or forgotten clauses (right). The second row shows the instance layout
after one, two, and three million learned or forgotten clauses (from left to
right). The third row shows the evolution of the magnified centers of the
updated layouts after one, and three million learned or forgotten clauses and
how it looks close to the of search at about 3.5 million learned or forgotten
clauses (from left to right).

##  0.5 Future Work and Conclusion

Little research has gone into incremental weight updates in a force-directed
layout algorithm. We conducted initial experiments with real-time incremental
adjustment of force-directed node placement. However, it was problematic to
regulate the extend of which node positions are allowed to change in such
incremental updates. Developing a stable incremental placement algorithm is a
research project on its own.

Our future work will mainly focus on the consolidation and extension of SATViz
in the following respects. SATViz should be able to process additional input
about _named equivalence classes of variables_. Such extra information can be
displayed, e.g., by using additional node colors or when hovering nodes with
the mouse. Equivalence classes of variables can also be used to specify
preferences for the contraction algorithm or to increase the node proximity
for such related variables with an adapted layout algorithm. With the
possibility to process and use classes of variables, SATViz can be even more
helpful in analyzing large proofs and in identifying patterns and stages in
proof generation.

We also want to increase the number of available graph reduction methods and
parameters to control them. As SATViz decouples clause producers via sockets,
SATViz can also be extended to monitor clause exchange of parallel SAT solvers
in the cloud.

Visually observing the effects of the methods and heuristics that control
clause learning and forgetting, or even clause sharing, can deepen our
understanding of these methods and how they interact. The visualizations of
SATViz can even leverage our understanding of the coarse structure of clausal
proofs. SATViz can help identify and analyze arguments that consist of a large
number of clauses. Recognizing patterns in clause proofs can help us better
understand, compress, and explain them.

## References

  * [1] Tomás Balyo, Armin Biere, Markus Iser, and Carsten Sinz.  SAT race 2015.  Artif. Intell., 241:45–65, 2016.
  * [2] Armin Biere, Katalin Fazekas, Mathias Fleury, and Maximillian Heisinger.  CaDiCaL, Kissat, Paracooba, Plingeling and Treengeling entering the SAT Competition 2020.  In Tomas Balyo, Nils Froleyks, Marijn Heule, Markus Iser, Matti Järvisalo, and Martin Suda, editors, Proc. of SAT Competition 2020 – Solver and Benchmark Descriptions, volume B-2020-1 of Department of Computer Science Report Series B, pages 51–53. University of Helsinki, 2020\.
  * [3] Armin Biere, Marijn Heule, Hans van Maaren, and Toby Walsh, editors.  Handbook of Satisfiability - Second Edition, volume 336 of Frontiers in Artificial Intelligence and Applications.  IOS Press, 2021.
  * [4] Marko Kleine Büning, Carsten Sinz, and David Faragó.  QPR verify: A static analysis tool for embedded software based on bounded model checking.  In Maria Christakis, Nadia Polikarpova, Parasara Sridhar Duggirala, and Peter Schrammel, editors, Software Verification - 12th International Conference, VSTTE 2020, and 13th International Workshop, NSV 2020, Los Angeles, CA, USA, July 20-21, 2020, Revised Selected Papers, volume 12549 of Lecture Notes in Computer Science, pages 21–32. Springer, 2020.
  * [5] Markus Chimani, Carsten Gutwenger, Michael Jünger, Gunnar W Klau, Karsten Klein, and Petra Mutzel.  The open graph drawing framework (ogdf).  Handbook of graph drawing and visualization, 2011:543–569, 2013\.
  * [6] Niklas Eén and Niklas Sörensson.  An extensible sat-solver.  In Enrico Giunchiglia and Armando Tacchella, editors, Theory and Applications of Satisfiability Testing, 6th International Conference, SAT 2003\. Santa Margherita Ligure, Italy, May 5-8, 2003 Selected Revised Papers, volume 2919 of Lecture Notes in Computer Science, pages 502–518, 2003.
  * [7] Nick Feng and Fahiem Bacchus.  Clause size reduction with all-uip learning.  In Luca Pulina and Martina Seidl, editors, Theory and Applications of Satisfiability Testing - SAT 2020 - 23rd International Conference, Alghero, Italy, July 3-10, 2020, Proceedings, volume 12178 of Lecture Notes in Computer Science, pages 28–45. Springer, 2020.
  * [8] Marijn J. H. Heule.  The DRAT format and drat-trim checker.  CoRR, abs/1610.06229, 2016.
  * [9] Marijn J. H. Heule.  Schur number five.  In Proc. AAAI, pages 6598–6606, 2018.
  * [10] Marijn J. H. Heule, Oliver Kullmann, and Victor W. Marek.  Solving and verifying the boolean pythagorean triples problem via cube-and-conquer.  In Proc. SAT, pages 228–245, 2016.
  * [11] T. C. Hu and K. Moerder.  Multiterminal Flows in a Hypergraph.  In T.C. Hu and E.S. Kuh, editors, VLSI Circuit Layout: Theory and Design, chapter 3, pages 87–93. 1985.
  * [12] Daniela Kaufmann and Armin Biere.  Amulet 2.0 for verifying multiplier circuits.  In Jan Friso Groote and Kim Guldstrand Larsen, editors, Tools and Algorithms for the Construction and Analysis of Systems - 27th International Conference, TACAS 2021, Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2021, Luxembourg City, Luxembourg, March 27 - April 1, 2021, Proceedings, Part II, volume 12652 of Lecture Notes in Computer Science, pages 357–364. Springer, 2021\.
  * [13] Donald E Knuth.  The art of computer programming, Volume 4, Fascicle 6: Satisfiability.  Addison-Wesley Professional, 2015.
  * [14] Zack Newsham, William Lindsay, Vijay Ganesh, Jia Hui Liang, Sebastian Fischmeister, and Krzysztof Czarnecki.  Satgraf: Visualizing the evolution of SAT formula structure in solvers.  In Marijn Heule and Sean A. Weaver, editors, Theory and Applications of Satisfiability Testing - SAT 2015 - 18th International Conference, Austin, TX, USA, September 24-27, 2015, Proceedings, volume 9340 of Lecture Notes in Computer Science, pages 62–70. Springer, 2015.
  * [15] Chanseok Oh.  Between SAT and UNSAT: the fundamental difference in CDCL SAT.  In Marijn Heule and Sean A. Weaver, editors, Theory and Applications of Satisfiability Testing - SAT 2015 - 18th International Conference, Austin, TX, USA, September 24-27, 2015, Proceedings, volume 9340 of Lecture Notes in Computer Science, pages 307–323. Springer, 2015.
  * [16] Ezequiel Orbe, Carlos Areces, and Gabriel G. Infante López.  isat: Structure visualization for SAT problems.  In Nikolaj S. Bjørner and Andrei Voronkov, editors, Logic for Programming, Artificial Intelligence, and Reasoning - 18th International Conference, LPAR-18, Mérida, Venezuela, March 11-15, 2012. Proceedings, volume 7180 of Lecture Notes in Computer Science, pages 335–342. Springer, 2012.
  * [17] Usha Nandini Raghavan, Réka Albert, and Soundar Kumara.  Near Linear Time Algorithm to Detect Community Structures in Large-Scale Networks.  Physical Review E, 76(3):036106, 2007.
  * [18] Dominik Schreiber.  Lilotane: A lifted sat-based approach to hierarchical planning.  J. Artif. Intell. Res., 70:1117–1181, 2021.
  * [19] Daniel G. Schweikert and Brian W. Kernighan.  A Proper Model for the Partitioning of Electrical Circuits.  pages 57–62, 6 1972.
  * [20] Carsten Sinz.  Visualizing SAT instances and runs of the DPLL algorithm.  J. Autom. Reason., 39(2):219–243, 2007.
  * [21] Johannes Werner, Tomás Balyo, Markus Iser, and Michael Klein.  Fast approximate calculation of valid domains in a satisfiability-based product configurator.  In Michel Aldanondo, Andreas A. Falkner, Alexander Felfernig, and Martin Stettinger, editors, Proceedings of the 23rd International Configuration Workshop (CWS/ConfWS 2021), Vienna, Austria, 16-17 September, 2021, volume 2945 of CEUR Workshop Proceedings, pages 24–32. CEUR-WS.org, 2021.


Modifications of FastICA in Convolutive Blind Source Separation

YunPeng  Li

YunPeng Li was with the Department of Automation, Tsinghua University, Beijing, China e-mail<EMAIL_ADDRESS>
Shell et al.: Bare Demo of IEEEtran.cls for IEEE Journals

Convolutive blind source separation (BSS) is intended to recover the unknown components from their convolutive mixtures. Contrary to the  contrast functions used in instantaneous cases, the  spatial-temporal prewhitening stage and the para-unitary filters constraint are difficult to implement in a convolutive context. In this paper, we propose several modifications of FastICA to alleviate these difficulties. Our method performs the simple prewhitening step on convolutive mixtures prior to the separation and optimizes the contrast function under the diagonalization constraint implemented by single value decomposition (SVD). Numerical simulations are implemented to verify the performance of the proposed method.

Convolutive blind source separation (BSS), FastICA, diagonalization constraint, singular value decomposition (SVD).

§ INTRODUCTION

Convolutive blind source separation (BSS) has received much attention in recent years. It has been successfully applied in many signal processing problems, such as sonar array processing, seismic exploration, and the "cocktail party problem". In the convolutive BSS, $n$ dimensional observation $\mathbf{x}(k)=\left(x_{1}(k),\cdots,x_{n}(k) \right)^{T}$ is the convolutive mixtures of the $m$ dimensional independent source $\mathbf{s}(k)=\left(s_{1}(k),\cdots,s_{m}(k) \right)^{T}$. The unknown mixing process $\mathbf{A}$ can be described by a Multi-Input Multi-Output (MIMO) linear time invariant (LTI) system $\left( \mathbf{A}(k)_{k\in Z}\right)$,

\begin{equation}
\label{eq:mixing}
\mathbf{x}(k)=\mathbf{A}(k)*\mathbf{s}(k)=\sum_{l\in Z}\mathbf{A}(l)\mathbf{s}(k-l)
\end{equation}

Given $N$ samples of observation $\mathbf{x}(k)$, recovering the sources' estimation $\mathbf{y}(k)=\left(y_{1}(k),\cdots,y_{m}(k) \right)^{T}$ from these mixtures is equivalent to find a filter banks $\left( \mathbf{B}(k)_{k\in Z}\right)$ to
inverse the mixing system $\left( \mathbf{A}(k)_{k\in Z}\right)$:

\begin{equation}
\label{eq:unmixng}
\mathbf{y}(k)=\mathbf{B}(k)*\mathbf{x}(k)=\sum_{l\in Z}\mathbf{B}(l)\mathbf{x}(k-l)
\end{equation}

each $y_{i}(k)$ is a scaling and filtering version of the unique source's component $s_{i^{'}}(k)$.
Although the ordering and scaling ambiguity in the instantaneous BSS can be efficiently handled by the prewhitening on $\mathbf{x}(k)$ and the orthogonal constraint of unmixing matrix $\mathbf{B}$, the filtering ambiguity caused by the time delay is difficult to deal with in convolutive context and these strategies (used in instantaneous BSS) become invalid.
These ambiguities lead to much more extreme points in convolutive BSS than in instantaneous BSS, making the convolutive BSS a tough problem to deal with.
In this paper, we shall assume that:

* The source signals $\mathbf{s}(k)$ are real-valued, zero-mean, and mutually statistically independent, at most one of them is Gaussian.

* The filter banks $\left( \mathbf{A}(k)_{k\in Z}\right)$, $\left( \mathbf{B}(k)_{k\in Z}\right)$ are stable, causal, and finite impulse response (FIR).

Many Methods[4][5] have been proposed to solve the convolutive BSS. Most of them can be classified into two groups: the frequency domain approaches and the time domain approaches. In frequency domain[6][7], convolutive BSS can be considered as instantaneous BSS for each frequency bin, where each bin has own scaling and ordering indeterminacy as mentioned before. Complex value after discrete Fourier transform (DFT) and circularity problem happen in frequency domian.
Time domain approaches include density matching methods and contrast function methods. Density matching methods apply the well-known InfoMax[8] proposed in ICA to the convolutive case[9]. The performance of the density matching approaches is highly dependent on
the prior knowledge on the unknown density distributions of $s_{i}$. It's important to determine whether the source $s_{i}(k)$ is super-Gaussian or
sub-Gaussian beforehand. Density matching methods linearly transform the observed mixtures $\mathbf{x}(k)$ with the demixing system  $\left( \mathbf{B}(k)_{k\in Z}\right)$, forcing the
$p_{y}(y_{i})$ close to a selected density $p_{s}(s_{i^{'}})$. Most of these methods are based on the gradient optimization, requiring appropriate choice of learning rate and step direction. Natural gradient[10] and relative gradient[11] have been proposed to alleviate the drawbacks of stability and convergence above. Contrast function methods make the advantage of the statistical independence or non-Gaussianity of the components of source to recover the unknown $s_{i}$. In addition, High order statistics (HOS)[12] of source's estimation $y_{i}(k)$ like cumulants, cross-cumulants, and cross-moments can work as contrast function in separation. Prewhitening stage and coefficents constraints are required to guarantee the uniqueness of the extracted components during the contrast function methods. For instantaneous case, prewhitening stages and orthogonal constraint can be implemented by Gram–Schmidt orthogonalization or singular value decomposition (SVD). While in the convolutive mixtures, the corresponding spatial-temporal prewhitening stage and
the para-unitary filters constraint are more difficult to realize.
FastICA[3] is a contrast function method in instantaneous BSS by maximizing the approximation of negentropy[15]. It works well in instantaneous case for its fast convergence and robustness. Unfortunately, it cannot directly be adapted to the convolutive BSS due to the filtering ambiguity caused by the time delay. Several convolutive extensions of FastICA have been proposed[13][14] in recent years. The method proposed in [13] conducts a convolutive prewhitening to the transformed observation $\mathbf{x}(k)$ at first, then removes the extracted estimations $y_{i}$ at previous steps from the mixtures $\mathbf{x}$ in a deflation mode, leading to an accumulation of estimation errors which may become excessive after a certain number of source extractions. In [14], a spatial-temporal prewhitening stage and the para-unitary constraint is implemented at the cost of high computation burden. The existing FastICA extensions in convolutive BSS are difficult and inefficient.

In this paper, we propose a novel extension of FastICA, under a simpler framework. The proposed modifications combine a convolutive prewhitening stage for observations with the diagonalization constraint in both deflation or symmetric mode.
The rest of the paper is organized as follows. In Section <ref>, we state several assumptions of our method and rearrange the convolutive mixtures into instantaneous mixtures. An optimization problem based on the FastICA is described in Section <ref>.
We describe the prewhitening strategy along with the diagonalization constraint in Section <ref>. Experiments' results are presented to verify our modifications in Section <ref>. We conclude our method in Section <ref>.

§ PROBLEM STATEMENT

We describe the mixing system $\mathbf{A}$ by the following FIR filter equation,

\begin{equation}
\label{x}
x_{i}(k) =\sum_{j=1}^{m}\sum_{l=0}^{P-1}a_{ij}(l)s_{j}(k-l) \quad \quad i=1,\cdots,n
\end{equation}

where the $\mathbf{a}_{ij}$ are the mixing filters. Without loss of generality, we assume all the mixing filters have the same filter order $P$.
The estimation of the demixing system $\mathbf{B}$ shares the same form,

\begin{equation}
\label{y}
y_{i}(k) =\sum_{j=1}^{n}\sum_{l=0}^{Q-1}b_{ij}(l)x_{j}(k-l) \quad \quad i=1,\cdots,m
\end{equation}

where the $\mathbf{b}_{ij}$ are the demixing filters and we assume all the demixing filters have the same filter order $Q$. $y_{i}(k)$ is a unique scaled, permuted, and filtered vesion of$s_{i^{'}}(k)$.
To reconstruct the contributions in each observation $x_{i}(k)$, another assumption is considered:

* Each signal source $s_{i}(k)$ is produced by an innovation process $u_{i}(k)$ via a stable FIR filters $F_{i}(k)$, where $u_{i}(k)$ is zero-mean, mutually independent and non-gaussian random process.

Then, signal source $s_{i}(k)$ can be expressed in the following form,

\begin{equation}
\label{innovation}
s_{i}(k)=\sum_{l=-R+1}^{R-1}F_{i}(l)u_{i}(k-l) \quad \quad i=1,\cdots,m
\end{equation}

The order of the non causal  FIR $F_{i}$ is $2R-1$, all the  filter banks   $\mathbf{F}(k)=diag(F_{1}(k),\cdots,F_{m}(k))$ are connected with
innovation process $\mathbf{u}(k)=\left(u_{1}(k),\cdots,u_{m}(k) \right)^{T}$ respectively. In order to individually extract $s_{i}$, we regard the $\mathbf{F}(k)$ as the coloring filters for the innovation process $\mathbf{u}(k)$.
To rearrange the convolutive BSS to instantaneous BSS, we created several variables as follow,

\begin{equation}
\label{s1}
\underline {\mathbf{s}}(k)=\left(\mathbf{s}_{1}^{T}(k),\mathbf{s}_{2}^{T}(k),\cdots,\mathbf{s}_{m}^{T}(k)\right)^{T}
\end{equation}

\begin{equation}
\label{s2}
\mathbf{s}_{i}(k) = \left(s_{i}(k),s_{i}(k-1),\cdots,s_{i}(k-P+1)\right)^{T}
\end{equation}

where $\underline{\mathbf{s}}(k)$ is  a $mP\times1$ column vector, and $\mathbf{s}_{i}(k)$ is a $P\times1$ column vecto. The mixing system $\mathbf{A}$ can be expressed in the matrix form,

\begin{equation}
\label{A1}
\mathbf{A}= \left(\begin{IEEEeqnarraybox*}[][c]{,c/c/c/c,}
\mathbf{a}_{11}^{T}&\mathbf{a}_{12}^{T}&\cdots &\mathbf{a}_{1m}^{T}\\
\mathbf{a}_{21}^{T}&\mathbf{a}_{22}^{T}&\cdots &\mathbf{a}_{2m}^{T}\\
\vdots&\vdots &\ddots &\vdots\\
\mathbf{a}_{n1}^{T}&\mathbf{a}_{n2}^{T}&\dots &\mathbf{a}_{nm}^{T}%
\end{IEEEeqnarraybox*}\right)
\end{equation}

\begin{equation}
\label{A2}
\mathbf{a}_{ij} = \left(a_{ij}(0),a_{ij}(1),\cdots,a_{ij}(P-1)\right)^{T}
\end{equation}

where $\mathbf{A}$ is a $n\times mP$ matrix, and $\mathbf{a}_{ij}$ is $P\times1$ column vector. The mixing system (<ref>) becomes,

\begin{equation}
\label{mix1}
\mathbf{x}(k)=\mathbf{A}\underline{\mathbf{s}}(k)
\end{equation}

The demixing process can be transformed in the same way.

\begin{equation}
\label{x1}
\underline {\mathbf{x}}(k)=\left(\mathbf{x}_{1}^{T}(k),\mathbf{x}_{2}^{T}(k),\cdots,\mathbf{x}_{n}^{T}(k)\right)^{T}
\end{equation}

\begin{equation}
\label{x2}
\mathbf{x}_{i}(k) = \left(x_{i}(k),x_{i}(k-1),\cdots,x_{i}(k-Q+1)\right)^{T}
\end{equation}

where $\underline{\mathbf{x}}(k)$ is  a $nQ\times1$ column vector, and $\mathbf{x}_{i}(k)$ is a $Q\times1$ column vector.

\begin{equation}
\label{B1}
\mathbf{B}= \left(\begin{IEEEeqnarraybox*}[][c]{,c/c/c/c,}
\mathbf{b}_{11}^{T}&\mathbf{b}_{12}^{T}&\cdots &\mathbf{b}_{1n}^{T}\\
\mathbf{b}_{21}^{T}&\mathbf{b}_{22}^{T}&\cdots &\mathbf{b}_{2n}^{T}\\
\vdots&\vdots &\ddots &\vdots\\
\mathbf{b}_{m1}^{T}&\mathbf{b}_{m2}^{T}&\dots &\mathbf{b}_{mn}^{T}%
\end{IEEEeqnarraybox*}\right)
\end{equation}

\begin{equation}
\label{B2}
\mathbf{b}_{ij} = \left(b_{ij}(0),b_{ij}(1),\cdots,b_{ij}(Q-1)\right)^{T}
\end{equation}

where $\mathbf{B}$ is a $m\times nQ$ matrix, and $\mathbf{b}_{ij}$ is $Q\times1$ column vector.

\begin{equation}
\label{demix1}
\mathbf{y}(k)=\mathbf{B}\underline{\mathbf{x}}(k)
\end{equation}

Owing to the concise expressions of the mixing process (<ref>) and demixng process (<ref>), we have converted
the $m$ signals $n$ observations covolutive BSS into a $m$ signals $nQ$ observations instantaneous BSS problem. For the sake of uniqueness in extraction, a prewhitening stage is required in contrast function methods:

\begin{equation}
\label{prew}
\mathbf{v}(k) =\mathbf{H}(\mathbf{x}(k))
\end{equation}

we represent the prewhitening stage in function $\mathbf{H}()$, and the prewhitening ouput $\mathbf{v}(k)$ should comply with some constraint. In instantaneous case, the particular step can be conducted via eigenvalue decomposition or PCA,

\begin{equation}
\label{prew2}
\mathbf{v}(k) =\mathbf{H}\mathbf{x}(k)
\end{equation}

where $\mathbf{v}(k)$ conforms to the below constraint.

\begin{equation}
\label{prew3}
\boldmath{E}\{\mathbf{v}(k)\mathbf{v}^{T}(k)\}=\mathbf{I}_{n\times n}
\end{equation}

Unfortunately, the above prewhitening strategy fails in the convolutive case, more details are provided in Section <ref>.

\begin{equation}
\label{sep}
\mathbf{y}(k)=\mathbf{B}\underline{\mathbf{x}}(k)=\mathbf{W}\mathbf{H}\underline{\mathbf{x}}(k)=\mathbf{W}\mathbf{v}(k)
\end{equation}

After prewhitening stage, our method (<ref>) adjusts coefficients in separation matrix $\mathbf{W}$ to recover the $y_{i}(k)$ as an estimation of a delayed scaled innovation process $u_{i^{'}}(k-l)$ .

Flow chart of the convolutive BSS.

The complete routine of the convolutive BSS can be concluded in Fig.<ref>. The source's component $\mathbf{s}(k)$ is the output of
the innovation process $\mathbf{u}(k)$ filtered by the coloring filters $\mathbf{F}$. Unknown mixing system $\mathbf{A}$ mixes the $\mathbf{s}(k)$ both in time and space. Given the observation
$\mathbf{x}(k)$, the proposed method conducts a demixing procedure $\mathbf{B}$ to recover $\mathbf{y}(k)$. The prewhitening $\mathbf{H}$ stages and seperating $\mathbf{W}$ stages are the keys in the demixing system.

§ FASTICA EXTENSION

In the instantaneous case, FastICA looks for a sequence of orthogonal projections to maximize the negentropy  $J(y_{i})$ [15], which amounts to seek components as independent as possible.

\begin{equation}
\label{neg}
\begin{aligned}
&\approx c\left[\boldmath{E}\{G(y_{i}) \}-\boldmath{E}\{G(z_{i})\}\right]^{2}
\end{aligned}
\end{equation}

where  $J(y_{i})$ is the negentropy, $\boldmath{H}()$ is the random variable's entropy, $c$ is an irrelevant constant, $G()$ is any non-quadratic function, $z_{i}$ is a Gaussian random variable with the same variance as $y_{i}$. Mututal information $\boldmath{I}(\mathbf{y})$ between the components of the random variable vector $\mathbf{y}(k)$ is a nature measure of dependence. It is always non-negative and becomes zero only when the components are statistically independent.

\begin{equation}
\label{info}
\begin{aligned}
\boldmath{I}(\mathbf{y})&=\sum_{i=1}^{m}\boldmath{H}(y_{i})-\boldmath{H}(\mathbf{y})\\
\end{aligned}
\end{equation}

where $C$ is an irrelevant constant, the minimization of the mutual information $\boldmath{I}(\mathbf{y})$
for independence (under the constraint of decorrelation) is equivalent to the maximization of the sum of the negentropies of the components $\sum_{i=1}^{m}J(y_{i})$, and the original FastICA in instantaneous can be modeled as the optimization problem[3] below,

\begin{equation}
\label{insopt}
\begin{aligned}
\max \quad & \sum_{i=1}^{m}J(y_{i}) & \\
\mathrm{s.t.} \quad & \boldmath{E}\{y_{i}(k)y_{j}(k)\} =\delta_{ij}, &\quad i,j=1,2,\ldots,m \,.
\end{aligned}
\end{equation}

where $\delta_{ij}$ is the item in identity matrix $\mathbf{I}_{m\times m}$. The objective function in (<ref>) aims at the maximization of independence. In both the deflation and symmetric mode, the optimization problem in (<ref>) is divided as single maximization of $J(y_{i})$, and there are 2$m$  extreme points to this problem due to the ordering and scaling ambiguities in instantaneous case. The constraint in (<ref>) is designed to avoiding extracting the same solution more than once.

While in the convolutive case, the filtering ambiguity introduces much more extreme points as the increasing of the demixing filters $\mathbf{B}$'s filter order $Q$, resulting in the failure of the same strategies in instantaneous BSS. It's straightforward to change the equation (<ref>) in the convolutive context,

\begin{equation}
\label{conopt}
\begin{aligned}
\max \quad & \sum_{i=1}^{m}J(y_{i}) & \\
\mathrm{s.t.} \quad & \boldmath{E}\{y_{i}(k)y_{i}(k)\} =1, &\quad i=1,2,\ldots,m \,,\\
& \boldmath{E}\{y_{i}(k)y_{j}(k-l)\} =0, &\quad i\neq j,-\infty<l<+\infty \,.
\end{aligned}
\end{equation}

The first constraint in (<ref>) is for the scaling ambiguity, and the second constraint is designed to tack the ordering and filtering ambiguities. In order to simplify the second constraint, we choose $-L\leq l\leq L$, where $L$ is a positive large enough integer. Prewhitening stage and diagonalization constraints in our modifications are used to satisfy these constraints.

§ PREWHITENING AND DIAGONALIZATION CONSTRAINTS

Many contrast function methods for convolutive mixtures have a prewhiteing stage (<ref>) in space and time[4][14].

\begin{equation}
\boldmath{E}\{\mathbf{v}(k)\mathbf{v}^{T}(k-l)\}=\delta_{l}\mathbf{I}_{n\times n} \quad \forall l \in Z
\end{equation}

The whitening filter $\mathbf{H}$ is not unique and hard to determine, the contrast functions are required to be  optimized under the constraint of complicate para-unitary filters $\mathbf{W}$. Although several methods have been proposed to alleviated these difficulties[16][17], they are both difficult and demands mass computing. We conduct a same prewhitening stage in instantaneous case[6] on $\underline {\mathbf{x}}(k)$ to produce the $nQ\times1$ column vector $\mathbf{v}(k)$, which is regarded as the convolutive prewhitening.

\begin{equation}
\label{precov}
\boldmath{E}\{\mathbf{v}(k)\mathbf{v}^{T}(k)\}=\mathbf{I}_{nQ\times nQ}
\end{equation}

After the above prewhitening stage, the main effort is to adjust the seperating filters $\mathbf{W}$ to optimize equation (<ref>),

\begin{equation}
\label{conopt2}
\begin{aligned}
\max \quad & \sum_{i=1}^{m}J(\mathbf{w}^{T}_{i}\mathbf{v}(k)) & \\
\mathrm{s.t.} \quad & \mathbf{w}^{T}_{i}\mathbf{w}_{i} =1, &\quad i=1,2,\ldots,m \,,\\
& \mathbf{w}^{T}_{i}\boldmath{E}\{\mathbf{v}^{T}(k)\mathbf{v}^{T}(k-l)\}\mathbf{w}_{j} =0, &\quad i\neq j,-L\leq l\leq L \,.
\end{aligned}
\end{equation}

where $\mathbf{w}^{T}_{i}$ is the $i$th row of $m\times nQ$ matrix $\mathbf{W}$. In pursuit of particular $\mathbf{w}_{i}$, the following iteration routine is carried out until convergence.

\begin{equation}
\label{fix}
\begin{aligned}
\mathbf{w}_{i}&=\boldmath{E}\{\mathbf{v}(k)g(\mathbf{w}^{T}_{i}\mathbf{v}(k))\}-\boldmath{E}\{g^{'}(\mathbf{w}^{T}_{i}\mathbf{v}(k))\}\mathbf{w}_{i}\\
\mathbf{w}_{i}&=\mathbf{w}_{i}/\|\mathbf{w}_{i}\|_{2}
\end{aligned}
\end{equation}

where $g()$ is the derivative of particular non-quadratic function $G()$, coefficients constraints are required during the process of (<ref>) in deflation and symmetric mode.
The validity of constraint in (<ref>)(<ref>) can also supported from [6][18] the fact: If the source signals have unique temporal structures or non-stationary, simultaneous diagonalization of ouput correlation matrices over multiple time lags can separate the independent sources from convolutive mixtures.

\begin{equation}
\label{Ry}
\begin{aligned}
\mathbf{R}_{y}(\tau)&=\boldmath{E}\{\mathbf{y}(k)\mathbf{y}^{T}(k-\tau)\}\\
\end{aligned}
\end{equation}

The output correlation matrix at time lag $\tau$ is represented as $\mathbf{R}_{y}(\tau)$, it is required to be a diagonal matrix, particularly,$\mathbf{R}_{y}(0)$ is the identity matrix $\mathbf{I}_{m\times m}$. For more intuitive explanation, we define several column vectors.

\begin{equation}
\label{y1}
\underline {\mathbf{y}}(k)=\left(\mathbf{y}_{1}^{T}(k),\mathbf{y}_{2}^{T}(k),\cdots,\mathbf{y}_{m}^{T}(k)\right)^{T}
\end{equation}

\begin{equation}
\label{y2}
\mathbf{y}_{i}(k) = \left(y_{i}(k+L),y_{i}(k+L-1),\cdots,y_{i}(k-L)\right)^{T}
\end{equation}

where $\underline{\mathbf{y}}(k)$ is a $m(2L+1)\times1$ column vector, and $\mathbf{y}_{i}(k)$ is a $(2L+1)\times1$ column vector, the equivalent expression in (<ref>) can be described below,

\begin{equation}
\label{Ry1}
\begin{aligned}
\mathbf{R}_{\underline{\mathbf{y}}}&=\boldmath{E}\{\underline{\mathbf{y}}(k)\underline{\mathbf{y}}^{T}(k)\}\\
\mathbf{R}_{\mathbf{y}_{1}\mathbf{y}_{1}}&\mathbf{R}_{\mathbf{y}_{1}\mathbf{y}_{2}}&\cdots &\mathbf{R}_{\mathbf{y}_{1}\mathbf{y}_{m}}\\
\mathbf{R}_{\mathbf{y}_{2}\mathbf{y}_{1}}&\mathbf{R}_{\mathbf{y}_{2}\mathbf{y}_{2}}&\cdots &\mathbf{R}_{\mathbf{y}_{2}\mathbf{y}_{m}}\\
\vdots&\vdots &\ddots &\vdots\\
\mathbf{R}_{\mathbf{y}_{m}\mathbf{y}_{1}}&\mathbf{R}_{\mathbf{y}_{m}\mathbf{y}_{2}}&\dots &\mathbf{R}_{\mathbf{y}_{m}\mathbf{y}_{m}}%
\end{IEEEeqnarraybox*}\right)
\end{aligned}
\end{equation}

\begin{equation}
\label{Ry2}
\mathbf{R}_{\mathbf{y}_{i}\mathbf{y}_{j}}=\boldmath{E}\{\mathbf{y}_{i}(k)\mathbf{y}^{T}_{j}(k)\}
\end{equation}

where $m(2L+1)\times m(2L+1)$ matrix $\mathbf{R}_{\underline{\mathbf{y}}}$ is the combination of correlation matrix concerning different time lags. According to the constraint in (<ref>) and the nonstationarity property of source signals, the $(2L+1)\times(2L+1)$ matrix $\mathbf{R}_{\mathbf{y}_{i}\mathbf{y}_{j}}$ becomes nonzero only in the diagonal position of $\mathbf{R}_{\underline{\mathbf{y}}}$. We consider the coefficents constraints above as diagonalization constraints.
Compared with the existing para-unitary constraint, the diagonalization constraints are relaxed, it can be efficiently imposed via singular value decomposition in both deflation and symmetric mode.

§.§ deflation mode

When one of the sources' estimation $y_{i}(k)$ has been extracted, it's necessary to substract its contribution from the observations to obtain the mixtures of $m-1$ sources, then we repeat this procedure to extract the remained sources one by one, until all the sources have been extracted.
Considering the constraints in (<ref>), we construct the block matrix $\mathbf{O}$ during the extraction of the $y_{i}(k)$ to guarantee uniqueness.

\begin{equation}
\label{block}
\mathbf{O} =\left(\mathbf{R}_{v}(-L)\mathbf{w}_{1},\mathbf{R}_{v}(-L+1)\mathbf{w}_{1},\cdots,\mathbf{R}_{v}(L)\mathbf{w}_{i-1}\right)
\end{equation}

\begin{equation}
\label{otrh}
\mathbf{w}^{T}_{i}\mathbf{O}=\mathbf{0}
\end{equation}

The $\mathbf{w}_{i}$ in the $i$th extraction is required to be orthogonal to the column space of the  block matrix $\mathbf{O}$, we represent the column space of $\mathbf{O}$ as $\boldmath{span}\{\mathbf{O}\}$.
If the block matrix $\mathbf{O}$ is overdetermined,the column space $\boldmath{span}\{\mathbf{O}\}$ is not full column rank, it's easy adjust the $\mathbf{w}_{i}$ via least square solution.

\begin{equation}
\label{lss}
\mathbf{w}_{i} =\mathbf{w}_{i} -\left(\mathbf{O}^{T}\mathbf{O}\right)^{-1}\mathbf{O}^{T}\mathbf{w}_{i}
\end{equation}

Unfortunately, $\mathbf{O}$ is often designed to be underdetermined due to choice of large $L$, so the strategy in (<ref>) cannot work any longer. The proposed method draws lessons from the principal component analysis based on singular value decomposition, only taking the directions with most variation in $\mathbf{O}$ into consideration.

\begin{equation}
\label{svd}
\mathbf{O} =\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{T}
\end{equation}

Here $\mathbf{U}$ and $\mathbf{V}$ are orthogonal matrices, with the columns of $\mathbf{U}$ spanning the column space of $\mathbf{O}$, and the columns of $\mathbf{V}$ spanning the row space. $\mathbf{\Sigma}$ is a diagonal matrix, with diagonal entries in decreasing order.

\begin{equation}
\sigma_{11}\ge\sigma_{22}\ge\cdots\ge0
\end{equation}

The algorithm picks the first $r$ columns vectors in $\mathbf{U}$ to obtain the most variance in original $\boldmath{span}\{\mathbf{O}\}$ based on the effective rank of $\mathbf{O}$, $r$ is the minimum integer when $\mu(r)$ is greater than particular threshold $\alpha$ (such as $0.99995$).

\begin{equation}
\label{effrank}
\mu(r) =\frac{\sqrt{\sigma_{11}^{2}+\cdots+\sigma_{rr}^{2}}}{\|\Sigma\|_{F}}
\end{equation}

The most variant $r$ column vectors in  $\mathbf{U}$ is represented as $\mathbf{U}_{(r)}$, and each column vector in $\mathbf{U}_{(r)}$ is orthogonal with each other, $\mathbf{U}_{(r)}$ is always high matrix,the (<ref>) can be described in simple form.

\begin{equation}
\label{def1}
\mathbf{w}_{i} =\mathbf{w}_{i} -\mathbf{U}_{(r)}\mathbf{U}_{(r)}^{T}\mathbf{w}_{i}
\end{equation}

The  FastICA convolutive algorithm in deflation mode is summarized in Alg.$\ref{alg:def}$, where $tol$ is the threshold in iteration (such as $10^{-7}$).

convolutive FastICA: deflation mode

$\mathbf{W}=\mathbf{0}_{m\times nQ}$
$\mathbf{w}_{i}=\mathbf{0}_{nQ\times 1}$
$\mathbf{O} =\left(\mathbf{R}_{v}(-L)\mathbf{w}_{1},\mathbf{R}_{v}(-L+1)\mathbf{w}_{1},\cdots,\mathbf{R}_{v}(L)\mathbf{w}_{i-1}\right)$
calculate $\mathbf{U}_{(r)}$ from (<ref>)(<ref>) with threshold $\alpha$
$\mathbf{w}_{i} =\mathbf{w}_{i} -\mathbf{U}_{(r)}\mathbf{U}_{(r)}^{T}\mathbf{w}_{i}$
$||\mathbf{w}^{T}_{i}\mathbf{w}^{'}_{i}|-1.0|\leq tol$

§.§ symmetric mode

Every extracted source has different priority in deflation mode, leading an accumulation of estimation errors. Symmetric mode is proposed to extract $\mathbf{w}_{i}$ equally. Supposing only the $\mathbf{w}_{i}$ is in process, other extracted sources are fixed. It's nature to construct the block matrix from (<ref>)(<ref>).

\begin{equation}
\label{block2}
\mathbf{O} =\left(\mathbf{R}_{v}(-L)\mathbf{w}_{1},\cdots,\mathbf{R}_{v}(L)\mathbf{w}_{i-1},\mathbf{R}_{v}(-L)\mathbf{w}_{i+1},\cdots\right)
\end{equation}

The FastICA convolutive algorithm in symmetric mode is similar to the deflation mode, and it can be summarized in Alg.<ref>.
convolutive FastICA: symmetric mode

$\mathbf{W}=\mathbf{0}_{m\times nQ}$
$\mathbf{O} =\left(\mathbf{R}_{v}(-L)\mathbf{w}_{1},\cdots,\mathbf{R}_{v}(L)\mathbf{w}_{i-1},\mathbf{R}_{v}(-L)\mathbf{w}_{i+1},\cdots\right)$
calculate $\mathbf{U}_{(r)}$ from (<ref>)(<ref>) with threshold $\alpha$
$\mathbf{w}_{i} =\mathbf{w}_{i} -\mathbf{U}_{(r)}\mathbf{U}_{(r)}^{T}\mathbf{w}_{i}$
$\||\mathbf{W}^{'}\mathbf{W}^{T}|-\boldmath{I}_{m\times m}\|_{2}\leq tol$

§.§ sources' reconstruction

After the extraction of $\mathbf{y}(k)$ in deflation or symmetric mode, we get the  estimation of the ordered, scaled, and delayed version of innovation process $\mathbf{u}(k)$ (<ref>), we then conduct a reconstruction process to recover the signals $\mathbf{s}(k)$'s contributions in the observation mixtures $\mathbf{x}(k)$. In this part,for the convenience of discussing, the ordering ambiguity is ignored.
For the purpose of achieve innovation $s_{i}(k)$'s contributions in observation $x_{j}(k)$, a $(2L+N)\times(2L+1)$ high matrix $\mathbf{T}$ is built to represent the column space via zero padding, including all the possible forward and backward time shift of $y_{i}$.

\begin{equation}
\label{uspace}
\mathbf{T} =\left(\begin{IEEEeqnarraybox*}[][c]{,c/c/c/c/c,}
y_{i}(0)&0&\cdots &0&0\\
y_{i}(1)&y_{i}(0)&\cdots &0&0\\
y_{i}(2)&y_{i}(1)&\cdots &0&0\\
\vdots&\vdots &\ddots &\vdots&\vdots\\
0&0&\cdots &y_{i}(N-2)&y_{i}(N-3)\\
0&0&\cdots &y_{i}(N-1)&y_{i}(N-2)\\
0&0&\cdots &0&y_{i}(N-1)%
\end{IEEEeqnarraybox*}\right)
\end{equation}

$N\times 1$ column vector $x_{j}$ is padding with $L$ zeros both forward and backward.

\begin{equation}
\label{xj2}
\end{equation}

The reconstruction is based on the regression opinion: regress $x_{j}^{(2L+N)}$ on the column space of $\mathbf{T}$, finding the closest $\hat{s}_{ij}$ in the $\boldmath{span}\{\mathbf{T}\}$ in least square sense. We describe the $s_{i}$ contribution to observation $x_{j}$ as $\hat{s}_{ij}$.

\begin{equation}
\label{regress}
\hat{s}_{ij} =(\mathbf{T}^{T}\mathbf{T})^{-1}\mathbf{T}^{T}x^{(2L+N)}_{j}
\end{equation}

§ EXPERIMENT RESULTS

We conducted two experiments in this section to explore the performance of the proposed algorithm.

[Source signals]

Source signals and observations in the 2$\times$2 case(simulation).

The first experiment mixed the given innovation process $u_{1}$ and $u_{2}$ in Fig.<ref>, both satisfying the referred assumptions. After the convolutive prewhitening stage, the algorithm in symmetric mode recovered the $y_{1}$ and $y_{2}$, and calculated their contributions on each observation.

[Innovation process]

[Sources' contirbutions]

Innovation process and sources' contributions in the 2$\times$2 case(symmetric mode).

In Fig.<ref> the recovered estimation $y_{1}$ and $y_{2}$ are similar to the original innovation process, while there are little deviations in the orignal zero amplitude regions. The contributions in Fig.<ref> happened to be the same as innovation process. In this simulation.the algorithm accomplished the mission of convolutive BSS.
A more difficult task was considered in real recorded source signals from the public data of Salk Institute[19].

* $s_{1}$ :speaker says the digits from one to ten in English.

* $s_{2}$ :loud music in the background.

[Source signals]

Source signals and observations in the 2$\times$2 case(record).

We conduct deflation mode in Alg.<ref> to produce the estimations and contributions. The results in Fig.<ref> showed the similarity between the calculated contributions and original source signals.

[Innovation process]

[Sources' contirbutions]

Innovation process and sources' contributions in the 2$\times$2 case(deflation mode).

§ CONCLUSION

In this paper, we have derived a novel extension of the FastICA for convolutive mixtures that enforces the diagnoalization constraints on the seperating filters for uniqueness. Our algorithm also has simple convolutive prewhitening stages and contributions' reconstruction procedure. Experiments are given to illustrate the performance of the proposed algorithm.
Our algorithm enjoys the robustness and fast convergence due to its fix-point iterations, no particular parameter tuning is required during the optimization. Compared with spatial-temporal prewhitening stages and para-unitary filter constraints in other contrast function methods[4], the corresponding procedures in our algorithm are much more straightforward and simpler.

[2]
P. Comon, “Independent component analysis, a new concept?” Signal
Process., vol. 36, no. 3, p. 287–314, Apr. 1994.

[3]
A. Hyvarinen, “Fast and robust fixed-point algorithms for independent
component analysis,” IEEE Transactions on Neural Networks, vol. 10,
no. 3, pp. 626–634, 1999.

[4]
M. Castella, A. Chevreuil, and J.-C. Pesquet, “Chapter 8 - convolutive
mixtures,” in Handbook of Blind Source Separation, P. Comon and
C. Jutten, Eds.1em plus 0.5em minus 0.4emOxford: Academic
Press, 2010, pp. 281 – 324.

[5]
M. S. Pedersen, “A survey of convolutive blind source separation methods,”
Springer Handbook on Speech Communication, 2007.

[6]
D. Yellin and E. Weinstein, “Criteria for multichannel signal
separation,” IEEE Transactions on Signal Processing, vol. 42, no. 8,
pp. 2158–2168, 1994.

[7]
P. Smaragdis, “Blind separation of convolved mixtures in the frequency
domain,” Neurocomputing, vol. 22, no. 1, pp. 21 – 34, 1998.

[8]
A. J. Bell and T. J. Sejnowski, “An information-maximization approach to
blind separation and blind deconvolution,” Neural Computation,
vol. 7, no. 6, pp. 1129–1159, 1995.

[9]
K. Torkkola, “Blind separation of delayed sources based on information
maximization,” 1996, pp. 3509–3512.

[10]
S.-I. Amari, “Natural gradient works efficiently in learning,” Neural
Comput., vol. 10, no. 2, p. 251–276, Feb. 1998.

[11]
J. . Cardoso, “Blind signal separation: statistical principles,”
Proceedings of the IEEE, vol. 86, no. 10, pp. 2009–2025, 1998.

[12]
P. Comon, “Contrasts for multichannel blind deconvolution,” IEEE
Signal Processing Letters, vol. 3, no. 7, pp. 209–211, 1996.

[13]
J. Thomas, Y. Deville, and S. Hosseini, “Time-domain fast fixed-point
algorithms for convolutive ica,” IEEE Signal Processing Letters,
vol. 13, no. 4, pp. 228–231, 2006.

[14]
S. C. Douglas, M. Gupta, H. Sawada, and S. Makino, “Spatio–temporal
fastica algorithms for the blind separation of convolutive mixtures,”
IEEE Transactions on Audio, Speech, and Language Processing, vol. 15,
no. 5, pp. 1511–1520, 2007.

[15]
A. Hyvärinen, “New approximations of differential entropy for independent
component analysis and projection pursuit,” in Proceedings of the 1997
Conference on Advances in Neural Information Processing Systems 10, ser.
NIPS ’97.1em plus 0.5em minus 0.4emCambridge, MA, USA: MIT
Press, 1998, p. 273–279.

[16]
P. Comon and L. Rota, “Blind separation of independent sources from
convolutive mixtures,” IEICE Transactions on Fundamentals of
Electronics, Communications and Computer Sciences, vol. E86A, pp. 542–549,
03 2003.

[17]
J. G. McWhirter, P. D. Baxter, T. Cooper, S. Redif, and J. Foster,
“An evd algorithm for para-hermitian polynomial matrices,” IEEE
Transactions on Signal Processing, vol. 55, no. 5, pp. 2158–2169, 2007.

[18]
M. Kawamoto, K. Matsuoka, and N. Ohnishi, “A method of blind separation for
convolved non-stationary signals,” Neurocomputing, vol. 22, no. 1,
pp. 157 – 171, 1998.

[19]
T.-W. Lee, “Blind source separation of recorded speech and music signals,”