Citation
A real-time implementation of fastica for blind source separation of audio signals

Material Information

Title:
A real-time implementation of fastica for blind source separation of audio signals
Creator:
Jackson, Anthony Roger
Publication Date:
Language:
English
Physical Description:
xii, 172 leaves : ; 28 cm

Subjects

Subjects / Keywords:
Principal components analysis ( lcsh )
Multivariate analysis ( lcsh )
Computer sound processing ( lcsh )
Computer sound processing ( fast )
Multivariate analysis ( fast )
Principal components analysis ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 171-172).
General Note:
Department of Electrical Engineering
Statement of Responsibility:
by Anthony Roger Jackson.

Record Information

Source Institution:
University of Colorado Denver
Holding Location:
Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
317586900 ( OCLC )
ocn317586900
Classification:
LD1193.E54m 2008m J32 ( lcc )

Full Text
A REAL-TIME IMPLEMENTATION OF FASTICA FOR BLIND SOURCE
SEPARATION OF AUDIO SIGNALS
by
Anthony Roger Jackson
B.S., Brigham Young University, 2004
A thesis submitted to the
University of Colorado Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Electrical Engineering
2008


This thesis for the Master of Science
degree by
Anthony Roger Jackson
has been approved
by
Miloje S. Radenkovic
u/iyfa>


Jackson, Anthony Roger (M.S., Electrical Engineering)
A Real-time Implementation of FastICA for Blind Source Separation of Audio
Signals
Thesis directed by Associate Professor Miloje S. Radenkovic
ABSTRACT
FastICA is a relatively new algorithm for performing independent compo-
nent analysis. Background information into the theory and operation of ICA
and FastICA is presented, particularly as applied to blind source separation of
audio signals. A real-time implementation of FastICA suitable for the separation
and playback of audio signals is also presented and discussed at length. Audio
hardware is interfaced in real-time via libjack, the Jack Audio Connection Kit.
Qualitative and quantitative analyses follow the implementation. Full source
code is included to the extent that licensing allows.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
Miloje S. Radenkovic


DEDICATION
I dedicate this work foremost to my wife Rachel. Without her loving support
this would not be a reality. I further dedicate this thesis to my daughter Sabrina,
and my new son Oliver for their special contributions, i.e., hugs and kisses.


ACKNOWLED GMENT
I would like to thank Dr. Radenkovic for his support over the last three years
as an exceptional instructor and as a patient mentor for this thesis.
I also owe a debt of gratitude to EchoStar Corporation for financial support
and for allowing me to follow a flexible schedule so that I could attend classes
during normal working hours.
Finally, I want to thank all of those who have contributed to the open
source sound server JACK, the ubiquitous libsndfile, and to the excellent JACK
front-end qjackctl, without which this thesis would not have been possible.


CONTENTS
Figures ............................................................... xi
Tables................................................................ xii
Chapter
1. Introduction........................................................ 1
1.1 Contributions....................................................... 2
1.2 Literature Review................................................... 3
2. Background Information.............................................. 5
2.1 Model............................................................... 5
2.2 Assumptions ........................................................ 6
2.2.1 Independence..................................................... 6
2.2.2 Linearity ....................................................... 7
2.2.3 Ambiguities...................................................... 7
2.2.4 Observations..................................................... 8
2.2.5 Nongaussianity................................................... 8
2.2.5.1 Kurtosis........................................................ 9
2.2.5.2 Negentropy..................................................... 10
2.2.5.3 Approximating Negentropy..................................... 11
2.3 Preprocessing.................................................. 12
2.3.1 Centering....................................................... 12
2.3.2 Principal Component Analysis ................................... 12
vi


2.3.3 Whitening...................................................... 14
2.4 FastICA Algorithm.............................................. 15
2.4.1 Single Component FastICA....................................... 15
2.4.2 Symmetric FastICA ............................................. 19
2.5 Selection of G(-) 20
3. Implementation.................................................... 22
3.1 JACK........................................................... 22
3.2 Logical Application Programs................................... 23
3.2.1 Back-end Application Programs.................................. 27
3.2.2 Front-end Application Programs ................................ 28
3.3 Implementation Details fasticax.............................. 29
3.3.1 General Matrix Functions....................................... 31
3.3.2 FastICA Specific Functions .................................... 34
3.4 Implementation Details jackutils............................. 35
4. Analysis.......................................................... 37
4.1 Qualitative Analysis........................................... 37
4.1.1 Auditory Analysis ............................................. 37
4.1.2 Visual Analysis................................................ 38
4.2 Quantitative Analysis ......................................... 38
5. Conclusion........................................................ 45
5.1 Analysis of Implementation..................................... 45
5.2 Future Research................................................ 45
5.3 Main Conclusions............................................... 46
Appendix
vii


A. Source Code fastica_c.................................................... 47
A.l fastica_c/api.c.............................................................. 47
A.2 fastica_c/cov.c.............................................................. 48
A.3 fastica_c/cov.h ............................................................. 51
A.4 fastica_c/cp.c............................................................... 51
A.5 fastica_c/cp.h............................................................... 52
A.6 fastica_c/det.c.............................................................. 52
A.7 fastica_c/det.h.............................................................. 54
A.8 fastica_c/eig.c.............................................................. 54
A.9 fastica_c/eig.h.............................................................. 55
A.10 fastica_c/fastica-c.c....................................................... 56
A.11 fastica_c/fastica.h......................................................... 61
A.12 fastica.c/fastica-music, c ................................................. 62
A.13 fastica_c/ica.c............................................................. 68
A.14 fastica_c/ica.h............................................................. 75
A.15 fastica_c/ident.c........................................................... 75
A. 16 fastica.c/ident .h......................................................... 76
A. 17 fastica_c/inv.c............................................................ 76
A.18 fastica_c/invdr............................................................. 79
A. 19 fastica.c/jacobi.h......................................................... 79
A.20 fastica_c/libfastica.c ..................................................... 79
A. 21 fastica_c/Makefile......................................................... 81
A.22 fastica_c/matrix.c.......................................................... 82
A.23 fasticax/matrix.h........................................................... 88
viii


A.24 fastica_c/matrix_utest.c................................................. 89
A.25 fastica_c/mean.c ....................................................... 106
A.26 fastica_c/mean.h ....................................................... 109
A.27 fastica_c/mult.c........................................................ 109
A.28 fastica_c/mult.h........................................................ 110
A.29 fastica_c/op.c.......................................................... 110
A.30 fastica_c/op.h.......................................................... Ill
A.31 fastica_c/pca.c......................................................... 112
A.32 fastica_c/pca.h ........................................................ 116
A.33 fastica_c/real.h........................................................ 116
A.34 fastica_c/sqrt.c........................................................ 116
A.35 fastica_c/sqrt.h........................................................ 118
A.36 fastica_c/trans.c....................................................... 118
A.37 fastica_c/trans.h ...................................................... 119
A.38 fastica_c/var.c......................................................... 119
A.39 fastica_c/var.h......................................................... 120
A.40 fastica.c/whiten.c...................................................... 121
A. 41 fastica_c/whiten.h.................................................... 123
B. Source Code jackutils................................................ 124
B.l jackutils/conn.c......................................................... 124
B.2 jackutils/conn.h......................................................... 125
B.3 jackutils/jackconnect.c.................................................. 125
B.4 jackutils/jackconnect.py................................................. 131
B.5 jackutils/jackconnect.sh ................................................ 133
IX


B.6 jackutils/jackfilestream................................................ 133
B.7 jackutils/.jackfilestream .............................................. 136
B.8 jackutils/jackfilestream.sh............................................. 136
B.9 jackutils/jackica.c..................................................... 136
B.10 jackutils/jackica.py................................................... 143
B.ll jackutils/jackica.sh................................................... 145
B.12 jackutils/jackmix.c.................................................... 145
B.13 jackutils/jackmix.py .................................................. 149
B.14 jackutils/jackmix.sh .................................................. 152
B.15 jackutils/jackrecord.c................................................. 152
B.16 jackutils/jackrecord.py................................................ 156
B.17 jackutils/jackrecord.sh................................................ 157
B.18 jackutils/jackstream.c ................................................ 157
B.19 jackutils/Makefile..................................................... 162
B.20 jackutils/mempool.c ................................................... 163
B.21 jackutils/mempool.h.................................................... 163
B.22 jackutils/mix_and_stream.sh............................................ 164
B.23 jackutils/rbuffer.c.................................................... 164
B. 24 jackutils/rbuffer.h.................................................. 168
C. Source Code icacorr.m ............................................... 169
References................................................................ 171
x


FIGURES
Figure
3.1 An example of the connections possible in JACK....................... 24
3.2 The expected connections between jackfilestream, jackmix, jackica,
and jackrecord........................................................ 27
3.3 Screenshot of the JACK connections between jackutils applications,
as presented by qjackctl, that are needed for BSS of audio signals. . 28
3.4 Screenshot of PyGTK GUIs for front-end jackutils applications and
qjackctl.............................................................. 29
4.1 Comparison of the first two components of s with corresponding com-
ponents of s. 5 seconds of each waveform is displayed................. 39
4.2 Comparison of the second two components of s with corresponding
components of s. 5 seconds of each waveform is displayed........... 40
4.3 Correlation coefficients between components of s and s................ 43
4.4 Waveform of 5 seconds of mixed sources.............................. 44
xi


TABLES
Table
2.1 Cumulants up to and including the 4th order, where yix E{x} and
x is a random variable that follows an arbitrary distribution. ... 10
3.1 Average elapsed time allocating and deallocating matrices over 3 tri-
als of 5000 iterations each............................................. 32
4.1 Mapping between s and s................................................. 38
4.2 Correlation coefficients between components of s and s................ 41
4.3 Correlation coefficients between components of s and x................ 42
xii


1. Introduction
Imagine attending a typical cocktail-party of a few dozen attendees in a
crowded room. In such a situation, several conversations would be independently
unfolding. Although the party goers would likely focus on one conversation at
a time, in reality, each would hear all the conversations simultaneously.
Now imagine that you want to record the conversations and separate the
superimposed voices into distinct audio tracks, one for each speaker. Doing so is
a difficult task known as the cocktail-party problem. Solving the cocktail-party
problem requires some form of Blind Source Separation or BSS. Blind means
that we know little or nothing about the signals that we wish to separate, audio
in this case. A Source is one of the original unmixed signals. In this example,
each speakers voice is a source.
A popular method for performing BSS is independent component analysis
or ICA. ICA represents a whole class of algorithms that are able to separate
statistically independent signals into individual components (in our case audio
tracks) under certain conditions. These conditions will be discussed throughout
this thesis. BSS is not the only application of ICA. It has been successfully
applied to several classes of problems including image feature extraction, digital
communications, and financial analysis.
FastICA is a relatively new ICA algorithm that is suitable for but not limited
to BSS. Therefore, FastICA can solve the cocktail-party problem as explained
above given that certain conditions are met. FastICA, though not necessarily
1


designed for real-time applications, can, in fact, be used to separate sources in
real-time. Returning to our cocktail-party illustration, doing BSS in real-time
would mean that the separated audio tracks (one for each speaker) would be
available for playback after only a small negligible delay.
Contributions of this author, and a brief literature review on FastICA are
provided. This thesis then introduces the theoretical background required to
generally understand ICA and to specifically implement FastICA. A real-time
implementation of FastICA, as applied to the cocktail-party problem, is included
in C and Python source code form. Finally, the implementation is explained in
detail and analyzed using real audio signals.
1.1 Contributions
The inventor of FastICA provides an implementation in the MATLAB
languagefl]. However, the real-time requirements of this thesis make the C
language a better choice. It is possible to mix MATLAB and C, but MATLAB
makes no attempt to optimize itself for real-time operation. Therefore, it was
necessary to re-implement the algorithm in C. Since C has no native support
for matrix operations, and since FastICA relies heavily on them, the library
fastica.c was written from scratch to provide real-time safe matrix operations.
This library, listed in appendix A for reference, can be licensed for use, from
the author, for future work in ICA or in any other problem space that requires
real-time matrix manipulations in C. Therefore, the fastica_c library is a major
contribution of this thesis.
To the knowledge of the author, no implementation of FastICA exists that
performs the algorithm in real-time, without the assistance of an FPGA (field-
2


programmable gate array). Proving that a non-FPGA assisted real-time imple-
mentation is possible is another contribution of this work. In fact, this implemen-
tation is written in highly portable C and only requires a way to interface with
the mixed sources in a real-time manner. Therefore, being written in portable
C is also a contribution since the code can run on virtually any platform that
has a C compiler be it an embedded system, a workstation, or a supercomputer.
1.2 Literature Review
Since FastICA was invented by Aapo Hyvarinen, most of the works discussed
in this section are partially or fully his work. For an exhaustive survey of inde-
pendent component analysis, and not specifically FastICA, see [2]. In addition
to [2], [3] is a very good introduction and overview of independent component
analysis and specifically FastICA. It explains, in introductory terms, important
concepts like nongaussianity, kurtosis, negentropy (and approximations of it),
ICA preprocessing, the algorithm itself, and applications of ICA.
FastICA was first introduced in [4]. This version of FastICA used kurtosis
only, a non-robust measure, as a measure of nongaussianity. A later paper,
namely [5], introduced an approximation of negentropy, a more robust measure
of nongaussianity, in addition to kurtosis. [6] justifies the new approximation of
negentropy used in [5] and later papers.
A later work, [7], extends the algorithm to complex valued signals although
only real valued signals are supported by fastica_c.
There are several different contrast functions (or cost functions) that can be
utilized by FastICA. The choice of contrast function can affect the performance
of FastICA. [8] provides a detailed explanation of the theoretical and practical
3


considerations that are made when choosing a contrast function.
Though not related to FastICA specifically, these ICA resources are quoted
too often to exclude from this review: [9], and [10]. [9] is a relatively early work
that is the basis for several more recent works on ICA. [10] describes a popular
method for ICA which utilizes gradient methods as opposed to FastICAs fixed-
point iteration. It too serves as the basis for several recent works on ICA.
[11] is also not FastICA specific, but gives additional insight into the fixed-
point iteration and robust measures of nongaussianity that FastICA uses.
Finally, [12] is a textbook that summarizes and aggregates many of the
previous works listed above into a single reference. It devotes an entire chapter to
FastICA specifically. It also provides a proof of the convergence of the FastICA
algorithm and a detailed derivation of it.
4


2. Background Information
To understand and implement FastICA there is a significant amount of
background information that must be understood. This chapter will establish
the mathematical model used throughout this thesis and in the supplied source
code, as well as explain the preprocessing required for and the workings of the
FastICA algorithm.
2.1 Model
A suitable model for the cocktail-party problem is:
Xi(t) = ailSi(f) + CLi2S2(t)+ +dinSn(t)
%2 (t) = 02lSl(t) + ^22^2 {t)+ +Cl2nSn{t)
(2.1)
Emit') (^ml^l(t) T dm2$2 (^)~b T rann(^)'
Si(t),s2(t), ,sn(t) are the original sources and can be represented in vector
form as s = (si(t), s2(t), ,sn(t))T. Similarly, x\(t), x2(t), ,xm(t) are the
observations of the mixed original sources, and is x = (xi(t), x2(t), ,xm(t))T
in vector form. a,j, where i=l..m, j=l..n, are the mixing coefficients and A
is the corresponding m x n so called mixing matrix representation of them.
Therefore, n is defined as the actual number of independent sources, and m is
defined as the actual number of observations. Employing the matrix and vector
forms of the components of equation (2.1), namely A, s, and x, allows us to
5


represent equation (2.1) in this concise form:
x = As. (2.2)
To solve an ICA problem, we attempt to determine s. If W = A-1, then
clearly
s = Wx (2.3)
and the problem is solved. Therefore, ICA algorithms attempt to find the demix-
ing matrix W given only x, under the assumptions made in section 2.2. How-
ever, in later sections it will become clear that we are only estimating s, therefore
equation (2.3) should be modified to read s = Wx. Or if we define W to be the
matrix that gives an estimate of s, s = Wx.
As an example, if our fictitious cocktail party had 11 attendees, all of whom
were speaking, and 20 microphones placed throughout the room, then n would
be 11 and m would be 20.
2.2 Assumptions
As discussed in the introduction, it is possible to solve the cocktail-party
and other related problems given that a set of assumptions hold. In this section
we will discuss what these assumptions must be in order to effectively perform
ICA.
2.2.1 Independence
ICA algorithms assume that the original sources are statistically indepen-
dent. This assumption distinguishes independent component analysis from de-
6


