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;cncol ; 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;rnrow ; r++) {

for ( c =0;cncol ; 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;rnrow ; r++) {

for ( c = 0;cncol ; 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;rnrow ; r++) {

for ( c =0;cncol ; 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