pendent component analysis. The best results will be achieved when this as-
sumption holds, nevertheless, when sources are not statistically independent
some form separation is usually still possible.
2.2.2 Linearity
ICA solves problems which are accurately modeled by equation (2.2). The
sources in equation (2.1) are mixed in a linear manner. If a mixture of sources
is not linear then ICA may not be able to separate the sources. Therefore,
throughout this thesis all mixed sources are assumed to be mixed linearly in
accordance with the model presented in section 2.1.
2.2.3 Ambiguities
Several factors will influence the order in which sources are separated. For
instance, depending on the implementation, PCA, as described in section 2.3.2
may reorder sources. It is a basic assumption that the order of sources cannot
be determined by ICA and is therefore ambiguous.
Similarly, it is impossible to tell if the sign of a separated source is the same
as the original or inverted. It is assumed that this is an acceptable ambiguity
when performing ICA.
As an example, say that n is 4 and therefore s = s2{t), s3(f), s4(t))T.
ICA could yield an s that is this vector s = (s2(t), S4(t), s'i(t), s3(t))T. In
other words, the estimated independent sources or components of s are in a
different order and have arbitrary sign inversions with respect to the original
components of s.
Fortunately, with audio signals, as is the case in this thesis, the sign of the
resulting signals is not important, nor is the order.
7


2.2.4 Observations
ICA needs at least as many observation signals (m), as there are number
of sources (n) we would like to separate. If under-determined (m < n) it is
believed that ICA will still separate the n sources into m signals, but they
will not correspond exactly to the original n sources. If over-determined (m > n)
a method such as principal component analysis, as discussed in 2.3.2, can reduce
the dimensionality such that m = n.
2.2.5 Nongaussianity
Of all the assumptions, the assumption of the nongaussianity of sources is
the most important, i.e., if this assumption does not hold, ICA is not possible.
All the components of s, save one, must not be Gaussian[3]. Again in [3],
Hyvarinen aptly explains why this is. Given two Gaussian random variables X\
and x2, ...one can prove that the distribution of any orthogonal transformation
of the Gaussian (x1,x2) has exactly the same distribution as (X\,x2), and that
x\ and x2 are independent. Thus...we can only estimate the ICA model up to
an orthogonal transformation.
That the distribution of the sum of arbitrarily distributed independent ran-
dom variables tends to the Gaussian (under certain circumstances), is the well
known central limit theorem [3] [13]. Intuitively then, the distribution of the
sum of any two arbitrarily distributed random variables leads to a more Gaus-
sian distribution than either variable is on its own[3].
It is therefore desirable that sources be nongaussian. Nongaussianity is a
measure of how different a random variables distribution is from the Gaussian.
8


Now It should be clear that if we find an algorithm that calculates the demix-
ing matrix W as in equation (2.3) such that it maximizes the nongaussianity of
each source, we will have identify the individual independent components, not-
ing that a single source will be more nongaussian than any mixture of sources,
again by the central limit theorem. The focus then is to find an algorithm that
can maximize the nongaussianity of the mixed sources. FastICA is one such
algorithm.
Because nongaussianity plays such a vital role in ICA, we will discuss dif-
ferent ways to measure and approximate it in the following sections.
2.2.5.1 Kurtosis
Before learning how kurtosis allows us to measure the nongaussianity of a
source, we must have a very basic understanding of cumulants. For our pur-
poses, cumulants can be defined as statistical quantities [14] that characterize
distributions. When we speak of cumulants we must always specify of which or-
der we speak. Informally, the first-order cumulant is the mean of a distribution.
The second-order is the variance. The third-order is a property called skewness,
which tells us if a distribution is skewed to the left or right. The fourth-order is
called kurtosis, the topic of this section. Table 2.1 summarizes cumulants, their
definitions, and their interpretations.
Kurtosis is a measure of how peaked or flat a particular distribution is. The
Gaussian distribution has the special property that its fourth order cumulant
(kurtosis) is zero. For our purposes, we assume that the Gaussian distribution is
the only one in which this is the case, even though there are some rare exceptions.
A distribution that is more peaked than Gaussian will have a positive kurtosis,
9


Table 2.1: Cumulants up to and including the 4th order, where px E{x}
and x is a random variable that follows an arbitrary distribution.
Cumulant Definition
Interpretation
Cj{x} E{x}
Mean
Variance
Skewness
C2{x} E{(x fj,x)2}
C3{x} E{(x-fixf}
CA{x} E{(x-/ix)A-
3 [E{(x px)2}}2} Kurtosis
conversely, a flatter-than-Gaussian distribution will have a negative kurtosis.
So, by measuring the absolute value of kurtosis, we can measure how far a
distribution is from the Gaussian, i.e., its nongaussianity.
Although kurtosis is easily understood and provides us a natural measure
of nongaussianity, unfortunately, outliers in the data can severely distort it.
Therefore it is rarely used, in practice, for ICA[3]. Since kurtosis is not robust,
a different measure of nongaussianity is often used called negentropy.
2.2.5.2 Negentropy
For a discrete random variable X, entropy is defined [15] as
Extending to continuous random variables, and for the continuous random
variable X with support set S, and probability density function/(a:), the differ-
ential entropy is defined [15] as
(2.4)
(2.5)
10


Since the Gaussian distribution maximizes the entropy over all distribu-
tions with the same variance, [15] entropy can be used as a measure of non-
gaussianity. Usually this is done by measuring negentropy instead of entropy.
Negentropy is defined as
J(y) = H{ygauss) H(y)
, where ygaUss is a Gaussian random variable of the same covariance matrix as
y. [3] It should be clear that J(y) will only be zero when y is a Gaussian random
variable, and that the farther y is from a Gaussian random variable, the larger
the negentropy, and consequently the nongaussianity, will be.
Although negentropy is a good measure of nongaussianity, like kurtosis,
it is not without its problems. Computing entropy is difficult since doing so
relies on knowing a priori (or estimating) probability mass or probability density
functions as required by equations (2.4) and (2.5).
2.2.5.3 Approximating Negentropy
Approximations based on higher order cumulants have been developed to
estimate negentropy without the need of a probability density function. Re-
call, however, that measuring nongaussianity through kurtosis (the 4th order
cumulant) is subject to difficulties with outliers. High-order cumulant based
approximations of negentropy are also subject to the same problem. For this
reason, Hyvarinen developed a new approximation of negentropy in [6]. The end
result of [6] as stated in [3] is restated here for reference without development
11


or proof
J(V) oc {E{G(y)} E{G(y)}f. (2.6)
Noting that v is a zero mean, unit variance Gaussian variable, and y is also
a zero mean, unit variance random variable, but of arbitrary distribution. Fi-
nally, G(-) is nearly any non-quadratic function, and will be discussed in section
2.5. This approximation has proven itself to be robust and still approximates
negentropy (and therefore measures nongaussianity) sufficiently well.
2.3 Preprocessing
In order to make FastICA and other ICA algorithms tractable and to reduce
their computational costs, some preprocessing of the input data is usually done.
FastICA, specifically, performs centering, principal component analysis (PCA),
and whitening. Each of these procedures will be explained in the following
sections.
2.3.1 Centering
Centering simply refers to removing the mean from each observed source.
Throughout this thesis it is assumed that all observations are centered. For
example, to center the vector x0 as x, this operation: x = x0 F'jxo}, would
be performed. Centering the data simplifies the mathematics involved in ICA
and does not harm its ability to separate sources.
2.3.2 Principal Component Analysis
Principal component analysis is a method to reduce the dimensionality of a
dataset. PCA is well suited as a preprocessing step in ICA since often m > n
yet we usually desire that m = n. Reducing the dimensionality also leads to
12


fewer calculations in the later stages of ICA[2]. This is an import consideration
when n is large, or when we desire to perform ICA in real-time.
PCA is performed by finding the eigenvalue decomposition of the covariance
matrix of the data, i.e., the eigenvalue decomposition of
£{xxr}. (2.7)
Let V be the eigenvectors and D be the matrix containing the eigenvalues of
equation (2.7) along its diagonal and zeros elsewhere. The diagonal elements
of D are the variances of the components of x projected onto the correspond-
ing column vectors of V. Intuitively, the larger the variance (eigenvalue), the
more interesting a particular direction (eigenvector) is. To actually reduce
the dimensionality one chooses to keep the largest, or in some cases the small-
est, eigenvalues and corresponding eigenvectors, and discards the rest. Thus
choosing to keep only the l largest eigenvalues, where l < m, of D and the cor-
responding vectors in V is the process of reducing the dimensionality to l and
simultaneously choosing the most interesting, or principal components. The
l principal components can be chosen manually, or algorithmically.
For ICA, we desire that n = m. Usually, n eigenvalues will have similar
relatively high values, and the rest will be an order of magnitude or more smaller.
The number of relatively high valued eigenvalues is n. The m n small
valued eigenvalues are the ones chosen for elimination and thus dimensionality
is reduced to the point that n = m.
13


From an implementation point-of-view, note that the eigenvalue decompo-
sition employed for PCA may be one optimized for symmetric matrices since
equation (2.7) will always be symmetric. This can improve the performance of
the algorithm. See the jacobi() function in the ubiquitous [16].
2.3.3 Whitening
A whitened data set, say x, is one where the following property holds:
£{x?iT} = I.
That is, the covariance matrix is the identity matrix. This leads to the desirable
property that ...[xs] components are uncorrelated and their variances equal
unity [3].
There are a number of different methods to whiten data. This one is pre-
sented in [3], and is used in the provided code:
1. Find the eigenvalue decomposition of the covariance matrix of x, (i?{xxT} =
VDV-1, but JF{xxr} is symmetric so V is symmetric, and therefore
V-1 = VT[17]), where again, V is the orthogonal matrix containing the
eigenvectors (in columns) and D is a matrix whose elements dij = A,, when
i = j, and 0 elsewhere, assuming that Awith i = l..m, are the eigen-
values. Note that when implemented, E{-} must be replaced by sample
averages since probability mass and density functions are seldom known a
priori.
2. Find x, where x = VD"^VTx with D ^ being the same as D except that
_ 1
the diagonal elements are replaced by A- 2. for i = l..m.
14


As applied to ICA, whitening x to yield x leads to a modified model; instead of
identifying A in equation (2.2), we now search for B in x = Bs, where B = QA,
if we let Q = VD-IVT, noting that Q 1 = VD2VT. Again our new model
becomes
x = Bs. (2.8)
Now finding an approximation of the original sources (s) involves solving
this equation:
s = Brx.
Noting that it is shown in [2] that B is orthogonal, and therefore Br = B-1
[17], somewhat simplifying the problem. Once B is found, A can be computed
as Q-1B and W can be computed as BTQ. Then W can be used to demix the
original observation matrix x giving s. The algorithms presented below are in
terms of W and its components wi5 i = 1, 2,..., n, where each w* is a column
vector. In practice, we find B, however, it is easy to compute W from B.
Following in the example of [3], we will assume that x is whitened to simplify
the notation.
2.4 FastICA Algorithm
Without further delay, we introduce two FastICA algorithms. The first
separates components one at a time, the other separates all components simul-
taneously, and is the method used in the provided implementation.
2.4.1 Single Component FastICA
This version of FastICA finds one component of W, w, at a time, where
w is a column of W. It is often referred to as the deflation approach. The
15


algorithm is:
1. Initialize w to random normalized values.
2. w0id = w
3. w+ = E{xg{wTx)} E{g'{wTx)}w
4. w = w+/||w+||
5- If (! lwLwl) < e is true, goto 8
6. If maximum number of iterations is reached, goto 8
7. Goto 2
8. Done.
Now follows a discussion of each step of the algorithm in detail. Step 1 initializes
w and ensures that it is normalized, for reasons that will become apparent. Step
2 caches the old (previous) value of w so that we can check for convergence
in step 5. Step 3 is the heart of FastICA. This step is the fixed point iteration
that maximizes the nongaussianity of w through a modified approximation of
negentropy, as in equation (2.6). Step 4 simply normalizes w+. Doing so makes
the convergence detection easier, and also emphasizes that it is the direction
of w that matters. Since both wm and w are normalized the maximum their
absolute dot product will be is 1. Further, an absolute dot-product of 1 (within
the small positive threshold e) signifies that the two vectors are aligned and have
therefore converged. Steps 6, 7, and 8 handle the branching, termination, and
avoidance of infinite loops.
16


The updating rule in step 3 will now be derived by following closely the
example of Hyvarinen in [5]. Recall that if nongaussianity is maximized, the
term wTx will represent one of the independent components. Nongaussianity
can be estimated by equation (2.6). So as in [5], for one component:
Jg{w) = [£{G'(wTx} E{G(v)}]2.
Next, ...note that the maxima of Jq{w) are obtained at certain optima of
£,{G(wTx)} [5]. Since FastICA assumes whitened data, ||w||2 = 1, and it can
be shown that the optima of E{G(wTx)} are found when
£{xp(wTx)} (3w = 0 (2.9)
assuming that w0 is the the optimum value, j3 is a constant and can be calculated
as [3 = E{\Vq xg(w^x)}, and g(-) is the derivative of G(-). Employing Newtons
method to solve equation (2.9), the Jacobian matrix JF(w) is
JF{w) = E{xxV(wtx)} pi. (2.10)
Newtons method requires a matrix inversion at every step, so to ease its compu-
tational cost the first term in equation (2.10) is approximated by E{(/(wTx)}I
leading to a diagonal, and thus easily invertible, matrix. Further, since w0 is
unknown in p, it is approximated by the most recent value of w. Employing
17


the approximations, the same as in [5], we arrive at this Newton iteration:
w+ = w [£'{x5r(wrx)} Pw]/[E{g'(wT-x.)} j3\
w = w+/||w+||.
Multiplying by 0 E{g'{wTx)} results in the simpler fixed-point algorithm:
w+ = £'{x^(wTx)} E{g'{ wrx)}w
w =w+/||w+||
which is steps 3 and 4 of the algorithm.
See [12] for a similar derivation, and [5] and [12] for proofs of convergence.
Running the above algorithm will, again, result in only one component of
W. In order to find more, or all, the components, the algorithm must be run
repeatedly. However, to keep the result of each run from converging to the same
fixed point as previous runs, we must decorrelate the outputs ... after every
iteration [5].
The Gram-Schmidt orthogonalization process proves useful for this purpose.
Assuming that p components are desired, the Gram-Schmidt process works by
subtracting all previous projections from the current. Doing so ensures that com-
ponents are orthogonal to one another. From [3] [5] the Gram-Schmidt method
of orthogonalization is:
Wp+l = wp+1 - wJ+iWjWi
Wp+1 = Wp+1/^wj+1 Wp+1.
18


2.4.2 Symmetric FastICA
Now that single component FastICA is understood, the algorithm can be
extended to multiple simultaneously obtained components, often known as the
symmetric approach.
1. Initialize each column of W to random normalized values, then orthogo-
nalize it.
2. WoW = W
3. For each i = 1,2, ...,m compute w+ = E{xg(wfx)} E{g'(wjx)}wi,
resulting in W+
4. W = ((W+)(W+)r)(-1/2)(W+)
5. If (1 min(|diag(WTW0/d)|)) < e is true, goto 8
6. If maximum number of iterations is reached, goto 8
7. Goto 2
8. Done.
The algorithm is very similar to the single component one, but a few steps
require some explanation. Step 4 extends normalization to a matrix. This
process is called symmetric orthogonalization. It requires the square root of a
matrix which is explained later in section 3.3.1. Step 5 extends the simple single
component convergence test to all the components in W. More specifically,
finding the minimum absolute value of the diagonal elements of WTWoW ensures
19


that all components of W have converged sufficiently (to within the value e).
To reiterate, this algorithm finds all m components simultaneously.
2.5 Selection of G(-)
As mentioned in 2.2.5.3, the function G(-) is almost any non-quadratic func-
tion. This section will briefly list common choices for G(-) and their specialized
uses as discussed in [5].
Gi{u) = jMogcosh(aiit)
g\{u) tanh(oiu)
G2(u) = A exp(-o2u2/2)
g2{u) =uexp(-a2u2/2)
(2.11)
(2.12)
G3{u) = \u4
g3(u) = u3
where 1 < a\ < 2 and a2 ~ 1 [5].
Equation set (2.11) has generally good performance properties, however if
the signals in use are highly super-Gaussian, equation set (2.12) may be better[5].
It can be shown that equation set (2.13) is equivalent to using kurtosis, so its
use suffers from problems with robustness unless the signals in use are sub-
Gaussian[5j. Nevertheless, it is straight forward to compute and was therefore
used in the implementation of FastICA that accompanies this thesis.
Density functions are rarely known in advance, but it should be noted that
if known, the optimal function G(-) for a given situation can be determined
given the density function of the source being estimated. If f{u) is the density
20


function, then Gopt(u), the optimal function, is Gopt(u) = log f(u).
21


3. Implementation
The author of this thesis implemented a real-time blind source separation
system for audio using the FastICA algorithm. The system is made up of several
logical application programs that work in conjunction to achieve the end result.
Indeed, each logical application is truly modular and only performs one piece
of the work required to simulate mixed signals, record signals, and separate
them into independent sources. The real-time audio data flows from program
to program in a sink/source type system. This is facilitated by the use of the
JACK audio connection kit, or simply JACK.
Each logical application program is actually divided into two real applica-
tion programs, one written in C, known as the back-end, and the other written
in Python, known as the front-end. The C-language applications do the actual
digital signal processing, interface with JACK, provide a simple socket based
control interface, and respond to stimuli from the Python front-end programs.
The Python application programs utilize the PyGTK Python module to build
simple GUIs which control (through a socket) the back-end application pro-
grams.
3.1 JACK
The JACK audio connection kit is a low-latency sound server for the
GNU/Linux and Mac OS X operating systems. Its main feature and virtue
is the ability to stream real-time audio to, from, and through separate audio
application programs. This leads to a very flexible sink/source/filter design.
22


Hardware sources and applications that synthesize sound serve audio to filter
applications which in turn feed any number of applications which eventually
feed the hardware sinks (loudspeakers) or recording applications. Connections
between applications are dynamic and can be created or destroyed at run-time.
As an example, think of a audio synthesizer that creates sound effects. It
is often desirable to filter sound effects by amplifying or attenuating certain fre-
quency ranges. If such a filter program existed, the synthesizer output could be
linked to the filter programs input at run-time. The output of the filter could
then be connected to the hardware (loudspeaker) sink device. At this point,
one would have a filtered synthesizer sound being output over the loudspeaker.
Suppose that a recording of the filtered audio is also desired. A recording ap-
plication could be started and the filter output added to its input. Figure 3.1
illustrates a block diagram of these example connections.
3.2 Logical Application Programs
In the authors real-time FastICA implementation there are 5 logical appli-
cations that work together to perform BSS. They are jackfilestream, jackconnect,
jackrecord, jackica. and jackmix. jackfilestream. jackmix, jackrecord, and jack-
connect are convenience applications used to setup an environment suitable for
jackica to do the demixing. Each application will be described here in some
detail.
jackica. jackica is the application that runs the FastICA algorithm to
perform BSS of the audio signals. It takes m mixed audio inputs and
produces the n demixed audio outputs. As currently implemented, jackica
23


Figure 3.1: An example of the connections possible in JACK.
assumes that m = n, though this could be changed with some additional
effort. The m audio inputs are assumed to be linear mixtures of the n
original audio sources. In real-time, at a sample rate of 8192Hz, jackica
reads the mixed audio data, demixes it using the FastICA algorithm, and
writes it to its outputs.
jackfilestream. jaekfilestreain is a convenience application that streams,
in real-time, an arbitrary number of audio files to its outputs. From the
point-of-view of other applications, jackfilestream appears to be an array
24


of microphone outputs. Indeed, in JACK there is no difference between
an application synthesizing, say 15, audio streams, and actually having 15
microphones streaming live data. In other words, the applications that
connect to the hypothetical 15 streams are completely unaware if the up-
stream audio is being synthesized or if it is from a real microphone. In
practice, jackfilestream is used to stream the n unmixed audio sources in a
controlled, reproducible environment, however, if live audio (microphones)
were used nothing in the rest of the system would have to change.
jackmix. jackica assumes that its audio inputs are linear mixtures of the
original n sources. A mixer is required to simulate the m observed audio
sources since jackfilestream does nothing but stream audio sources. Ob-
viously in a real-life experiment utilizing microphones, jackmix would not
be necessary, jackmix provides a way to mix the n sources in a completely
random way, or to choose from 3 preset mixtures. The preset mixtures,
again, facilitate the reproducibility of experiments. A GUI displays the
components of A, which is the mixing matrix as defined in equation (2.2).
Summarizing, jackmix takes n audio inputs (usually from jackfilestream)
and mixes them according to a preset or random mixing matrix and out-
puts m audio sources. Again, m = n is assumed, but with a small effort
this could be changed to allow for m to be greater than or equal to n.
jackrecord. jackrecord is a simple audio recording application. It is used
to record the audio pipeline at any stage so that the real-time signals can
be analyzed off-line. Multiple instances of jackrecord can be run simulta-
25


neously and all of the recordings from each instance can be synchronized
through JACKS transport control API. In other words, every jackrecord
instance will start its respective recording at the same time as every other
instance making comparisons between recordings trivial.
jackconnect. JACK is very flexible in that it allows for arbitrary con-
nects between the inputs and outputs of audio applications. At times
this flexibility is difficult and tedious to manage, jackconnect manages
some of the complexity by automatically connecting jackfilestream, jack-
mix, jackrecord, and jackica in the expected manner which is shown in
figure 3.2. For performance reasons, m and n where kept to 4 as in the
figure, but this can be changed at compile time, jackrecord is strictly op-
tional, and as already mentioned, can be made to record any part of the
pipeline. In the figure, si, s2, s3, and s4 signify Si(t), s2(t), s3(t), s4(f). xl,
x2, x3, x4 signify x\(t), x2(t), x3(t), x{t), and sl, s2, s3, and s4 signify
Usually only one element of s will be output to a
speaker at a time. This is the case in figure 3.2 where s2(t) is the one
being output. Figure 3.3 is a screenshot of the actual JACK connections,
as presented by qjackctl, that are required to make the jackutils/fastica.c
system simulate and perform ICA on the audio data, qjackctl is an open
source application that can display, create, and destroy arbitrary connec-
tions between JACK applications.
The logical applications are divided into two real applications, the back-end
and front-end.
26


Figure 3.2: The expected connections between jackfilestream, jackmix, jackica,
and jackrecord.
3.2.1 Back-end Application Programs
The back-end applications are written in C for performance reasons. They
depend on libjack, the library portion of JACK, libjack is released under the
LGPL license which allows non open source applications to link with it without
the need to open source themselves. Additionally, jackica depends on a new
library fastica.c which is explained in section 3.3. jackfilestream depends on
libsndfile in order to read several popular audio file formats, libsndfile is also
an LGPL licensed library. The fastica_c library is, at this time, proprietary
but listed in appendix A for reference purposes. The applications connect to
27


Figure 3.3: Screenshot of the JACK connections between jackutils applications,
as presented by qjackctl, that are needed for BSS of audio signals.
JACK and then listen on a BSD style socket for commands. Each command is
processed in turn to perform the requested actions.
3.2.2 Front-end Application Programs
The front-end applications are written in Python and depend on the PyGTK
package for their graphical user interfaces. Python was chosen because it is
simple, yet powerful, and lends itself to rapid GUI development with PyGTK.
These front-end applications draw their GUIs and then send commands over
BSD style sockets to back-end applications in response to user interaction. See
figure 3.4 for a screenshot of the front-end applications.
28


Figure 3.4: Screenshot of PyGTK GUIs for front-end jackutils applications
and qjackctl.
3.3 Implementation Details fastica_c
Implementation details of the back-end applications are presented in this
section, jackica requires calls into the new fastica.c library. fastica_c is divided
into two main sections, general matrix handling functions, and FastICA specific
functions. The FastICA specific functions make heavy use of the general matrix
functions.
fastica.c is very modular and could be easily ported to and used in systems
requiring matrix manipulations but not ICA. Its implementation uses the soft-
ware engineering practice design by contract. Meaning that before a function will
29


perform any actions it checks to make sure that all of its preconditions are met.
If they are not, an error is returned before any action takes place. A set of C lan-
guage macros was developed to assist in the sometimes tedious task of checking
every precondition. See MATRIX_VERIFY_PRECOND(), MATRIX_VERIFY_GOTO0,
and MATRIX_VERIFY_SET_GOTO() in matrix.h in appendix A. While adding pre-
condition checks to every function might seem to slow development, the cost
was well justified since the precondition checks caught many non-obvious and
difficult-to-find software bugs in early testing of fastica_c.
fastica_c functions will, for the most part, fit on oik; screen of 25 lines of text.
This is by design, since smaller functions are often easier to debug, understand,
and maintain.
fastica_c also makes heavy use of goto, goto has often been criticized as
leading to unmaintainable code. However, in fastica.c goto is used strictly as a
way to avoid deeply nested if statements, which in the opinion of the author,
are far more difficult to maintain than well constructed functions that utilize
goto.
Most functions in fasticam have this general structure:
1. Assert statements. Passing a NULL matrix_t pointer into a fastica_c is
considered a major violation and will therefore cause an assertion to fail
leading to the demise of the application using fastica_c.
2. Precondition checks. As discussed above, checks guarantee that all pre-
conditions (or assumptions) that a particular function is making are valid.
If any are not valid an error is returned immediately.
30


3. Action statements. Once the preconditions are met the desired actions
can be performed. If a function has many preconditions, a static helper
function might be called to do the actions instead. These helper functions,
for performance reasons, do not check any parameters but facilitate in
keeping functions to less than 25 lines of text. The names of these functions
usually begin with an underscore^). Once again, goto statements are used
liberally to avoid nested if statements.
4. Cleanup statements. This part of a typical fastica.c function involves the
cleanup of temporary storage and final processing of the return value.
Some good examples of this general structure can be found in matrix_cp() in
cp.c, matrix_eig() in eig.c, and matrix_ica() in ica.c.
3.3.1 General Matrix Functions
At the heart of fastica.c is the matrixjt type, matrixjt is a simple container
structure for a 2D array. It is made up of three fields matrix, nrow, and ncol.
matrix is the 2D array and nrow, and ncol are the dimensions of matrix. These
are stored so that runtime precondition checks can be performed.
matrixjt type matrices must be allocated before use. matrixmewO allo-
cates a matrix of the supplied dimensions on the heap. Once a matrix is no
longer needed matrix_free() returns the memory to the heap. Heap alloca-
tions, while very user friendly, are notoriously bad for real-time performance
because allocation times can vary drastically given the state of the heap at the
time of allocation or deallocation.
31


Table 3.1: Average elapsed time allocating and deallocating matrices over 3
trials of 5000 iterations each.
Dimensions Heap Pool
500x500 108928 ms 1 ms
8192x4 3615 ms 1 ms
4x8192 8 ms 1 ms
Often it is possible to know how many matrices of a particular dimension
are required to solve a problem a priori. Certainly that is the case with fas-
tica_c. In such situations, a matrix pool can drastically reduce the allocation
and deallocation times. A matrix pool works by preallocating all the matrices
a program or problem might require at any given time. The operation of al-
locating a matrix reduces to finding a free matrix in the pool, and freeing a
matrix reduces to returning it to the pool. When there are a very large number
of allocations and deallocations in a short period of time, a matrix pool can
improve performance dramatically as can be seen in table 3.1. Each time value
represents the average number of milliseconds it took to execute matrix_new()
and matrix_free() 5000 times over 3 trials on the authors system. It should
be clear that matrix pools lead to nearly constant allocation and deallocation
times no matter the dimensions, while heap allocations/deallocations are much
slower in general, but are particularly slower the larger the number of rows, fas-
tica_c enables and disables matrix pools through the calls matrix_pool_init ()
and matrix_pool_destroy(). See appendix A for details, jackica makes use
of matrix pools to avoid allocation and deallocation performance problems in
real-time.
32


It may be difficult to calculate the number of matrices of each dimension
that are needed a priori. For this reason, fastica_c tracks the number of alloca-
tions made with matrix_new() and provides the function matrix_trace_dump()
which reports the number of allocations broken down by dimensions.
fastica_c requires the use of matrix.t type matrices, however, tradition-
ally matrices in C are implemented as simple 2D arrays. To avoid unnec-
essary copying from one type to the other, another set of allocation and
deallocation functions are provided namely matrixmew_from_buffer() and
matrix_free_from_buffer(). matrix_new_from_buf f er () takes a standard 2D
array and converts it into a matrix_t type matrix without copying.
matrix_free_from_buffer() clearly, frees such an allocated matrix, jackica
employs matrix_new_from_buff er0 and matrixJree_from.buffer() to share
buffers between JACK and fastica_c without copying.
Once a matrix.t is allocated it can participate in several matrix opera-
tions provided by fastica.c. Copying, transposing, multiplying, adding, sub-
tracting, elementwise-multiplication, and elementwise-division, are provided by
matrix_cp(), matrix_transpose(), matrix_mult(), matrix_add(), matrix_sub(),
matrix_elem_mul(), and matrix_elem_div() respectively. In addition to these
simple functions, functions for computing the mean and variance of vectors,
covariance, square root, determinant, inverse, and eigenvalue decomposition
(jacobi()) are provided. It should be noted that matrix inversion, although
provided, is not required for the FastICA algorithm.
Of special interest is matrix.sqrt (). It computes the so called square root
of a semi-positive definite symmetric matrix by way of the eigenvalue decom-
33


position. Doing so for a matrix say Z, means finding another square matrix P
such that Z = P P. P is found by decomposing Z into V and D by eigenvalue
decomposition, then finding the square root of each eigenvalue in D. Therefore,
P = VD2 V-1 if Z = VDV1, in the general case, and P = VD 2 VT since Z is
symmetric.
All of the source code for fastica_c save one file, jacobi.c (taken from [16]),
is the sole work of the author of this thesis. The main function in jacobi.c is
jacobi(), the eigenvalue decomposition of a symmetric matrix. Since permission
to list jacobi.c was not granted nor requested, it is not listed in appendix A.
Free replacements for jacobi() do exist but have not been tested with the rest
of fastica.c and jackutils.
3.3.2 FastICA Specific Functions
The ICA specific portions of fastica.c are based heavily on Aapo Hyvarinens
free FastICA MATLAB source code [1] [18]. matrix_pca() implements principal
component analysis as discussed in section 2.3.2. matrix_whiten() whitens
a dataset as well as provides a whitening matrix (Q, as in section 2.3.3)
and dewhitening matrix (Q~\ also as in section 2.3.3) for later use. Finally,
matrix_ica() performs the FastICA algorithm.
For this thesis, equation (2.13) was chosen as the function for G(-) because
of its straight forward and simple implementation, despite the fact that kurtosis
based approaches can suffer from robustness issues. The symmetric approach, as
discussed in section 2.4.2, as opposed to the deflation approach, was implemented
since it potentially lends itself to more parallel processing. The desire for highly
parallel processing was a primary concern to the author when first implementing
34


fastica_c since a GPU, or graphics processing unit, was originally going to be
used for the calculations. This plan was abandoned when it was determined
that a modern PC was faster at separating 4 independent sources than a GPU.
However, if more than say 20 sources are to be separated, a GPU implementation
would very likely be faster.
Finally, fastica_decompose_sources() is a convenience function that en-
capsulates all of the preprocessing and ICA specific calls into an easy to use
API. jackica uses this function to separate sources in real-time.
3.4 Implementation Details jackutils
jackfilestreain. jackmix, jackica, jackrecord, and jaekconnect are collectively
known as jackutils. Appendix B contains the complete source of jackutils. Of
these the most important is obviously jackica.
JACK audio applications are built on a callback model. When the JACK
audio server, known as jackd, is ready for a particular application to process
audio data, it will call that applications registered callback (known as its pro-
cess() callback). The processQ callback is given a fixed number of audio samples
to process and output, jackutils is designed to work at a sample rate of 8192Hz,
with one second of data being processed per process() callback.
FastICA requires long sequences of historical data to work well, therefore,
jackicas processQ callback buffers the last 6 seconds worth of audio input and
calls FastICA on all of it, of which only the last one second is output. A circular
buffer is utilized to maintain the last 6 seconds.
6 seconds was chosen empirically, since any amount of audio data smaller
than that suffered from a problem where components would suddenly change
35


the order in which they were output from the FastICA algorithm. This usually
happened when one of the n components in the mixtures went to silence briefly.
When one component is silent n effectively changes from 4 to 3 and an essentially
new problem is solved. By processing more data per process() callback, the
length of the silence becomes small when compared to the 6 seconds worth of
audio data and therefore no order change results. This order change problem
can be seen as a limitation of FastICA. FastICA assumes stationary sources,
but in general audio signals are not stationary processes. To keep output orders
from changing more data can be processed per process() callback, and it will be
less likely that order changes will occur. This resilience to order changes comes
at the price of more computations per period, and slower response times to long
lasting changes in source statistics.
36


4. Analysis
fastica_cs goal is to decompose and playback the independent components
of mixed signals in real-time with negligible delay. There are qualitative and
quantitative methods for determining if fastica_cs goal has been reached. The
analysis will show that fastica.c does a reasonable job of attaining its goal.
4.1 Qualitative Analysis
There are two main forms of qualitative analysis that will be used in this
section, auditory and visual.
4.1.1 Auditory Analysis
There are many psychological factors outside the scope of this thesis that
contribute to human perception of audio signals. A study unto itself would be
required to determine if audio signals separated by FastICA do indeed sound like
the original. Suffice it to say that, in the authors opinion, the test signals, when
scaled properly, mixed according to the canned mixtures of jackmix, demixed
with jackica, and output over conventional loudspeakers, are indistinguishable
from the originals.
JACK requires that all audio signals be scaled between -1 and 1 before
playback, or else clipping will result. Early versions of fastica_c suffered from a
problem where the output from jackica was in the range of -10 to 10. Even with
severe clipping, the audio signals were unmistakably identifiable, at least by the
author.
4.1.2 Visual Analysis
37


Table 4.1: Mapping between s and s.
s(l) PS -5(1)
s(2) -5(4)
s(3) rs i(2)
s(4) H3)
Modern systems allow us to see audio by way of computer graphics and
audio software. By comparing visual representations of recordings we have a
reliable, even if not scientific, method for determining how well fastica_c sepa-
rates sources. As an example, some real audio recorded by jackrecord is shown
in figures 4.1 and 4.2. They show 8, 5-second signals in all. The first signal
(cleanTJ. or s(l)), and every other after it that begin with clean_, represent
the original sources or s(l), s(2), s(3), and s(4) in our general model equation
(2.1). Immediately below each clean signal is its corresponding output signal
from jackica, that begins with demixedJ and are also referred to as s(l), s(2),
s(3), and s(4). Notice that s(l) and s(2)s corresponding output signals are
flipped with respect to the originals.
By visual inspection, it should be clear how the author was able to identify
which input signals went with which output and if an output was flipped with
respect to the original or not. This leads to the mapping in table 4.1. The
signals are very similar and seemingly only differ by a scaling factor which could
be negative. Recall that section 2.2.3 explains that this is one of the ambiguities
inherent to FastICA.
4.2 Quantitative Analysis
38


Figure 4.1: Comparison of the first two components of s with corresponding
components of s. 5 seconds of each waveform is displayed.
Qualitative observations led us to table 4.1, but based only on subjectivity.
We need quantitative support that table 4.1 is correct. We will use a measure
called correlation to do so.
For our purposes, given random variables X and Y, with means fix and //y
respectively, covariance is defined as
cov(X, Y) = Y, ^ ~ '
i=1
39


Figure 4.2: Comparison of the second two components of s with corresponding
components of s. 5 seconds of each waveform is displayed.
Covariance can also be denoted as ctxy, be., cov(Al, Y) = axr- Utilizing the
new notation, correlation between X and Y is defined as
Qxy
yJOXXOYY
Correlation will be specifically used as a measure of how similar two signals are
to one another. A positive correlation near 1 means that two signals are very
similar in the sense that as one rises the other does so also, and as one falls so
does the other. Correlation near 0 means that two signals are vastly dissimilar.
A negative correlation near -1 means that the signals are very similar, but are
cor(Af, Y)
40


Table 4.2: Correlation coefficients between components of s and s.
s(l) s(2) (3) *(4)
s(l) -0.9998 -0.0055 0.0072 -0.0073
s(2) -0.0046 0.0134 0.9991 0.0011
(3) -0.0071 0.0085 -0.0135 0.9999
s(4) 0.0098 -0.9999 0.0011 -0.0099
still different by a factor of -1, meaning that if one signal goes up the other
goes down and vice versa. Now given the definition of correlation, we should
be able to roughly predict what the correlation between matched original and
decomposed sources should be.
Since s(l) looks to be the mirror image of s(l) we should expect a strong
negative correlation. The same should be true for s(2) and s(4). s(3) and s(2)
rise and fall in step with one another so we should expect a strong positive
correlation. Again, the same should be true for s(4) and s(3). All 16 actual
correlations between s and s are computed in table 4.2. The guesses based
on qualitative visual observations seem to hold. Our visually mapped signals
have correlations -0.9998, -0.9999, 0.9991, and 0.9999. As expected the first
two are strongly negatively correlated while the last two are strongly positively
correlated. Figure 4.3 (a graph of table 4.2) makes it very clear how conclusive
table 4.2 is.
It is not enough to show that we can visually and quantitatively determine
which original sources go with the decomposed sources, we must also show that
the components in the assumed mappings have stronger correlations than the
decomposed signals as compared to all the mixed signals. Table 4.3 shows that
41


Table 4.3: Correlation coefficients between components of s and x.
x(l) x(2) x(3) x(4)
5(1) 0.8677 -0.8716 0.7225 -0.8842
8(2) 0.2826 0.4700 -0.3262 -0.1631
8(3) 0.4035 0.0855 0.6087 0.4156
5(4) 0.0327 0.1099 0.0176 0.1359
this is indeed the case. In other words, none of the correlations in table 4.3
are as strong as -0.9998, -0.9999, 0.9991, and 0.9999, which are the correlations
between mapped signals. Figure 4.4 is the waveforms of the mixed original
signals or x. Comparing it by visual inspection to figures 4.1 and 4.2 shows why
this result is not surprising; the signals truly are mixtures and do not strongly
resemble any original source in s.
The MATLAB source code used to generate the correlations in tables 4.2,
4.3 and figure 4.3 can be found in appendix C.
Given the qualitative and quantitative results, it is clear that the FastICA
is working and that the linearly mixed original signals are recoverable to within
a scale factor.
42


Figure 4.3: Correlation coefficients between components of s and s.
43


Figure 4.4:
Waveform of 5 seconds of mixed sources.
44


5. Conclusion
A detailed explanation of the cocktail-party problem, general blind source
separation, independent component analysis, the fast FastICA algorithm, and
finally a real-time implementation and analysis of the FastICA algorithm as
applied to blind source separation of audio signals were presented in this the-
sis. Several pieces of supporting theory, including principal component analysis,
whitening, centering, negentropy, cumulants, kurtosis, nongaussianity, and as-
sumptions for use of ICA, were also presented.
5.1 Analysis of Implementation
The analysis of the real-time implementation of FastICA strongly indicated
that source separation was successful and that real-time ICA is practical on
todays hardware. The analysis of the implementation used a combination of
qualitative and supporting quantitative metrics (correlations of nearly 1 and -1)
to finally conclude that source separation was indeed happening.
5.2 Future Research
One of the most challenging aspects of real-time FastICA is its need to pro-
cess relatively large amounts of historical data. Brute force calculation proved
to work in this implementation, but more efficient methods for dealing with it
should be researched.
The initial intent was to implement FastICA on a GPU utilizing a language
such as Brook, Brook+, OpenCL, or the CUDA framework, not a general pur-
pose CPU. Such languages and frameworks make it easy to exploit the massive
45


parallelism that GPUs provide. However, initial testing indicated that to over-
come the high cost of transferring data between the CPU and GPU many more
sources, say greater than 20, were required. The main purpose of this thesis was
real-time performance of BSS for audio signals, so 3 to 8 signals were sufficient
so the GPU was abandoned. Nevertheless, many of the preprocessing steps of
FastICA and maybe even the fixed point iteration itself could be tailored to
exploit the parallelism provided by GPUs. More research should be done in this
area. Similarly research into DSP and FPGA based implementations would also
be beneficial. Generally, where to divide the work between a main processor
and a specialized one (DSP, FPGA, GPU, etc.) would be a challenging research
topic.
FastICA could be extended to video source separation. The current FastICA
implementation could be used if the 2D video signal in question is somehow
pieced together into a ID one. However, keeping video in its 2D form should
also be investigated.
5.3 Main Conclusions
The main conclusion of this thesis is that under the assumptions presented,
FastICA is a viable option for real-time blind source separation of audio signals.
Some secondary conclusions are that modern off-the-shelf hardware is sufficient
to perform real-time ICA. Conversely, FastICA would likely benefit from special-
ized hardware like DSPs, FPGAs, cell processors, and GPUs, but more research
is needed. Another secondary conclusion is that design-by-contract can be useful
for catching software bugs early in the development process.
46


APPENDIX A. Source Code fastica_c
This appendix contains the source code for fastica_c.
A.l fastica.c/api.c
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
^include
#include
#include
#include
#include
#include
^include
#include
^include
< s t d i o.h>

matrix h
whiten h
mult. h
mean. h
ica h
pea.h
cp.h
#include fastica. h
int verbose = 0;
/ *
* fastica.decompose.sources ()
*
* Easy to use wrapper for matrix_ica() et al .
* Applications will most likely want to use this directly.
*
* mix mixed sources (each source is a row vector)
* demix independent components are placed into each row of demix
* maxiter maximum number of iterations to perform
* actualiter number of iterations actually performed or 1 if didnt
converge
* numsources number of sources to keep.
* A negative number means keep | numsources | but keep those with
the least power
* A positive number means keep numsources but keep those with
the most power
*
*/
int fastica_decompose_sources (matrix_t *mix, matrix_t *demix ,
int maxiter, int ^actualiter ,
int numsources)
{
int ret ;
matrix.t A, W;
matrix.t E, D; /* E eigenvectors D eigenvalues */
matrix.t s; /* sources */
matrix.t wm, dwm, ws; /* wm whitening matrix, dwm dewhitening matrix, ws
whitened source */
matrix.t means; /* means of each row in s */
matrix.t restored.s ; /* restored sources */
47


MATRIXVERIFYERECOND ((mix != NULL) &&
(demix != NULL) -1) ;
MATRIXVERIFYERECOND( (mix>nrow = demix>nrow) SiSi
(mix>ncol = demix>ncol) 1);
MATRIXVERIFYERECOND( abs ( numsources ) <= mix>nrow 1);
MATRIX_VERIFY_PRECOND( matrix.new (&s mix>nrow mix>ncol) =
MATRIXERROR-NONE, -1);
MATRIXVERIFYERECOND) mat rix _cp ( mix &s ) =
MATRIX ERROR-NONE, -1);
MATRJX_VERIFY_PRECOND( matrix.new (&ws s nrow s.ncol) =
MATRIXERROR-NONE, -1);
MATRIXVERIFYERECOND) matrix.new(&wm, s nrow s.nrow) =
MATRIX_ERROR_NONE, -1);
MATEDC_VERIFYEEECOND( matrix.new(&dwm, s.nrow, s.nrow) =
MATRIXERROR-NONE, -1);
MATRIXVERIFYERECOND( matrix_new(&E, s.nrow, numsources)
MATRIX ERROR. N( )\K. -1);
MATRIXVERIFYERECOND) matrix.new(&D, s.nrow, 1) =
MATRIXERROR-NONE, -1);
MATRIXVERIFYEEECOND( matrix.new (fcmeans s.nrow, 1) =
MATRIXERROR-NONE, -1);
MATRDC_VERIFYERECOND( matrix.new (&A, s.nrow, s.nrow) =
MATRIXERROR-NONE, -1);
MATRIXVERIFYERECOND ( matrix .new (&W, s.nrow, s.nrow) =
MATRIXERROR-NONE. -1);
MATRIXVERIFYERECOND( matrix_new(&restored_s s.nrow, s.ncol) =
MATRIXERROR-NONE, -1);
/* remove mean */
ret = matrix_remove.mean_transpose(&s femeans) ;
/* pea */
ret = matrix_pca(&s &E, &D, numsources);
/* whiten */
ret = matrix.whiten(&s &D, &E, &ws, &wm, &dwm) ;
/* ica */
ret = matrix_ica(&ws &wm, &dwm, maxiter MATRIXJCAEOW3, &A, MV, actualiter
);
/* restore sources */
ret = matrix.mult (MV, Sis, Sit es t o r e d _s ) ;
/ restore mean */
ret = matrix_restore_mean_transpose (Ssmeans Sires t or ed _s ) ;
ret = matrix_cp(&restored.s demix);
MATRIXVERIFYERECOND( matrix_free(&s)
MATRIXVERIFYERECOND ( matrix_free(&ws)
MATRIXVERIFYERECOND ( matrix_free(&wm)
MATRLX.VERIFYERECONDf matrix.free (&dwm)
MATRIXVERIFY_PRECOND( mat rix .free (&E)
MATRIXVERIFYERECOND( mat rix.free (&D)
MATRIXVERIFYERECOND) mat rix.free (fcmeans)
MATRIXVERIFYERECOND( mat r ix .free (&A)
MATRIXVERIFYERECOND( mat r ix .free (MV)
MATRIXVERIFYERECOND) matrix_free(&restored.s )
= MATRIXERROR-NONE,
= MATRIXERROR_NONE,
= MATRIXERROR-NONE.
= MATRIXERROR-NONE,
= MATRIXERROR-NONE,
= M ATRJXERROR-NONE.
= MATRDCERROFLNONE,
= MATRIXERROR-NONE,
= MATRIXERROR-NONE,
= MATRIXERROR-NONE,
-i);
-i);
-i);
-i);
-i);
-i);
-i);
-i);
-i);
-i);
return ret ;
}
A.2 fastica_c/cov.c
/a*****************************************************************
* Copyright (c) 2008 Anthony R. Jackson
*
48


* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
#include matrix, h
#include mean.h
#include var.h
/* computes the covariance between two column vectors (a, b ) and puts
result into result */
static int vector.cov (matrix.t *a, matrix.t *b, real_t ^result unsigned char
n_not.n_minus.one)
{
unsigned int r ;
r e a 1 t. sum ;
int ret ;
/* statically allocate lxl matrices for mean of a and b */
r e a 1 _t mu.a.r ;
r e a 1 -1 *mu_a_p ;
matrix.t mu.a;
real.t mu.b.r;
real.t tmu.b.p;
matrix.t mu.b;
/* asserts */
MATRIX ASSERT ( a != NULL);
MATRECASSERT(b != NULL);
MATRIX ASSERT (result = NULL) ;
/* preconditions */
MATRIX_VERIFY_PRECOND( (a>nrow != 0) &&:
(a> n c o 1 = 1) &&
(b>nrow != 0)
(b> n c o1 = 1) ,
MATRIXJ1RRORJNVALIDDIMENSIONS) ;
/* statically allocate 2 lxl matrices */
mu.a nrow = 1;
mu.a ncol = 1;
mu.a.p = &mu.a.r ;
mu.a. matrix = &mu_a_p ;
mu.b nrow = 1;
mu.b .ncol = 1;
mu.b.p = &mu_b_r ;
mu.b. matrix = femu.b.p ;
/* compute mean of a */
ret = matrix.mean (a &mu_a) ;
MATRDC.VERIFYPRECOND( ( ret = MATRDCERRORJMONE) ret ) ;
/* compute mean of b */
ret = matrix.mean (b femu.b) ;
MATRDC.VERIFYPRECOND(( ret = MATRIXPRRORNONE) ret ) ;
/* cov(a,b) = 1/M sum (( amu.a) (bmu.b) ) where M=N or N1 depending on
n.not.n.minus.one */
sum = 0.0;
for (r=0; r < a>nrow; r++) {
sum += ( a>matrix [ r ] [ 0 ] mu.a matrix [ 0 ][ 0 ] ) (b>matrix [ r ] [ 0 ] mu.b.
matrix [ 0 ] [ 0 ] ) ;
}
49


sum /= ( n_not_n_minus_one) ? a>nrow : ((a>nrow==l) ? a>nrow : a>nrow 1) ;
result = sum;
return MATRKJ5RROR-NONE;
}
/*
* matrix.cov ()
*
* Takes an NxM matrix and produces the MxM cov. matrix
* where N is the number of observations and M is the number
* of variables .
*
* matrix computes covariance of this
* result places results into this
* n.not.n.minus.one see matrix_var() for explanation
*
*/
int matrix_cov ( matrix.t *matrix matrix_t *result unsigned char
n_not_n_minus_one)
int ret = MATRDCERROFLNONE;
unsigned int r;
unsigned int rr ;
unsigned int c;
matrix.t columna;
matrix_t columnb ;
r o u I -f vrocnlt
/* asserts */
MATRK_ASSERT( matrix =
MATRK_ASSERT( result =
/* preconditions */
NULL):
NULL) ;
matrix>nrow = 0) Ml
matrix>ncol = 0) Ml
matrix>ncol = result
result > n c o 1 result

MATRIX .ERROR JNVALID_DIMENSIONS)
/* create a vector Nxl */
MATRIX_VERIFY_SET_GOTO( matrix.new (&columna matrix>nrow , 1) =
MATRDCERRORJSfONE, ret, MATRDCERROR_NO_MEMORY, free.none);
/* create another vector Nxl */
MATRK_VERIFY_SET_GOTO( matrix_new(&columnb matrix>nrow , 1) =
MATRDGERROR-NONE, ret, MATRIX_ERROR-NO-MEMORY, free.one);
/* compute diagonal */
for (c = 0; c < matrix>ncol ; C++) {
/* fill up the column vector */
for (r = 0; r < matrix>nrow; r++) {
columna matrix [ r ][ 0 ] = matrix>matrix [ r ][ c ] ;
}
/* compute variance */
matrix.var(fecolumna, fevresult n.not.n.minus.one) ;
/* fill in the diagonal */
result >matrix [c ][ c ] = vresult;
}
/* compute upper triangle and copy to lower triangle */
for (r = 0; r < result >nrow1 ; r++) { /* for each row except last one ...
*/
for (c = r + 1; c < result>ncol ; C++) { /* for each col in row r whose
index is greater than diagonal index ... */
50


/* fill up the column vectors a and b */
for (rr = 0; rr < matrix >nrow; rr++) {
columna. matrix [rr ] [0] = matrix>m at rix [ r r ] [ c ] ;
columnb matrix [ rr ][ 0 ] = matrix>matrix [ rr ][ r ] ;
}
/* compute covariance */
vector.covj&columna &columnb &vresult n.not.n.minus.one ) ;
result >matrix [ r ][ c ] = vresult; /* fill in upper triangle ... */
result >matrix [ c ][ r ] = vresult; /* and lower triangle */
}
}
mat r ix.free (fccolumnb ) ;
freefone:
mat r ix.free (fccolumna) ;
free.none :
return ret ;
}
A.3 fastica_c/cov.h
/*****************************************************************
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
#ifndef COVJH ^
#define COVJJ
int matrix.cov ( matrix.t *matrix matrix-t ^result unsigned char
n.not.n.minus.one) ;
#endif
A.4 fastica_c/cp.c
*********************
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
*****************************************************************
#include
/
^include matrix.h
/ *
* matrix_cp ()
*
* Copies contents of a into result Dimensions of a and
* result must match.
*
* a src matrix
* result dst matrix
*
*/
int matrix-cp ( matrix.t *a, matrix.t *result)
{
int r ;
51


/* asserts */
MATRDCASSERT(a .'= NULL);
MATRIXASSERT (result = NULL) ;
/* preconditions */
MATRJX_VERIFY_PRECOND( (a->nrow != 0) &&
(a> n c o1 = 0) &&
(a>nrow = result >nrow) &&
(a> n c o 1 = result > n c o 1 ) ,
MATRIX_ERRORJNVALID_DIMENSIONS) ;
for ( r =0;rnrow ; r++) {
memcpy( result >matrix [ r ] a>matrix [ r ] a>ncol sizeof ( rea 1 _t ) ) ;
}
return MATRDCEIRRORJSIONE;
}
A.5 fastica_c/cp.h
/*****************************************************************
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail. com.
*
#ifndef CP_H
#define CP-H
int matrix_cp ( matrix_t *a, matrix.t ^result);
#endif
A.6 fastica_c/det.c
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony .jackson@gmail com .
*
*****************************************************************/
#include matrix.h
#include det.h
/*
* matrix.det ()
*
* find the det of a square matrix only supported upto 4x4
*
* matrix find the det of this
* result place result here
*
*/
int matrix.det ( matrix_t *matrix real.t *result)
{
/* asserts */
MATRDCASSERT (matrix != NULL):
MATRDCASSERT (result = NULL) ;
52


/ preconditions /
MATRDCVERIFYJRECOND(( matrix->nrow != 0)
(matrix>ncol != 0) &&
(matrix>nrow <= 4) &&
( matrix>nrow = matrix>ncol) ,
MATRJX_ERRORJNVALIDJDIMENSIONS) ;
matrix>matrix [ 1 ] [ 1 ]
matrix>matrix [1] [0]
bad /
matrix>matrix [ 1 ] [ 1 ]
matrix>matrix [ 1 ] [ 2 ]
matrix>matrix [ 0 j [ 1 j
matrix>matrix [0] [2]
matrix>matrix [ 0 ] [ 1 ]
matrix>matrix [0] [2]
straight from matlab
matrix>matrix [ 1 ] [ 1 ]
switch ( matrix>nrow) {
case 1:
result = matrix>matrix [ 0 ][ 0 ] ;
break;
case 2:
result = / adbe /
matrix>matrix [ 0 J [ 0 ]
matrix>matrix [ 0 ][ 1 ]
break;
case 3:
result = / big but not
matrix>matrix [ 0 ][ 0 ]
matrix>matrix [ 0 ][ 0 ]
matrix>matrix [ 1 ] [0 ]
matrix>matrix [ 1 ][0]
matrix>matrix [ 2 ][ 0 ]
matrix>matrix [ 2 ][ 0 ]
break ;
case 4:
result = / yikes!, but
matrix>matrix [ 0 ][ 0 ]
matrix>matrix [ 3 ] [ 3 ]
matrix>mat rix [ 0 ][ 0 ] matrix>matrix [ 1 ][ 1 ]
matrix>matrix [ 3 ] [ 2 ]
matrix>matrix [ 0 ][ 0 ] matrix>matrix [ 2 ][ 1 ]
matrix>matrix [ 3 ] [ 3 ] +
matrix>matrix [ 0 ][ 0 ] matrix>matrix [ 2 ][ 1 ]
matrix>matrix [ 3 ] [ 2 ] -f
matrix>matrix [ 0 ] [ 0 ] matrix>matrix [ 3 ][ 1 ]
matrix>matrix [ 2 ] [ 3 ]
matrix>matrix [ 0 ][ 0 ] matrix>matrix [ 3 ][ 1 ]
matrix>matrix [ 2 ] [ 2 ]
matrix>matrix [ 1 ] [0 ] matrix>matrix [ 0 j [ 1 ]
matrix>matrix [ 3 ] [ 3 ] +
matrix>matrix [ 1 ][ 0 ] matrix>matrix [ 0 ][ 1 ]
matrix>matrix [ 3 ] [ 2 ] +
matrix>matrix [ 1 ][ 0 ] matrix>matrix [ 2 ][ 1 ]
matrix>matrix [ 3 ] [ 3 ]
matrix>matrix [ 1 ][ 0 ] matrix>matrix [ 2 ][ 1J
matrix>matrix [ 3 ] [ 2 ]
matrix>matrix [ 1 ][ 0 ] matrix>matrix [ 3 ][ 1 ]
matrix>matrix [ 2 ] [ 3 ] +
matrix>matrix [ 1 ][ 0 ] matrix>matrix [ 3 ][ 1 ]
matrix>matrix [ 2 ] [ 2 ] +
matrix>matrix [ 2 ][ 0 ] matrix>matrix [ 0 ][ 1 ]
matrix>matrix [ 3 ] [ 3 ]
matrix>matrix [ 2 ][ 0 ] matrix>matrix [ 0 ][ 1 ]
matrix>matrix [ 3 ] [ 2 ]
matrix>matrix [ 2 ][ 0 ] matrix>matrix [ 1 ][ 1 ]
matrix>matrix [ 3 ] [ 3 ] +
matrix>matrix [ 2 ][ 0 ] matrix>matrix [ 1 ][ 1 ]
matrix>matrix [ 3 ] [ 2 J +
matrix>matrix [ 2 ][ 0 ] matrix>matrix [ 3 ][ 1 ]
matrix>matrix [ 1 ][ 3 ]
matrix>matrix [ 2 ][ 0 ] matrix>matrix [ 3 ][ 1 ]
matrix>matrix [ 1 ][ 2 ]
matrix>matrix [ 2 ] [ 2 ]
matrix>matrix [ 2 ][ 1 ]
matrix>matrix [ 2 j [ 2 j +
matrix>matrix [ 2 ][ 1 ] +
matrix>matrix [ 1 ][ 2 ]
matrix >matrix {1 ] [ 1 ] ;
*/
matrix>matrix [ 2 ] [ 2 ]
matrix >matrix [ 2 ] [ 3 ]
matrix>matrix [ 1 ][ 2 ]
matrix>matrix [ 1 ][ 3 ]
matrix>matrix [ 1 ][ 2 ]
matrix>matrix [ 1 ][ 3 ]
matrix>matrix [ 2 ] [ 2 ]
matrix>matrix [ 2 ] [ 3 ]
matrix>matrix [ 0 ][ 2 ]
matrix>matrix [ 0 j [ 3 ]
matrix>matrix [ 0 ][ 2 ]
matrix>matrix [ 0 ][ 3 ]
matrix>matrix [ 1 ][ 2 ]
matrix>matrix [ 1 ][ 3 ]
matrix>matrix [ 0 ][ 2 ]
matrix>matrix [ 0 ][ 3 ]
matrix>matrix [ 0 ][ 2 ]
matrix>matrix [ 0 ][ 3 ]
53


matrix>matrix [ 3] [ 0] matrix>matrix [ 0 ][ 1 ]
matrix>matrix [ 2 ] [ 3 ] +
matrix>matrix [ 3 ][ 0 ] * matrix>matrix [ 0 ][ 1 ]
matrix>matrix [ 2 ] [ 2 ] +
matrix>matrix [ 3 ][ 0 ] * matrix>matrix [ 1 ][ 1 ]
matrix>matrix [ 2 ] [ 3 ]
matrix>matrix [ 3 ][ 0 ] * matrix>matrix [ 1 ][ 1 ]
matrix>matrix [ 2 ] [ 2 ]
matrix>matrix [ 3 ][ 0 ] * matrix>matrix [ 2 ][ 1 ]
matrix>matrix [ 1 ][ 3 ] +
matrix>matrix [ 3 ][ 0 ] * matrix>matrix [ 2 ][ 1 ]
matrix>matrix [ 1 ] [ 2 ] ;
break ;
default :
/* should never get here */
return MATRIX_ERRORJNVALID_DIMENSIONS;
* matrix>matrix [ 1 ][2]
* matrix>matrix [ 1 ][ 3 ]
* matrix>matrix [ 0 ][ 2 ]
* matrix >matrix [ 0 ][ 3 ]
* matrix >matrix [ 0 ][ 2 ]
* matrix>matrix [ 0 ][ 3 ]
*
*
*
*
*
*
return MATKDCERRORJSfONE;
}
A.7 fastica_c/det.h
/*****************************************************************
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony .jackson@gmail com .
*
#ifndef DET_H ^
#define DET_H
int matrix.det ( matrix.t *matrix real_t ^result);
#endif
A.8 fastica_c/eig.c
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail. com.
*
#include
#i n c1u d e
#include
^include matrix.h
#include eig.h
^include jacobi.h
/*
* matrix.eig ()
*
* Compute the eigenvalues and eigenvectors of a asymmetric* matrix.
* Relies heavily of NR in Cs jacobi() function.
*
* mat matrix from which to generate eigenvalues and eigenvectors
* dmat vector of eigenvalues
54


* vmat matrix of eigenvectors
*/
int matrix.eig (matrix.t *mat, matrix.t *dmat, matrix.t *vmat)
{
float **a;
float v ;
float *d;
int r,c ;
int nrot ;
/* asserts */
MATRDC_ASSERT( mat != NULL);
MATRDLASSERT(dmat != NULL) ;
MATRJDCASSERT(vmat != NULL);
/* preconditions */
MATRIX.VERIFY_PRECOND((mat->nrow != 0) &&
(mat> n c o 1 = 0) &&
(dmat>nrow != 0)
(dmat>ncol != 0)
(vmat>nrow != 0) &&
(vmat>ncol != 0) &&
(mat>nrow = mat>ncol) &&
(dmat>nrow = mat>nrow)
(dmat>ncol = 1) &fe
(vmat>nrow = mat>nrow) &&
(vmat>ncol = mat>ncol) ,
MATRIX_ERRORJNVALID_DIMENSIONS) ;
a = matrix (1, mat>nrow + l, 1, mat>ncol+ 1) ;
v = matrix(l, mat>nrow + l, 1, mat>ncol+ 1) ;
d = vector (1, mat>nrow+l) ;
/* fill in a with mat */
for ( r =0;rnrow ; r++) {
for ( c =0;cncol ; C++) {
/* NR in C uses 1 based indices for some dumb reason */
a[r + l][c + l] = mat>matrix [ r ] [ c ] ;
}
}
jacobi(a, mat>nrow d, v, &nrot) ;
/* fill in dmat and vmat ... */
for ( r=0;rnrow ; r++) {
dmat>matrix [ r ] [ 0 ] = djr + 1]; /* NR in C uses 1 based */
}
for ( r =0:rnrow ; r++) {
for (c = 0;cncol ; C++) {
vmat>matrix [ r ] [ c ] = v[r + l][c + l];
}
}
free_vector(d, 1, mat>nrow +1) ;
free.matrix (v 1, mat>nrow + l, 1, mat>ncol+ 1) ;
free_matrix (a 1, mat>nrow + l, 1, mat>ncol+1) ;
return MATRIX.KRROR_NONH;
}
A.9 fastica.c/eig.h
/*****************************************************************
* Copyright (c) 2008 Anthony R. Jackson
55


*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony .jackson@gmail .com.
*
#ifndef EIG_H
#define EIG_H
int matrix.eig (matrix.t *matrix, matrix_t *d, matrix.t *v) ;
#endif
A.10 fastica.c/fastica-c.c
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
#include
#include
#include
#include
#include
#include
< u n i s t d.h>


< f c n 11.h>
#include
#i n c1u d e
^include
#include
^include
#include
#include
#include
#include
#include
#include
#include
f as t i c a.h
matrix h
cp.h
i c a h
mean. h
pea.h
whiten h
mult h
#define SOURCESJSOJMBER 9
#define MDCTUFUiNUNfBER 4
#define NFRAMES 10000
struct source.s {
SNDFILE *sfp :
SFJNFO sfi;
matrix_t mat;
};
struct source.s sources [SOURCES.NUMBER] ;
struct source.s mixture [X1IXT1JREJNTJMBER] ;
matrix.t demix ;
matrix.t mix;
static jack.client.t ^client ;
static jack.port.t *output.port[2];
static jack.nframes.t sample.rate :
56


static int (* user _cb )(j ack .defau It _aud io .sample, t *out[2] jack_nframes_t
nframes void ^private) ;
/*
* This function leaks memory on error !
*/
int open.source(struct source.s *source, int id)
{
int i c ;
int ret ;
unsigned int frames.left;
unsigned int frames.read ;
unsigned int total.frames.read = 0;
double buffer ;
char fname [18 0];
if ((id < 0) || (id >= SOURCESJWMBER))
return 1;
snpr i nt f ( fname , 179, ../clips / source%d wav id +1) ;
source>sfi .format = 0;
source>sfp = sf.open (fname SFM-READ, ^source >s f i ) ;
if (! source>sfp )
return 1;
buffer = (double *) malloc ( sizeof ( double ) source >s fi frames ) ;
if ( buffer = NULL)
return 1;
for (frames.left = source >sfi .frames; frames.left != 0;) {
frames.read = sf.read.double ( source>sfp buffer, frames.left);
if (frames.read = 0) /* eof */
break ;
if (frames.read < 0) /* error */
return 1;
frames.left -=frames_read ;
total_frames_read+=frames_read ;
}
ret = matrix.new (^source >mat, source>sfi .channels total.frames.read)
i f ( ret != 0)
return 1;
for ( i=0;i s f i .channels) {
for ( c =0;c sfi channels ; C++) { /* c for channel */
source >mat. matrix [ c ][ i+c ] = buffer [ i+c ] ;
}
}
free ( buffer ) ;
sf.close(source > s f p);
return 0;
}
int mix.sources.OO ( void )
{
int i id ;
int mixture-ids [MIXTUR1LNUMBER] = { 1, 4, 6, 8 };
int nframes ;
nframes = sources [mixture, ids [0]]. sfi. frames;
for ( id =0;id mixture [id ].sfp = sources [ mixture-ids [ id ]].sfp;
mixture [id], sfi = sources [mixture.ids [id]], sfi;
57


matrix_new(&:mixture [ id ] mat, sources [mixture.ids [id ]] mat. nrow sources [
mixture_ids[id]]. mat. n c o 1) ;
}
matrix_new(&mix MKTUKEJ'JUMBER, nframes);
#define SELEM(id, i) sources
printf(mix = [\n);
printf( 0.25 0.25, 0.25,
printf(0.10, 0.30, 0.10,
printf(0.10, 0.20, 0.30,
printf(0.15, 0.10, 0.10,
printf ( ] \ n ) ;
[mixture.ids [id ] ] mat. matrix [0] [ i ]
0.1 5 ; \ n ) ;
0.10; \ n ) ;
0.2 5 ; \ n ) ;
0.10 \n ) ;
for ( i =0; i mixture [ 0 ]. mat. matrix [ 0 ][ i ] =
0.25 SELEM(2 i ) + 0.15
mixt ure [ 1 ]. mat. matrix [ 0 ][ i ] =
0.10 SELEM(2 i ) + 0.10
mixt ure [ 2 ]. mat. matrix [ 0 ][ i ] =
0.30 SELEM(2 i ) + 0.25
mixture [ 3 ]. mat. matrix [ 0 ][ i ] =
0.10 SELEM(2 i ) + 0.10
0.25 SELEM(0 i)
* SELEM (3 i);
0.10 SELEM (0 i )
* SELEM (3 i) ;
0.10 SELEM(0 i )
* SELEM (3 i);
0.15 SELEM(0 i)
* SELEM (3 i):
+ 0.25 * SELEM (1 , i) +
+ 0.30 * SELEM (1, i) +
+ 0.20 * SELEM (1 , i) +
+ 0.10 * SELEM (1 , i) +
mix matrix [ 0 ] [ i ]
mix matrix [ 1 ] [ i ]
mix matrix [ 2 ] [ i ]
mix matrix [ 3 ] [ i ]
mixture [ 0 ] . mat. m at r i x [ 0 ] [ i ] ;
mixture [ 1 j . mat. matrix [ 0 ] j i ] ;
mixture [ 2 ] . mat. matrix [ 0 ][ i ] ;
mixture [ 3 ] . mat .matrix[0][i];
return 0;
}
void error.cb(const char *desc)
{
}
void shutdown.cb ( void *arg)
{
}
int samplerate.cb (jack_nframes_t nframes, void *arg)
{
sample_rate = nframes;
return 0;
}
int process_cb (j ack _n fr ames _t nframes, void *arg)
{
jack_default-audio_sample.t *out [2] ;
out[0] = (jack.default-audio_sample_t *) j ack .port-get-buffer ( output _port [ 0 ] ,
nframes ) ;
out[l] = (j ack .default .audio.sample.t *) j ac k _p o r t _ge t _b u f fe r ( out put .port [ 1 ] ,
nframes) ;
i f ( user.cb )
return user_cb(out, nframes, arg);
return 0;
}
int j ac k i n i t ( const char *name, int (* ucb) (j ack _de f au 11 .aud i o _s amp le t *out[2] ,
jack.nframes.t nframes, void *priv), void ^private)
{
user.cb = ucb;
58


jack.set.error.function(error.cb) ;
client = j ack _c 1 i e n t n e w (name) ;
if (client = NULL)
return 1;
jack.set.process_cal lback ( c 1 ient process.cb private);
j ack_set _samp 1 e _rat e _cal 1 back ( c 1 ient samplerate.cb 0) ;
jack.on.shutdown ( client shutdown.cb 0);
output.port [0] = j ac k _po rt _r eg is t e r ( client out put .1 ,
JACK_DEFAULTJ^UDIO_TYPE, JackPortlsOutput 0) ;
output.port [ 1 ] = j ack _po rt .regist e r ( client output _r ,
JACK_DEFAULT^UDIO_TYPE, JackPortlsOutput 0) ;
return 0;
int us er _p r oce ss (j ac k _d efau 11 _au d i o _s a m p le _t *out[2] jack.nframes.t nframes
void private )
{
static int frames;
static int id =0;
static int phase = 0;
int max.id =4; /* 8 for sources 4 for mixture */
jack.nframes.t i ;
jack.default.audio.sample.t sample ;
for ( i =0;i if ((frames+i) >= sources >mat. ncol) {
i d Hh;
printf(now playing... id: %d\n id % 4);
frames = 0:
}
if (id >= max.id) {
id = 0;
phase++;
if (phase = 1) {
printf(now playing demixed sources\n);
}
if (phase = 2) {
printf(done.\n);
kill (0 SIGINT); /* kill thyself */
}
}
switch (phase) {
case 0:
sample = mix matrix [ id ][ frames+i ] ;
break ;
case 1:
sample = demix. matrix [ id ][ frames+i ] 0.25; /* scale it down, too loud
*/
break;
default :
sample = 0.0; /* silence */
break;
}
out [0 ][ i ] = sample ; /* left */
out j 1 j [ i ] = sample : /* right */
}
frames+=nframes ;
59


return 0;
}
void j ack .st art ( void )
{
int i ;
const char **ports;
if (jack .activate(c1ient) != 0)
return ;
ports = jack.get.ports ( client NULL, NULL, JackPortlsPhysical | JackPortlsInput
);
for ( i =0; port s [ i ]! =NULL; i++) {
if (( i % 2) = 0)
jack_connect(client, jack_port.name(output.port[0]), ports[i]);
else
jack_connect(client, jack_port_name(output_port[l]), ports[i]);
}
free ( ports ) ;
}
void jack.done(void)
{
if (client)
j ac k _c 1 i e n t _c 1 o s e (client);
}
void s i gi n t _h a n d le r ( i n t signum)
{
p r i n t f ( killing JACK\ n ) ;
jack.done();
matrix.pool.destroy () ;
exit (0) ;
}
int setup_memory_pool ( int nframes)
{
/ *
50000x1: 2 (max at once: 2) depth: 0
50000x4: 11 (max at once: 1) depth: 0
4x1: 13 (max at once: 3) depth: 2
4x4: 142 (max at once: 17) depth: 5
4x50000: 15 (max at once: 6) depth: 5
1x50000: 13 (max at once: 13) depth: 13
*/
matrix.pool.t pool[] = {
{ nframes , 1} ,
{ nframes 4 } ,
{4, 1},
{4, 4},
{4 nframes } ,
{ 1 nframes } ,
};
matrix-pool_init(pool sizeof(pool)/sizeof(pool [0]) ) ;
return 0;
}
int main(int argc char *argv[])
{
int i ;
60


int piecewise = 0;
i f (argc > 1) {
piecewise = 1;
printf( piecewise is TRUE\n ) ;
}
setup.memory.pool (NFRAMES) ;
for ( i =0;i open.source (&sources[ i ] i ) :
}
mix_sources_00 () ;
matrix.new (&demix mix nrow mix.ncol);
if (piecewise) { /* only process nframes at a time */
int r c ;
int nframes = NFRAMES;
matrix.t submix;
matrix.t subdemix;
matrix_new(&submix mix nrow nframes);
matrix.new (&subdemix mix nrow nframes);
for ( i =0; i pr i n t f (% f\n , ( float ) i / ( float ) mix ncol 100.0);
for ( r =0; r for ( c = 0;c submix matrix [ r ][ c ] = mix matrix [ r ][ i+c ] :
}
}
fas t i c a.d ec om pose .sou rces (&submix fesubdemix, 1000, NULL, 4);
for ( r =0;r for ( c = 0;c demix matrix [ r ] [ i+c ] = subdemix matrix [ r ] [ c ] ;
}
}
}
mat r ix .free (&submix) ;
mat r ix_free (&subdemix ) ;
} else {
fast ica.decompose.sour ces (&mix &demix , 1000, NULL, 4);
matrix.trace.dump () ;
jack.init ( fastica c user.process NULL) ;
signal (SIGINT sigint .handler ) ;
jack.start():
w h i 1 e (1) {
sleep (60) ;
}
return 0;
}
A. 11 fastica_c/fastica.h
/*****************************************************************
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
61


* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
*****************************************************************/
#ifndef FASTICA_H
#define FASTICA_H
#include matrix.h
#include ica.h
typedef struct fastica.s fastica.t;
int fastica.decompose.sources ( matrix.t *mix, matrix.t *demix ,
int maxiter int *actualiter ,
int numsources) ;
#endif /* FASTICAJI */
A.12 fastica_c/fastica-music.c
/
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
#i n c 1 u d e
#include
^include
#include
#include
#include
#include
#include
#include
#include
#i nclude
#include
#i nclude
^include
#include
#i nc1u de
#i nclude
"matrix h
cp h
i c a h
mean. h
pea.h
whiten h
mult h
#define SOURCES-NUMBER 13
^define MIXTURE-NUMBER 2
struct source.s {
SNDFILE *sfp ;
SF-INFO sfi ;
matrix_t mat:
};
struct source.s sources [SOURCESJSTJMBER] ;
struct source.s mixture [MIXTUREJCUMBER] ;
matrix.t demix ;
matrix.t mix :
static jack.client.t ^client;
62


static jack.port.t output.port [2];
static jack.nframes.t sample.rate ;
static int (* user.cb ) (jack_default.audio_sample_t *out[2], j ack.nframes.t
nframes void ^private);
/ *
* This function is leaks memory on error!
. */
int open.source(struct source.s ^source int id)
{
int i ,c ;
int ret ;
unsigned int frames.left;
unsigned int frames.read ;
unsigned int totaLframes.read = 0;
float buffer ;
char fname[180];
if ((id < 0) || (id >= SOURCES-NUMBER) )
return 1;
s npr i nt f ( fname , 179, .. / c 1 i p s / sou r ce9cd wav id + 1);
source>sfi .format = 0;
source >sfp = sf.open ( fname SFM_READ, &source >s f i ) ;
if (! source >sfp )
return 1;
buffer = (float *) malloc ( s i z eo f ( f 1 oa t ) source > s f i frames ) :
if (buffer = NULL)
return 1;
for (frames_left = source>sfi .frames; frames.left != 0;) {
frames.read = sf.read.float ( source>sfp buffer, f r ame s 1 e f t ) ;
if (frames.read = 0) /* eof */
break;
if (frames.read < 0) /* error */
return 1;
frames.left =frames.read ;
total_frames_read+=frames.read ;
}
ret = matrix.newffcsource >mat, source >s f i channels tot a 1 .frames .read )
if (ret != 0)
return 1;
for ( i =0; i s fi channels ) {
for (c=0;cs f i channels ; C++) { /* c for channel */
source >mat. matrix [ c ][ i+c ] = buffer[i+c];
}
}
free ( buffer ) ;
return 0;
}
int mix.sources.OO ( void )
{
int i id ;
int mixture-ids [MKTURE1JVUMBER] = { 11, 12 };
int nframes ;
nframes = sources[mixture_ids [ 0 ] ] sfi .frames;
for (id =0;id<^!IXTUREJSnjMBER; id++) {
mixture [ id ]. sfp = sources [ mixture.ids [ id ]]. sfp :
mixture[id].sfi = sources [mixture-ids [id]], sfi;
63


matrix_new(&;mixture [ id ] mat, sources [mixture-ids [id]], mat. nrow sources [
mixture.ids [ id ] ] mat. nco 1) ;
}
matrix_new(&;mix MIXTURE.NUMBER, nframes) ;
#define SELEM( id, i) sources [mixture-ids [id]], mat. matrix [0] [ i ]
for ( i =0; i // mixt ure [ 0 ] mat. matrix [0][i] = 0.5 SELEM( 0 i ) + 0.5 SELEM( 1 i ) ;
// mixture [ 1 ]. mat. matrix [ 0 ][ i ] = 0.10 SELEM(0, i) + 0.30 SELEM(1, i)
/* try the original unaltered sources */
mixture [ 0 ] mat. matrix [ 0 ] [ i ] = SELEM(0 i ) ;
mixture [ 1 ] mat. matrix [ 0 ][ i ] = SELEM( 1 i ) ;
mix .matrix[0][i] = mixture [0]. mat. matrix [ 0 ] [ i ] ;
mix matrix [ 1 ][ i ] = mixture [ 1 ] mat. matrix [ 0 ] [ i ] ;
return 0;
}
int decompose.sources ( void )
{
int ret ;
matrix.t A, W;
matrix_t E, D; /* E eigenvectors D eigenvalues */
raatrix.t s; /* sources */
matrix.t wm, dwm, ws; /* wm whitening matrix dwm dewhitening matrix ws
whitened source */
matrix.t means; /* means of each row in s */
matrix.t restored.s; /* restored sources */
matrix_new(&demix mix.nrow, mix.ncol);
MATf!IX_VERIFY_PRECOND( matrix _new (&s mix nrow mix.ncol) =
MATRIX JFRROILNONE. -1);
MATRDCVERIFYJRECOND( matrix-cp(&mix &s ) =
MATRDCERRORRIONE, -1);
MATRIX_VERIFY_PRECOND( matrix.new (fcws s nrow s.ncol) =
MATRIX_ERROR_NONE, -1);
MATRIX_VERIFY_PRECOND( matrix.new (&wm, s nrow s.nrow) =
MATRIX_ERROR_NONE, -1);
MATRIX_VERIFY_PREXI!OND( matrix_new(&dwm, s nrow s nrow) =
MATRIX JERROR.NONE, -1);
MATRJX_VERIFY_PRIXIOND( matrix_new(&E, s.nrow, s.nrow) =
MATRIX _ERROR_\'ONE. -1);
MATRIX_VERIFYJ:>RECOND( matrix .new (&D, s.nrow, 1)
MATRDCERROR.NONE, -1);
MATRIX_VERIFY_PRECOND( matrix-new (fcmeans s.nrow, 1)
MATRJX_ERROR_NO\E, -1);
MATRIX_VERIFY_PRECOND( matrix.new (&A, s.nrow, s.nrow) ^
MATRJX_ERROR_NONE, -1);
MATRIX_VERIFY_PRECOND( matrix .new (&W, s.nrow, s.nrow) =
MATRK_ERROR_NONE, -1);
MATRIX_VERIFYJ:>RECOND( matrix.new(&restored_s s.nrow, s.ncol) ^
MATRIX-ERRORJXONE, -1);
/* remove mean */
ret = matrix.remove.mean.transpose(&s fcmeans) ;
printf( ret : %d\n ret);
pr in t f (( remove mean)\n) ;
//mat r ix_pr i nt (&s ) ;
printf ( \n ) ;
64


/* pea */
ret = matrix_pca(&s &E, &D, 2) ;
printf( ret : %d\n ret);
printf((eigenvectors)\n);
matrix.print (&E) ;
printf((eigenvalues)\n);
matrix.print (&D) ;
printf(\n) ;
/* whiten */
ret = matrix_whiten(&s &D, &E, &ws &wm, &dwm) ;
printf(ret: %d\n ret);
printf( (whitened sources)\n) ;
//matrix_print (&ws) ;
printf(\n) ;
printf ( (whitening matrix)\n) ;
matrix.print (fain) ;
printf(\n);
printf ( (dewhitening matrix)\n) ;
matrix.print (Mwm) ;
printf(\n) ;
/* ica */
ret = matrix_ica(&ws, &wm, &dwm, 1000, MATRIXJCA_POW3, &A, &^V, NULL);
printf( ret : %d\n ret);
printf ( (A) \n ) ;
matrix.print (&A) ;
printf(\n) ;
printf( (W) \n ) ;
mat r ix.print (MV) ;
printf(\n) ;
/* restore sources */
ret = matrix.mult (&W, &s &r est or ed _s ) ;
printf(ret: %d\n ret);
printf((restored sources)\n);
//matrix_print(&restored_s) ;
printf(\n) ;
/* restore mean */
ret = m at r ix_r est ore .me an _t r anspose (fcmeans &restored_s ) ;
printf(ret: %d\n ret);
printf((restored sources)\n);
//matrix_print (&: rest ored _s ) ;
printf(\n) ;
ret = mat r ix_cp (& r es t or ed _s &demix);
printf(ret: %d\n ret);
return ret ;
}
void error.cb(const char *desc)
{
}
void shutdown.cb ( void *arg)
{
}
int samplerate.cb (j ac k _nfr ames _t nframes void *arg)
{
sample.rate = nframes;
return 0;
65


int process.cb (jack_nframes_t nframes void *arg)
{
jack-default.audio.sample.t *out [2] ;
out[0] = (j ack .default .audio.samp le _t *) j ac k _po r t _ge t _b u ffe r ( out put .port [ 0 J ,
nframes ) ;
out[l] = (jack.default.audio.sample.t ) j ac k .port-get _bu ffer ( out put-port [ 1 ] ,
nframes ) ;
i f ( user.cb )
return user_cb(out, nframes, arg);
return 0;
}
int jack_init(const char *name int (* ucb) (j ack .defau 11 _audio _samp 1 e_t out[2],
jack.nframes.t nframes, void *priv), void *private)
{
user.cb = ucb;
jack_set_error.function(error.cb ) ;
client = j ack.client .new (name) ;
if (client = NULL)
return 1;
jack_set_process_ca 11 b ack ( c 1 ient process_cb private);
j ack.set _sam pie _r at e.callbac k ( c li en t samplerate.cb , 0);
jack.on.shutdown ( client shutdown.cb 0);
output.port[0] = jack_port.register(c1 ien t output_l,
JACK JJEFAULT-AUDIO.TYPE, JackPortlsOutput 0);
output-port [1 ] = jack_port_register(client output.r,
JACK_DEFAULT_AUDIO_TYPE, JackPortlsOutput 0) ;
return 0;
}
int user_process (j ack _defau It _aud io.sample.t *out[2] j ack. nframes _t nframes,
void private )
{
static int frames;
static int id=0;
static int phase = 0;
int max.id =2; / 8 for sources 4 for mixture */
j ac k .nframes _t i;
jack .default .audio, sampled sample ;
for ( i =0; i if ((frames+i) >= sources >mat. ncol) {
id H
printf(now playing... id: %d\n id % 2);
frames = 0;
}
if (id >= max.id) {
id = 0;
phase+-P;
if (phase = 1) {
printf(now playing demixed sources\n):
}
if (phase = 2) {
printf(done.\n);
kill (0 SIGINT) ; /* kill thyself /
66


}
}
switch (phase) {
case 0:
sample = mix matrix [ id ] [ frames+i ) ;
break ;
case 1:
sample = demix. matrix [ id ][ frames+i ] 0.25; /* scale it down, too loud!
*/
break ;
defauIt ;
sample = 0.0; /* silence */
break ;
}
out [0 ] [ i ] = sample ; /* left */
out [ 1 ] [ i] = sample ; /* right */
}
frames+=nframes ;
return 0;
void j ack.st art ( void )
{
int i ;
const char **ports;
if (jack_activate(c1 i e nt) != 0)
return ;
ports = j ack.get _ports ( client NULL, NULL, JackPortlsPhysical | JackPortlsInput
);
for ( i =0; ports [ i ]! =NULL; i++) {
if (( i % 2) = 0)
jack_connect(client, jack-port.name(output_port[0j), ports [ i ]) ;
else
jack_connect(client jack_port_name(output-port[l]) ports [i]) ;
}
free(ports);
}
void jack-done(void)
{
if (client)
jack_client_close (client);
}
int out put .wave _fi les ( void )
{
int i j ;
int nframes = demix. ncol;
SNDFILE sfp ;
SFJNFO *sfi;
char *fnames[2] = { out_000 wav out _111 wav } ;
for ( i =0; i <2; i++) {
s f i = ^sources [ i 4-11 ]. s f i ;
sfp = sf.open ( fnames [ i ] SFM_WRITE, sf i ) ;
for (j =0; j enframes ; j ++) {
demix.matrix[i ] [ j ] *= 0.25:
sf.write.float (sfp &;demix .matrix[i][j]. 1):
}
67


sf.close(sfp ) ;
}
return 0;
}
void sigint .handler ( int signum)
{
pri nt f ( ki 11 i ng JACK\n);
j ack.done() ;
exit (0) ;
int main(int argc char *argv[])
{
int i ;
for ( i = 0; i open.sourcej&sources [ i ] i ) ;
}
mix_sources_00 () ;
decompose.sources () ;
matrix_trace_dump () ;
output_wave_files () ;
printf(done writing files\n);
printf( w ai ting on jackd\n);
jack_init(fas tic ac user.process NULL) ;
signal (SIGINT sigint .handler ) ;
jack_start () ;
while (1) {
sleep (60) ;
}
return 0;
}
A. 13 fastica_c/ica.c
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
#include
#include
#inelude matrix h
#include ica.h
^include sqrt.h
#inelude trans.h
#inelude mult h
^include ident.h
^include cp.h
#include op.h
68


/ *
* compute: A = dewhiteningmatrix B;
* W = B whiteningmatrix ;
*/
static int f i 11 i c a f i 11 e r s ( matrix.t *dwm, matrix.t *wm, matrix.t *B, matrix.t
*A, matrix.t *W)
{
int ret = MATRDCERROFLNONE;
matrix.t Bt;
MATRIX_VERIFY_GOTO(( ret = matrix_new(&;Bt, B>ncol B>nrow) )
MATRIX_ERROR^NONE, out);
NlATRIX_VERIFY_GOTO(( ret matrix.trans (B, &Bt) )
MATRIX_ERROR_NONE, out .free .one ) ;
MATRIX_VERIFY_GOTO (( r e t = matrix.mult (dwm, B, A))
MATRIXJEFiROR_NONE, out .free .one ) ;
MATRIX_VERIFY_GOTO(( ret = matrix.mult (&Bt wm, W) )
MATRIX J5RKOR_NONE, out.free.one) ;
out.free.one :
mat r ix_fr ee (&Bt) :
out :
return ret ;
}
/* perform elementwise division in place (i e clobbers a) */
static int d i v _i n _p 1 ace ( matrix.t *a, real.t scalar)
{
int r, c ;
for (r=0;rnrow; r++) {
for ( c =0;c
ncol ; C++) {
a>matrix [ r ] [ c ] /= scalar;
}
}
return MATRDCERROR-NONE;
}
/* perform elementwise cubing in place (ie, clobbers a) */
static int cube.in.place ( matrix.t *a)
{
int r, c ;
for ( r=0;r
nrow ; r++) {
for ( c =0;c
ncol ; C++) {
a>matrix [ r ] [ c ] = a>matrix [ r ] [ c ] a>matrix [ r ] [ c ] a>matrix [ r ] [ c ] ;
}
}
return MATRDC_ERROR_NONE:
}
/* performs element wise scalar multiplication in place ( ie clobbers a?) */
static int scalar.mult (matrix.t *a, real.t scalar matrix.t ^result)
{
int r, c ;
for ( r=0;r
nrow ; r++) {
for ( c = 0;c
ncol ; C++) {
result >matrix [ r ][ c ] = a>matrix [ r ] [ c ] * scalar;
}
return MATRIX_ERROR_NONE;
69


}
/* Implements the P0W3 contrast function which is g(x) = x3 */
/* (B* whitesig ) = whitesig B
* compute: B = (whitesig ((whitesig B) 3)) / numsamples 3 B or
* B = (whitesig ((B whitesig) 3)) / numsamples 3 B (
should give better cache performance since multiply will be column major
order)
*/
static int contrast_function_MATRIX_ICA_POW3 ( matrix.t *B, matr ix_t whitesig )
{
/ *
* tempa = 3*B
* tempb = B'T
* tempc = tempb whitesig
* tempd = tempc T
* tempd = tempd 3
* tempe = whitesig tempd
* tempe = tempe ./ whitesig >ncol
* B = tempe tempa
*/
int ret = MATRDCEFIROFLNONE;
matrix-t tempa, tempb, tempc, tempd, tempe;
MATRDCVERIF Y-GOTO ((ret
MATRDCEFIROFLNONE,
MATRDCVERIFY.GOTO ((ret
MATRECERROR-NONE,
MATRIX ATiRIFY _GOTO ((ret
MATRECERROFLNONE,
MATRIX_VERIFY_GOTO ((ret
MATRIX_ERROR_NONE,
MATRIX_VERIFY_GOTO ((ret
MATRDCERROFLNONE,
= matrix_new(&tempa
out) ;
= matrix_new (&tempb
out.free.one ) ;
= matrix_new (fctempc
out_free_two) ;
= matrix_new (&tempd
out_free_three);
= matrix_new (fctempe
out_free_four) ;
B>nrow B>ncol))
A ' , B>nrow) )
v, whitesig >ncol) ) =
tempc ncol tempc nrow) ) =
whitesig>nrow tempd. ncol)) =
MATRDCVERIFY_GOTO ((ret
MATFUX-ERROFLNONE,
MATRIX-VERIFY-GOTO ((ret
MATRDCERROR_NONE,
MATRDC_VERIFY-GOTO ((ret
MATRDCERROFLNONE,
MATRDCVERIFY-GOTO ((ret
MATRDCERROFLNONE,
MATRDCVER.IFY_GOTO ((ret
MATRDCERROFLNONE,
MATRDCVERIFY-GOTO ((ret
MATRIX ERROR NONE,
MATRDCVERIFY_GOTO ((ret
MATRIX ERROR. NONE,
MATRDCVERIFY_GOTO ((ret
MATRDCERROFLNONE,
= scalar.mult (B, 3.0, faempa) )
out.free_five ) ;
= matrix-trans (B, &tempb))
out-free-five ) ;
= matrix-mult (&tempb whitesig &tempc))
out_free_five ) ;
= matrix-transffetempc &tempd))
out_free_five ) ;
= c ube-i n _place (&tempd) )
out-free_five ) ;
= matrix.mult ( whitesig &tempd fetempe) )
out_free_five ) ;
= d iv_i n .place (&tempe whitesig>ncol) )
out-free_five ) ;
= matrix_sub(&tempe faempa, B) )
out_free_five ) ;
out-free.five :
mat r ix.fr ee (fetempe) ;
out.free.four :
matrix-free (&tempd) ;
out _free_t hree:
matrix-free (fetempc) ;
out.free.two :
mat r ix_fr ee (&tempb) ;
out _free_one :
mat r ix_fr ee (fctempa) ;
out :
return ret ;
70


}
/ *
* compute: minAbsCos = min(abs(diag(B BOld));
*
* Useful for checking for convergence of fixed point iteration .
*
*/
static int con verged ( matrix.t *Bt matrix.t *BM. real.t epsilon int *convergd)
{
int i ;
int ret = MATRI5CERROR J'JONE;
real.t minAbsCos;
matrix.t tempa:
*convergd = 0;
MATRIX_VERIFY_GOTO(( ret = matrix.new (&tempa Bt>nrow BM->nco 1) ) =
MATRIX JiRROFLNONE, out);
MATRDCVERIFY.GOTO ((ret = matrix.mult (Bt, BM. fetempa)) = MATRDCERRORJMONE,
out.free.one ) ;
minAbsCos = fabs (tempa matrix [0][0]) ;
for ( i =0; i if ( fabs (tempa matrix [ i ][ i ] ) < minAbsCos) {
minAbsCos = fabs (tempa .matrix[i][i]);
}
}
if ((1 minAbsCos) < epsilon) {
*convergd = 1;
}
out.free.one :
mat r ix.fr ee (&;tempa) ;
out :
return ret ;
}
/* Symmetric orthogonalization done in place ( ie clobbers B) */
/* B = B real(inv(B B)'(1/2)); */
/* B = B (BT B) (-1/2) */
int matrix.symmetric.orthogonalization ( matrix.t *B)
{
/ *
* tempa = B ~ T B
* tempb = (tempa) (1/2)
* tempc = B tempb
* B = tempc
*/
matrix.t Bt;
matrix.t tempa;
matrix.t tempb;
matrix.t tempc;
int ret = MATRTCERRORJVONE ;
/* asserts */
MATRDCASSERT(B != NULL);
/* preconditions */
MATRIX_VERIFY.PRECOND ((B> nrow != 0) &&
(B>ncol != 0) &&
(B>nrow >= B>ncol) /* otherwise at least one
eigenvalue is 0 */
MATRIX_ER,RORJNVALID_DIMENSIONS) ;
71


MATRIX.VERIFY-GOTO(( ret = matrix_new (&;Bt, B>ncol B>nrow) ) =
MATRIX_ERROfLNONE, out);
MATRIX_VERIPY_GOTO ((ret matrix.new (&tempa B>ncol B>ncol))
\L-\THIX_ERROR_NONE, out.free.one ) ;
MATRDC_VERIFY_GOTO(( ret = matrix_new(&tempb B>ncol B>ncol ) )
MATRIX_ERROR_NONK. out.free.two ) ;
MATRIX_VEPUFYGOTO(( ret matrix.new(&tempc B>nrow B>ncol ) ) =
MATRIX_ERROR_NONE, ou t .free _t h ree ) ;
MATRIX_VERIFY_GOTO(( ret = matrix.trans (B, &Bt) ) =
MATRJX_ERROR_NONE, ou t _fr ee _fou r ) ;
MATRIX.VERIFY_GOTO((ret = matrix.mult (&Bt B, &tempa) )
MATRIX_ERRORJVONE, out_free_four) ;
MATRIX-VERIFY.GOTO (( r e t = matrix.sqrt (fctempa &tempb 1)) =
MATRjX_ERKOR_NONE, o u t _f r ee.fou r ) ; /* tempa(1/2) */
MATRJX_VERIFY_GOTO(( ret = matrix.mult (B, &tempb fctempc) ) =
MATRJX_ERKOR_NONE, out.free.four ) ;
MATRIX_VERIFY_GOTO (( r e t = matrix_cp(&;tempc B) ) =
MATRIX _ERROR_XONE. out.free.four ) ;
out.free.four :
mat rix .free (&tempc) ;
out.free.three :
matrix.free (&;tempb) ;
out.free.two :
matrix.free (&tempa) ;
out.free.one :
mat r ix.free (&Bt) ;
out :
return ret ;
}
/* fill a 2D array with value */
static int f i 11 m a t r i x ( i n t nrow int ncol real.t **matrix, real.t value)
{
int r c ;
for ( r =0;r for ( c =0;c matrix[r][cj = value;
}
}
return 0;
}
/* Initializes B and BM to and initial guess, currently the guess is always the
identity matrix */
static int initialize.B.BM (matrix.t *B, matrix.t *BM, matrix.t whiteningmat)
{
int ret = MATRDCERROFLNONE;
matrix.t guess ;
MATRIX_VERIFY_GOTO(( ret = matrix.new (&guess whiteningmat >nrow whiteningmat
>ncol) ) = MATRJX-ERRORJSTONE, out) ;
/* BOld = zeros ( size (B) ) ; */
fi 11 _matrix (BM->nrow BM->ncol BM->matrix , 0.0);
matrix_fi 11.ident (&guess ) ; /* guess the identity matrix for now */
/* B = whiteningmatrix guess; */
MATRIX_VERIFY_GOTO(( ret = matrix.mult ( whiteningmat &guess B) ) =
MATRIX J7.RR0R_\0.\E, out);
matrix.freef&guess ) ;
out :
return ret ;
72


}
/ *
* _matrix_ica ()
*
* Perform the actual fast ica fixed point algorithm.
*
* whitesig whitened signals (each row is a separate data set)
* whiteningmat whitening matrix
* dewhiteningmat dewhitening matrix
* maxiter maximum number of iterations to perform before giving up
* c o n t r as t _f u n c t i o n contrast function to use, only POW3 is currently
implemented
* A mixing matrix
* W demixing matrix
* iters number of iterations actually performed (may be NULL)
*
*/
static int _matrix_ica ( matrix_t *whitesig ,
matrix.t whiteningmat ,
matrix.t *dewhiteningmat ,
unsigned int maxiter ,
mat r ix.cont r as t _func t ion _t con t ras t-funct ion ,
matrix.t *A,
mat rix.t *W,
int liters)
{
int i ;
matrix.t B; /* contains ws */
matrix.t BM; /* previous B */
matrix.t Bt; /* B transpose */
int ret = MATRIJCERROFLNONE;
real.t epsilon = 0.000001;
i n t conv ;
MATRDCVERIF'Y.GOTO (( r e t = matrix_new(&B, whiteningmat>nrow whiteningmat>
ncol) ) = MATRDCERROFLNONE, out) ;
MATFUX.VEFIIF'Y.GOTO(( ret = matrix_new(&BM, B.nrow, B.ncol)) =
MATRDCERROFLNONE, out .free .one ) ;
MATFUX.VERIF Y.GOTO (( r e t = init ialize.B.B M (&B, &BM, whiteningmat)) =
MATRDCERROFLNONE, out.free.two ) ;
MATRIX.VERIFY_GOTO ((ret = matrix.new (&Bt, B.ncol, B.nrow)) =
MATRDCERROFLNONE, out.free.two);
for (i = 0; i < maxiter; i++) {
/* Symmetric orthogonalization => B = B* real(inv(B B)"(1/2)) */
MATFRX.VEFilFY.GOTO (( r e t = matrix.symmetric.orthogonalization (&B) ) =
MATRFXJ5RROFLNONE, out.free.three);
MATFIDCVEFIFFY.GOTO((ret = matrix.trans (&B, &Bt) ) = MATRDCERROFLNONE,
out.free.three) ;
/* Check for convergence */
MATRJX.VEFUFY.GOTO(( ret = converged(&Bt, &BM, epsilon faonv) ) =
MATFFDCERROFLNONE, out.free.three);
if (conv) {
/* pri ntf ( converged after: %d i t e r a t i o n s\n , i);*/
break ;
}
/* BOld = B; */
MATRIX.VERFFY.GOTO((ret = matrix.cp(&B, &BM) ) = MATRDCERROFLNONE,
out.free.three) ;
switch(contrast_function) {
case MATRDCJCA_POW3:
//B = (whitesig ((whitesig B) 3)) / numsamples 3 B
73


MATRJX_VERIFY.GOTO(( ret = contrast.function_MATRIX_ICA_POW3(&B, whitesig)
) = MATRIX-ERROFLNONE, out.free.three) :
case MATRIXJCA-TANH:
//hyptan = tanh(al whitesig * B)
//B = whitesig hyptan / numsamples ones ( size (B, 1) ,1) sum (1
hyptan 2) B / numsamples al ;
break ;
}
}
/* number of iterations it took to converge */
if (iters) {
/* 1 means 'didnt converge before max iterations */
if ( i = maxiter) {
liters = 1;
} else {
liters = i +1:
}
/* Calculate ICA filters */
MATRIX_VERIFY_GOTO(( ret = f i 11. i c a f i 11 e r s ( dewhiteningmat whiteningmat &B,
A, W)) = MATRDGERRORJSrONE, out.free.three);
out.free.three:
mat r ix_free (&Bt) ;
out.free.two :
matrix.free (&BM) ;
out.free.one :
matrix.free (&B) ;
out :
return ret ;
}
/*
* matrix.ica ()
*
* Peforms the FastICA fixed point symmetric approach (not deflation)
* algorithm. Mostly checks parameters, the real work is done in
* _mat r ix_ica () .
*
* whitesig whitened signals (each row is a separate data set)
* whiteningmat whitening matrix (matrix that whitened original signal)
* dewhiteningmat dewhitening matrix (matrix needed to dewhiten)
* maxiter maximum number of iterations to perform before giving up
* contrast-function contrast function to use, only POW3 is currently
implemented
* A mixing matrix
* W demixing matrix
* iters number of iterations actually performed (may be NULL),
or
* 1 if didnt converge before maxiter iterations
*
*/
int matrix.ica ( matrix.t *whitesig , /* inputs */
matrix.t whiteningmat ,
matrix.t *dewhiteningmat ,
unsigned int maxiter ,
mat r ix.con t r ast .fu nct ion.t cont r ast _funct ion ,
matrix.t *A, /* outputs */
matrix.t *W,
int *iters)
{
int ret = MATRIXJ^RRORJMONE;
74


/* asserts */
MATRDC_ASSERT ( w h i t e s i g != NULL);
MATRIX_ASSERT ( whiten ingmat != NULL);
MATRIX ASSERT ( dewhiteningmat != NULL);
/* preconditions */
MATRIX_VERIFY_PRECONE)(( whitesig >nrow != 0) &&
( w h i t e s i g > n c o 1 = 0)
( whiteningmat >nrow = whitesig >nrow)
( whiteningmat>nrow = whiteningmat>ncol) &&
( dewhiteningmat >nrow = whiteningmat >nrow)
( dewhiteningmat>ncol = whiteningmat >ncol)
(A>nrow = dewhiteningmat>nrow)
(A>ncol = whiteningmat >ncol) &&
(W->nrow = A>nrow) &&
(W->ncol = A>ncol ) ,
M ATRIX_ERROR JN VALID_DIMEN SIONS) ;
ret = _matrix_ica ( whitesig ,
whiteningmat ,
dewhiteningmat ,
maxiter ,
contrast_function ,
A, W,
iters ) ;
return ret ;
}
A.14 fastica.c/ica.h
/*****************************************************************
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
*****************************************************************/
#ifndef ICA.H
#define ICA_H
typedef enum {
MATRIX JCA_POW3,
MATRDUCA.TANH
} mat r ix.con t rast _funct ion _t ;
int mat r ix _sy m met r ic _ort hogon al iz at ion ( mat r i x _t *B) ;
int matrix_ica ( matrix_t *whitesig , /* inputs */
matrix_t whiteningmat ,
matrix_t ^dewhiteningmat ,
unsigned int maxiter ,
mat r ix.cont rast _fu net ion _t contrast_function ,
matrix_t *A, /* outputs */
mat r ix_t *W,
int iters);
#endif
A.15 fastica.c/ident.c
* Copyright (c) 2008 Anthony R. Jackson
*
75


* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
#include matrix.h
/*
* mat r i x _f i 11 _i d e n t ()
*
* Fill a square matrix a with the identity matrix.
*
* a square matrix to fill
*
*/
int matrix_fi 11 _ident ( matrix.t *a)
int r c ;
/ asserts */
MATRK_ASSERT(a != NULL);
/* preconditions */
NfATRIX_VERIFYJ>RECOND((a->nrow != 0) &&
(a> n c o1 = 0)
(a>nrow = a>ncol) ,
MATRIX_ERRORJNVALID_DIMENSIONS) ;
for ( r =0;r
nrow ; r++) {
for ( c =0;c
ncol ; C++) {
if (r = c) {
a>matrix [ r ] [ c ] = 1.0;
} else {
a>matrix [ r ] [ c ] = 0.0;
}
}
}
return MATRDC_ERROR_NONE;
}
A.16 fastica_c/ident.h
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
#ifndef IDENT_H ^
#define IDENTJI
int matrix_fi 11 _ident ( matrix_t *a) ;
#endif /* IDENTJI */
A.17 fastica_c/inv.c
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
76


* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
^include matrix.h
#include det.h
#inelude inv.h
/ *
* matrix_inv ()
*
* Compute the inverse of a square matrix given that it is lxl, 2x2,
* 3x3, or 4x4.
*
* matrix square matrix to invert
* result inverse of matrix placed here
*
*/
int matrix.inv (matrix.t ^matrix matrix.t ^result)
{
re a 1 _t det ;
/* asserts */
MATRDCASSERT( matrix != NULL):
MATRIX ASSERT (result = NULL) ;
/* preconditions */
MATRIX_VERIFYJ,RECOND((matrix->nrow != 0) &fc
( matrix>ncol != 0)
( matrix>nrow = matrix>ncol) ,
MATRIX_ERRORJNVALID_DIMENSIONS) ;
/* compute determinant */
matrix.det ( matrix &det) ;
if (det = 0.0) {
return MATRDCERROR-SINGULAR;
}
switch ( matrix>nrow) {
case 2:
result >matrix [ 0 ][ 0 ] = matrix>matrix [ 1 ][ 1 ]
matrix>matrix [ 0 ][ 1 ] / det;
result >matrix [ 1 ] [ 0 ] = matrix>matrix [ 1 ] [ 0 J
matrix>matrix [ 0 ][ 0 ] / det;
break ;
case 3:
#define aa matrix>matrix [ 0 ][0 ]
#define bb matrix>matrix [ 0][ 1 ]
#define cc matrix>matrix [ 0][ 2 ]
#define dd matrix >matrix [ 1 ][0 ]
#define ee matrix>matrix [ 1 ][ 1 ]
matrix>matrix [1] [2]
matrix>matrix [2] [0]
matrix>matrix [ 2 ] [ 1 ]
matrix>matrix [2] [2]
#define
#define
#define
^define
ff
gg
hh
result >matrix [ 0 ] [ 0 ]
ii cc*hh) / det ;
result >matrix [ 1 ] [0 ]
i i cc *gg) / det ;
result >matrix [ 2 ] [0 ]
hhbb*gg) / det;
#undef aa
#undef bb
#undef cc
#undef dd
= ( ee* ii ff *hh) / det;
result >matrix [ 0 ][ 2 ] =
= (dd* ii ff *gg) / det;
result >matrix [ 1 ] [ 2 ] =
= (dd*hh+ee*gg) / det;
result >matrix [ 2 ] [ 2 ] =
/ det; res ult >matrix [ 0 ] [ 1 ] =
/ det; result>matrix [1 ] [1] =
result >matrix [0] [1]
(bb* ff+cc*ee ) / det;
result> mat rix [1] [1]
(aa*ff+cc*dd) / det;
result > mat rix [ 2 ] [ 1 ]
(aa*ee-bb*dd) / det;
= (bb*
= (aa*
= (aa*
77


#undef ee
#undef ff
#undef gg
#undef hh
#undef i i
break ;
case 4:
#def i ne aa
#def i ne bb
#define cc
#def i ne dd
^define ee
^define ff
#define gg
#define hh
#define i i
#define j j
^define kk
^define 11
^define mm
^define nn
#define oo
#define PP
matrix>matrix [0] [0]
mat rix >mat r ix [ 0 ] [ 1 ]
matrix>matrix [0] [2]
matrix>matrix [0] [3]
matrix>m at rix [ 1J [0 ]
matrix>m at rix [ 1 ] [ 1 ]
matrix>matrix [ 1 ] [ 2 ]
matrix>m at rix [ 1 ] [3 ]
matrix >matrix [2] [0]
matrix >m at rix [ 2 ] [ 1 ]
matrix>matrix [2] [2]
matrix>matrix [2] [3]
matrix>matrix [3] [0]
matrix >matrix [3] [1]
matrix >matrix [3] [2]
matrix>matrix [3] [3]
result >matrix [ 0 ] [ 0 ] = (ff *kk*pp+f f 11 *oo+j j *gg*pp j j *hh*oonn*gg* 11+nn*
hh*kk) / det ; resu It >matrix [ 0 ] [ 1 ] = (bb*kk*ppbb* 11 *ooj j *cc*pp+j j *dd
*oo+nn* cc 11 nn*dd*kk)/det ; resu 11 >matrix [ 0 ] [ 2 ] = (bb*gg*ppbb*hh*oo
ff *cc*pp+ff *dd*oo+nn*cc*hhnn*dd*gg)/det ; resu 11 >matrix [ 0] [3 ] = ( bb*
gg* 1 l+bb*hh*kk+f f *cc* 11 ff *dd*kkj j *cc*hh+j j *dd*gg) / det ;
result >matrix [ 1 ] [ 0 ] = ( ee *kk*pp+ee 11 oo+ i i *gg*pp i i *hh*oo-mm*gg 11 -fmm*
hh*kk) / det; result >matrix [ 1 ][ 1 ] = ( aa*kk*ppaa* 11 *oo i i cc *pp+ i i *dd
*oodmm* cc 11 mm*dd*kk)/det ; result >matrix [ 1 ][ 2 ] = (aa*gg*ppaa*hh*oo
ee*cc*pp+ee*dd*oo-fmm*cc*htnnm*dd*gg) / det ; result >matrix [ 1 ] [3 ] = (aa*
gg* 1 l+aa*hh*kk+ee*cc* 11 ee*dd*kki i *cc*hh+i i*dd*gg)/det;
result >matrix [ 2 ] [ 0 ] = (ee j j *pp+ee 11 *nn+ i i f f *pp i i *hh*nn-mm* f f 1l-f*nm*
hh*jj) / det; result >matrix [ 2 ][ 1 ] = (aa* j j *ppaa* 11 *nni i *bb*pp+i i *dd
*nn-ftnm*bb* 11 -mm*dd* j j ) / det ; result >matrix [ 2 ] [ 2 ] = (aa*ff *ppaa*hh*nn
ee *bb*pp+ee *dd*nn-ffnm*bb*hlnnm*dd* f f ) / det ; result >matrix [ 2 ] [ 3 ] = ( aa*
f f 11+aa *hh* j j +ee bb *11 ee *dd* j j i i *bb*hh+ ii*dd*ff)/det;
result >matrix [ 3 ] [ 0 ] = ( ee j j *oo+ee *kk*nn+i i f f *oo i i *gg*nn-nm* f f *kk-fnm*
gg*jj) / bet; result >matrix [ 3 ][ 1 ] = ( aa* j j *ooaa*kk*nn i i *bb*oo+i i *cc
*nnd*rim*bb*kk-mm* cc j j ) / d e t ; result >matrix [ 3 ] [ 2 ] = (aa* f f *ooaa*gg*nn
ee*bb*oo+ee cc *nn-ffnm*bb*gg^nm*cc f f ) / det ; result >matrix [3] [3] = (aa*
f f *kk+aa*gg* j j -fee *bb*kkee *cc j j bb*gg* ii + ii*cc*ff)/det;
aa
#undef
#undef
#undef
#undef
#undef
#undef
#undef
#undef
#undef
#undef
#undef
#undef
#undef nm
#undef nn
#undef oo
#undef pp
break ;
default :
return
}
bb
cc
dd
ee
ff
gg
hh
i i
j j
kk
11
MATRIX .ERRORJNVALID DIMENSIONS;
return MATRIX_ERROR_NONE;
78


}
A.18 fastica_c/inv.h
/*****************************************************************
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
*****************************************************************/
#ifndef INV.H
#define INV.H
int matrix.inv ( matrix.t ^matrix matrix.t ^result);
#endif
A.19 fastica.c/jacobi.h
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
#ifndef JACOBLH ^
#define JACOBLH
float *vector(long nl, long nh);
void free_vector(float *v, long nl, long nh);
float **matrix (long nrl long nrh long ncl long nch);
void free.matrix ( float **m, long nrl, long nrh, long ncl, long nch);
void jacobi(float **a, int n, float d[] float **v, int *nrot);
#endif
A.20 fastica.c/libfastica.c
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail.com.
*
/* This file is deprecated and will be removed */
int decompose.sources ( void )
{
int ret ;
matrix.t A, W;
matrix.t E, D; /* E eigenvectors D eigenvalues */
matrix.t s; /* sources */
matrix.t wm, dwm, ws; /* wm whitening matrix dwm dewhitening matrix ws
whitened source */
matrix.t means; /* means of each row in s */
79


matrix.t restored.s; /* restored sources */
matrix_new(&demix mix nrow mix.ncol);
MATRIX_VERIFYPRECOND( matrix.new (&s mix nrow mix.ncol)
MATRIX ERROR. NOXH. -1);
MATRIX_VERIFY_PRECOND( matrix-cp (&mix &s ) ==
MATRIX_ERROR_NONE, -1);
MATRIX_VERIFY_PRECOND( matrix .new (&ws s nrow s ncol ) =
MATRDCERROR_NONE, -1);
MATRDC_VERIFYJ>RECOND( matrix.newffewm, s nrow s.nrow) =
MATMX_ERROFLNONE, -1) ;
MATRIX_VERIFY_PRECOND( matrix_new(&dwm, s.nrow, s.nrow) =
MATRIXJERROR-NONE, -1) ;
MATRIX.VERIFY_PRECOND( matrix.new (&E, s.nrow, s.nrow) =
MATKDCERRORJ'JONE, -1);
MATRJX_VERIFYJ:,RECOND(matrix_new(&D, s.nrow, 1) =
MATRIX ERROR. NONE, -1);
MATRIX_VERIFY_PRECOND( matrix.new (femeans s.nrow, 1) =
MATRIX .F.RRORXOXE. -1);
MATRIX_VERIFY_PRECOND( matrix.new (&A, s.nrow, s.nrow) =
MATRDCERRORXIONE, -1);
MATRIX_VERIFY.PRECONr)( matrix.new (&W, s.nrow, s.nrow) =
MATRIX _ERROR_NONE, -1);
MATRIX_VERIFY_PRECOND( matrix. new(&restored_s s.nrow, s.ncol) =
MATRIX_ERROR_NONE, -1);
/* remove mean */
ret = matrix_remove.mean_transpose(&s &means);
pr i n t f ( ret : %d\n ret);
pr i n t f (( remove mean)\n) ;
/ / matrix.print (&s) ;
printf(\n) ;
/* pea */
ret = matrix.pca(&s &E, &D, 4);
printf(ret: %d\n ret);
printf((eigenvectors)\n);
matrix.print (&E) ;
printf((eigenvalues)\n);
matrix.print (&D) ;
printf(\n) ;
/* whiten */
ret = matrix.whiten(&s &D, &E, &ws, &wm, &dwm) ;
printf(ret: %d\n ret);
printf ( (whitened sources)\n) ;
//matrix.print (&ws) ;
printf)\n);
printf) (whitening matrix)\n) ;
matrix.print (tan) ;
printf)\n) ;
printf ( (dewhitening matrix)\n) ;
matrix.print (&dwm) ;
printf)\n) ;
ret = matrix_ica (fcws &wm, &dwm, 1000, MATRIXJCA_POW3, &A, &W) ;
printf) ret : %i\n ret);
printf ( (A) \n) ;
matrix.print (&A) ;
printf)\n) ;
printf( (W) \n ) ;
matrix.print (&W) ;
80


printf(\n) ;
/* restore sources */
ret = matrix.mult (&W, printf( ret : %d\n ret);
printf((restored sourees)\n);
/ / matrix.print(&restored_s ) ;
printf(\n);
/ restore mean */
ret = matrix_restore.mean _transpose (&means &rest ored _s ) ;
pr i n t f(ret : %d\n ret);
printf((restored sources)\n);
/ / matrix_print(&restored.s) ;
printf(\n) ;
ret = matrix.cp(&restored.s &demix);
printf) ret : %d\n ret);
return ret ;
A.21 fastica.c/Makefile
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
# Copyright (c) 2008 Anthony R. Jackson
#
# This software is proprietary but listed in A Real-time
# Implementation of FastICA for Blind Source Separation of Audio
# Signals for evaluation purposes and for the fulfillment of its
# requirements. To license this software contact the author at
# anthony.jackson@gmail.com.
#
////////////////////////////////////M////////////////////////////////////////////////////////////////////////////////////////////
CFLAGS=Wall -pg
ifeq ($(PROD) 1)
CFLAGS+=02
else
CFLAGSF=-g DVERIFYPRECOND
endif
LDFLAGS=1 m -lsndfile
OBJS=mean. o \
var o \
cov.o \
det o \
inv o \
eig o \
jacobi.o \
trans.o \
op.o \
mult o \
pea.o \
whiten o \
cp.o \
i c a o \
sqrt.o \
ident.o \
matrix o \
api o
all: matrix.utest fasticac fast icamusic libfastica.a
libfastica.a: $(0BJS)
ar ruv S@ $?
matrix_utest : $(OBJS)
81


fasticac: S(OBJS) fastica c.c
$(CC) $(CFLAGS) pkgconfig cflags libs jack gtkH2.0 $(LDFLAGS)
fasticac.c o $@ $(OBJS)
fastica music : S(OBJS) fasticamusic c
$(CC) S(CFLAGS) pkgconfig cflags libs jack gtk-j2.0 $(LDFLAGS)
fast icamusic c o $@ $(OBJS)
.PHONY: clean all
clean :
rm f o
rm f ~
rm f matrix_utest
rm f fastica c
rm f fast icamusic
rm f libfastica.a
A.22 fastica_c/matrix.c
* Copyright (c) 2008 Anthony R. Jackson
*
* This software is proprietary but listed in 'A Real-time
* Implementation of FastICA for Blind Source Separation of Audio
* Signals for evaluation purposes and for the fulfillment of its
* requirements. To license this software contact the author at
* anthony.jackson@gmail. com.
*
*****************************************************************/
#include
#include
#include
^include real.h
^include matrix.h
typedef struct pool.elem.s {
unsigned int nrow;
unsigned int ncol ;
unsigned int count ;
unsigned char *free;
matrix.t ^matrix;
} pool.elem.t ;
static pool.elem.t *matpool = NULL;
static unsigned int matpool.count ;
/*
* matrix.print ()
*
* Print to stdout, in a MATLAB friendly way, the contents of matrix.
*
* matrix matrix to dump
*
*/
int matrix.print ( matrix.t ^matrix)
{
unsigned int r ;
unsigned int c;
MATRIX-ASSERT( matrix != NULL);
MATRJX_VERIFYJPRECOND((matrix->nrow != 0) &&
(matrix>ncol != 0) ,
MATRIX_ERRORJNVALID.DIMENSIONS) ;
printf ( [ ) ;
82


for ( r=0;rnrow ; r++) {
for (c =0;cncol ; C++) {
print f ( %g matrix>matrix [ r ] [ c ] ) ;
}
if (r != ( matrix>nrow 1)) {
printf(:\n) ;
}
}
p r i n t f ( ] \ n ) ;
return MATRJX.ERROFLNONE;
/ *
* _matrix_new ()
*
* The actual matrix allocator. An array of pointers is allocated
* first for each row, then space for each row is allocated.
*
*/
static int _matrix_new (matrix.t *matrix unsigned int nrow unsigned int ncol)
{
unsigned int r ;
r e a 1 -1 * mat;
MATRIX ASSERT (matrix != NULL);
MATRIX_VERIFY_PRECOND (( nrow != 0)
(ncol != 0).
MATRIX_ERRORJNVALID_DIMENSIONS) ;
mat = (real_t *) malloc ( nrow s i z e o f ( r e a 1 t *));
if (mat = NULL) {
return MATRDt-ERRORJ'JO.MEMORY;
}
for (r = 0; r < nrow ; rT+) {
mat [ r ] = malloc(ncol s i z eo f ( r e al _t ) ) ;
if (mat = NULL) {
return MATRIX ERRORJNOMEMORY; /* memory leak! */
}
}
matrix>nrow = nrow ;
matrix>ncol = ncol ;
matrix>matrix mat;
return MATRIX_ERROR_NONE;
}
/*
* _matrix_free ()
*
* Actually deallocate a matrix. Dual of _matrix_new () .
*
* matrix matrix to deallocate.
*
*/
static int .mat rix .free ( mat rix_t ^matrix)
{
unsigned int r ;
MATRIX JVSSERT( matrix != NULL);
for (r = 0; r < matrix>nrow; r++) {
free ( matrix> mat rix [ r ] ) ;
matrix>matrix [ r ] = NULL;
}
free ( matrix>m at rix ) ;
83


matrix>matrix = NULL;
matrix>ncol = 0;
matrix>nrow = 0;
return MATRIX_ERR,OR_NONE;
}
typedef struct matr ix .1 ist_s {
int maxdepth ;
int depth;
int count;
int nrow ;
int ncol;
struct m at r i x _1 i s t _s *next;
} matrix_list_t ;
static m at r i x _1 i s t _t *head = NULL;
/* Finds a matrix of the correct size on the pool and returns it or
error if one cant be found. */
static int matrix_pull_from_pool ( matrix.t *matrix unsigned int nrow
int ncol)
{
int 1 j ;
for ( i =0; i if (( matpool [ i ] nrow = nrow) &&
( matpool [ ij ncol = ncol)) {
for (j =0; j if ( matpool [ i ]. free [ j ] ) {
matpool [i ] free [j ] = 0;
matrix>nrow = matpool [ i ] nrow ;
matrix>ncol = matpool [ i ] ncol ;
matrix>matrix = matpool [ i ] matrix [ j ] matrix ;
return MATf!lX_ERR,OR_NONE:
}
}
}
}
return MATRIX_ERRORJNVALID_DIMENSIONS;
}
/* Puts a pooled matrix back into the pool. */
static int matrix_return.to.pool ( matrix.t ^matrix)
{
int 1 j ;
for ( i =0; i if (( matpool [ i ]. nrow = matrix>nrow) &&
( matpool [ i ]. ncol = matrix>ncol ) ) {
for (j =0; j if ( matpool [ i ]. matrix [ j ]. matrix = matrix>matrix ) {
MATRDCASSERT( matpool [ i ]. free [j ] = 0);
matpool [ i ] free [ j ] = OxFF ;
matrix>nrow = 0;
matrix>ncol = 0;
matrix>matrix = NULL;
return MATR.IX_ERROR_NONE;
}
}
}
}
return MATRIXJ1RRORJNVALIDDIMENSIONS;
}
returns an
unsigned
84


/*
* matrix.new ()
*
* Allocates a matrix matrix from the matrix memory pool, if
* intialized otherwise from the heap. If the pool is initialized
* and the requested size is not available a warning is printed .
*
* matrix matrix structure to populate once a matrix is allocated
* nrow number of desired rows
* ncol number of desired columns
*
. */
int matrix.new ( matrix.t *matrix, unsigned int nrow, unsigned int ncol)
{
int found = 0;
m a t r i x _1 i s t _t *iter;
mat r ix _1 is t _t *nitem;
/* keep statistics on requests */
for (iter = head; iter; iter=iter>next) {
if ((iter>nrow = nrow) &&
(iter > ncol = ncol)) {
iter >count++;
iter >depth++;
if (iter>depth > iter >maxdepth) {
iter >maxdepth = iter>depth;
}
found = l;
break;
}
}
if (! found ) {
nitem = (matrix_list_t *)malloc(sizeof(matrix_list_t));
nitem>nrow = nrow ;
nitem>ncol = ncol ;
nitem>count = 1;
nitem>maxdepth = 1;
nitem>depth = 1;
nitem>next = head ;
head = nitem ;
}
/* check pool first */
if ( matrix_pull_from_pool (matrix nrow, ncol) = MATREX_ERROR_NONE) {
return MATRDCERROFLNONE;
}
if (matpool != NULL) {
print f (WARNING slow malloc increase pool size for %ux%u matrices.\n
nrow ncol) :
}
/* make new matrix on the heap */
return _matrix_new ( matrix nrow, ncol);
}
/ *
* matrix_free()
*
* Deallocate a matrix. If it can be returned back to the pool this is
* done, otherwise it is free() d from the heap.
*
* matrix matrix to return to pool or free() from heap
*
*i
int matrix.free ( matrix.t ^matrix)
85


{
int found = 0;
m a t r i x 1 i s t _t *iter;
MATRIX ASSERT (matrix != NULL);
/* adjust statistics */
for (iter = head; iter; iter=iter>next) {
if ((iter>nrow = matrix>nrow)
(iter>ncol = matrix>ncol) ) {
iter >depth ;
found = l;
break ;
}
}
if (! found ) {
return 1;
}
/* try to return to pool first */
if ( mat rix.ret urn_to_pool ( matrix ) = MATRDC-ERROR.NONE) {
return MATRK-ERRORNONE;
}
/* back to the heap */
return .matrix.free (matrix) ;
}
/*
* matrix.trace.dump))
*
* Dump the gathered statistics on matrix.new() requests.
*
*/
int matrix.trace.dump ( void )
{
m at r i x _1 i s t t iter;
for (iter = head; iter; iter = iter>next) {
printf (%dx%d: %d (max at once: %d) depth; %d\n, iter>nrow, iter>ncol,
iter >count iter >maxdepth iter >depth ) ;
}
printf(\n) ;
return 0;
}
/*
* mat rix.new-from-buffer ()
*
* Create a matrix from user supplied memory. Useful for avoiding
* memcpy() s and double nested loops to populate matrices from 2D
* arrays.
*
* matrix matrix to create
* buffer buffer to use for matrix elements. Its size should conform to
nrow and ncol
* nrow number of rows in buffer and hence number of rows in matrix
* ncol number of columns in buffer and hence number of columns in
matrix
*
*/
int matrix.new.from.buffer (matrix.t ^matrix float **buffer unsigned int nrow.
unsigned int ncol)
{
MATRTX_ASSERT( matrix != NULL);
86


N1ATRJX_VERIFYJ>REC0ND (( nrow != 0) &&
(ncol != 0).
MATRIX-ERRORJNVALID_DIMENSIONS) ;
matrix>nrow = nrow ;
matrix>ncol = ncol ;
matrix>matrix = (real.t **)buffer;
return MATRIX_ERROR_NONE;
}
/*
* matrix.free_from_buffer ()
*
* Dual of matrix_new_from_buffer () .
*
* matrix matrix to free.
*
. */
int matrix_free_from_buffer ( matrix.t *matrix)
{
MATRDC_ASSERT( matrix != NULL);
matrix>ncol 0;
matrix>nrow = 0;
matrix>matrix = NULL;
return MATRIX_ERROR_NONE;
}
/ *
* Actually create the memory pool matrices .
*/
static int _matrix_poo 1.init ( matrix.poo 1 _t *pool unsigned int count)
{
int 1 j ;
matpool = malloc ( s i z e o f ( s t r u c t pool_elem_s) count);
matpool.count = count ;
for ( i =0; i matpool [ i ] nrow = pool [ i ] nrow ;
matpool [ i ] ncol = pool [ i ] ncol ;
matpool[i ] count = pool[i ] count;
matpool [ i ]. matrix = (matrix_t *) malloc ( s i z eof ( matrix_t ) * pool [ i ]. count)
matpool [ i ]. free = (unsigned char *) malloc ( pool [ i ]. count) ;
memset ( matpool [ i ] free OxFF pool [ i ]. count) ;
for (j =0;j _matrix_new (&mat pool [i].matrix[j], pool [ i ] nrow pool[i].ncol);
}
}
return MATRIX_ERROR_NONE;
}
/ *
* matrix_pool_init()
*
* Create (on the heap) all of the desired matrices in
* pool 5 .
*
* pool sizes and counts of matrices to allocate
* count number of unique sizes (length of pool)
*
*/
int matrix.pool.init ( matrix.pool.t *pool, unsigned int count)
{
87


h
J!P3#
(o(* P!OA )) TlflN au!J3P#
THIN J 3 P uJ !#
apnpui#
H~XIHXVM su'jap#
irxraxvM jspuj!#
raoD jiBUiSgiUosipBf' Auoqjub
jb ioi])nt aqj p)aoo aiBMjjos siqj asuaaij oj_ sjuamaiinbaj
sjt jo juaui[[ij[nj aqj joj puB sasodind uoijBnjBAa joj pjBuSig
oipny jo uoijBJBdag aojnog pui]g joj yQjjsBg J uoi jb juauiajduij
auiijpay y.^ ui pajsi] jnq XiBjaijdoid si ajBMjjos siqx

uosqoi3f y Xuoqiuy §00S (3) ^qSu^doQ *
:*****************************************/
q-xij jbui/d_o; jsbj £2y
{
!3NON^IOHaa_5aaXVK uanjaj
:TTT1N = podjBui
! (podjBui ) aa JJ
!0 = junoa-pod jbui
{
irTCTLN = xiijbui [ i ] podjBui
TITLN = a8JJ [ ] podjBui
! ( xijjbuj [ i ] podjBui) aajj
i ( aajj j i j podjBui) aajj
{
i ([ f ] xij^^ui [ i ] jood!j'Biu2^)88jj"xta^'eui-
} (++f! ^unoo [ i ] ioodi^ui>po = f) ioj
} (++i t iunoDtoodi'eui>pQ= ] ) J0J
!f 4 i }ui
}
(p!OA)Xojjsap-[ood-xijjBm j u i
/*
*
[ood Xjouiaui aqj u; saaiijBui jjb aaJjj *
*
() Xo i jsap-]ood~xi i jbui *
*/
{
i(junoa [ood ) jiupjood-xii jbui- ujnjaa
(sNoiSNawicranvANrHOHHaxraxvpv p^q i) oNO03HcTAaraaA"xrHxvw
{
{
!quajq
l = p^q
} ((0 = 13U [ ]Id )
j | (o = *om j i ] [ood )
|| (0 = junoa [i]|Ood)) ji
} (+4-i : junoo>i!o= ) JOj
/* sjajaui'eaBd pauoo s'eq food Xjijqa ^siij */
(sNOiSNaivicranvANrHOHHa'xraxviM
o =i ;unoa)aNOoaad_AxraaAxraxvi'M
; (THIN =i I ood )XH3SSVXIHXVW
! 1ui
0 = pvq ;ui