# How to perform Blind Source Separation on an audio mix?

There have been several scholarly publications on the Blind Source Separation (BSS) problem. Some are aimed at audio signals. I wish to start with a simple case.

I have generated a custom sound mix of the "C" notes of 4 different instruments.

soundmixed =
Sound[{SoundNote["C", {0, 4}, "Oboe"],
SoundNote["C", {0, 4}, "SynthVoice"],
SoundNote["C", {0, 4}, "Organ"],
SoundNote["C", {0, 4}, "VoiceAahs"]}];


Before I apply Fourier transform onto it I need to find the corresponding audio data of the mixed sound and the sampling rate.

Q1. How can I get this information?

Let data be my sample data and fs be the sampling rate, I should be able to extract the frequency components in the following way:

nl=fs/2;
ft = Fourier[data, FourierParameters -> {-1, -1}];
temp = Drop[Drop[Chop[ft], -nl], 1];
freq = Flatten[SparseArray[Abs[temp]]["NonzeroPositions"]];
amp = SparseArray[Abs[temp]]["NonzeroValues"]*2;
pha = SparseArray[Arg[temp]]["NonzeroValues"];


Q2. After Fourier transform, is it possible to recognise the frequency spectra for these four different instruments?

Edit 1: A great and illustrative answer is given by Anton Antonov using Independent Component Analysis (ICA). Following his answer, I was looking into the spectrograms of the two instruments he has considered for mixing and I found that the frequency values are spread all over the range. However, a specific note should have a specific set of harmonics - fundamental and the overtones. Interestingly, the overtones should be multiples of the fundamentals. Each instrument has its unique set of harmonics.

I prefer to separate the sounds based on the analysis of the harmonics (or timbre) similar to this (https://www.projectrhea.org/rhea/index.php/Fourier_analysis_in_Music) instead of PCA.

• It should be possible to look at higher harmonics of your C note and then solve system of linear equations to figure out what instruments were contributing? Somewhat related Jan 9, 2017 at 17:33

In this post the first question is fully answered, the second question is answered using Independent Component Analysis, not Fourier transform.

For further details on using Independent Component Analysis (ICA) in Mathematica see:

## Answer of the first question

The answer of the first question:

Q1. How can I get this information?

is to use Audio[_,SampleRate->_] and AudioData.

soundmixed =
Sound[{SoundNote["C", {0, 4}, "Oboe"], SoundNote["C", {0, 4}, "SynthVoice"],
SoundNote["C", {0, 4}, "Organ"], SoundNote["C", {0, 4}, "VoiceAahs"]}];

fs = 44000;
smAus = Audio[soundmixed, SampleRate -> fs];


Extract data:

smData = AudioData[smAus];
Dimensions[smData]
(* {2, 177088} *)


## Answer of a generalized and simplified version of the second question

Consider the recovery of the instruments notes using two different microphones (and recorders) in the room -- i.e. we consider the cocktail party effect.

Below it is shown how to recover the signals using Independent Component Analysis (ICA) using two instruments. Note that no Fourier analysis is used.

I did apply this procedure to other sets of mixed instruments (mostly pairs). The ICA separation was very successful fair amount of times, and also not so good fair amount of times. Re-calculating ICA might give better results. Experiments with different mixture ratios (specified with SoundVolume) produce noticeably diffrent results.

### Set-up

Time length and sampling frequency parameters:

tl = 1.5;
fs = 44000;


Instrument data:

instruments = {"Oboe", "SynthVoice", "Organ", "VoiceAahs", "Guitar"};
aInstrumentAus =
Association[# -> Audio[SoundNote["C", {0, tl}, #], SampleRate -> fs] & /@
instruments];


Select couple of instruments and make a cocktail party mix :

mixInstruments = {"Organ", "SynthVoice"};
cpSoundmixes = Map[Sound[{
SoundNote["C", {0, tl}, mixInstruments[[1]], SoundVolume -> #[[1]]],
SoundNote["C", {0, tl}, mixInstruments[[2]], SoundVolume -> #[[2]]]
}] &, {{1, 3/4}, {1/2, 1}}]


Sample the cocktail party mixtures:

cpAus = Audio[#, SampleRate -> fs] & /@ cpSoundmixes;


Extract data:

aInstrumentData = AudioData /@ aInstrumentAus;
Dimensions /@ aInstrumentData

(* <|"Oboe" -> {2, 67072}, "SynthVoice" -> {2, 67072}, "Organ" -> {2, 67072},
"VoiceAahs" -> {2, 67072}, "Guitar" -> {2, 67072}|> *)

cpData = AudioData /@ cpAus;
Dimensions[cpData]

(* {2, 2, 67072} *)


### Apply ICA

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/IndependentComponentAnalysis.m"]
SeedRandom[342];
X = Transpose[cpData[[All, 1]]];
icaRes = IndependentComponentAnalysis[X, 2, PrecisionGoal -> 5];
Norm[X - Dot @@ icaRes]

(* 2.10679*10^-14 *)


Plot ICA results:

Block[{data = Transpose[icaRes[[1]]], step = 100},
opts = {PlotRange -> All, ImageSize -> Medium,
PlotTheme -> "Detailed"};
ListLinePlot[#[[1 ;; -1 ;; step]], opts] & /@ data]


### The first instrument

Next we plot/show the original, ICA extracted, and mixed signals for each instrument.

ind = 1; icaInd = 1; cpInd = 1;
Grid[{{mixInstruments[[1]] <> " channel 1", "ICA extracted " <> ToString[icaInd], "recorder " <> ToString[cpInd] <> " channel 1"},


{Sound[SampledSoundList[aInstrumentData[mixInstruments[[ind]]][2],fs]], Sound[SampledSoundList[ Rescale[icaRes[[1, All, icaInd]], MinMax[icaRes[[1, All, icaInd]]], {-1, 1}], fs]], Sound[SampledSoundList[ Rescale[cpData[[cpInd, 1]], MinMax[cpData[[cpInd, 1]]], {-1, 1}],fs]]}}]

### The second instrument

ind = 2; icaInd = 2; cpInd = 2;
Grid[{{mixInstruments[[ind]] <> " channel 1", "ICA extracted " <> ToString[icaInd], "recorder " <> ToString[cpInd] <> " channel 1"},


{Sound[SampledSoundList[aInstrumentData[mixInstruments[[ind]]][2],fs]], Sound[SampledSoundList[-Rescale[icaRes[[1, All, icaInd]], MinMax[icaRes[[1, All, icaInd]]], {-1, 1}], fs]], Sound[SampledSoundList[ Rescale[cpData[[cpInd, 1]], MinMax[cpData[[cpInd, 1]]], {-1, 1}],fs]]}}]

Listen to the results -- the ICA extractions are noticeably cleaner than the recorded ones.

• Audio is new to MMA 11. Are you aware how to get the data out of Sound in v10? Jan 9, 2017 at 22:02
• @BlacKow I am sorry, I don't know. I wonder are there some hints how to do that in WTC talks for few years back... Jan 10, 2017 at 1:51
• @AntonAntonov : You have directly used the value of the sampling frequency as 44000. I know this is a standard value for audio signals. However, it may be possible (I'm not sure) to sample the signal at some other frequency. In that case, is there any generalised way to determine the sampling frequency programatically? Jan 10, 2017 at 12:43
• @AntonAntonov : While creating the cocktail party mix you have used different SoundVolumes. Is it only to recognise the sounds in the mix better or it has some other purpose? Jan 10, 2017 at 13:18
• @Majis "While creating the cocktail party mix you have used different SoundVolumes. Is it only to recognise the sounds in the mix better or it has some other purpose?" -- Generally speaking, if the microphones have the same volumes of the instruments they record, then ICA cannot be applied. (We have multiple recordings of the same signal.) With the sound files produced in this discussion if the volumes are the same for all recordings we can apply ICA, though, since we have two channels for each instrument, and ICA can be used to find a "mean" signal for one of the instruments. Jan 10, 2017 at 17:03

I recently discovered this interesting problem, known as cocktail party problem. My understanding is that spectacular source separations, such as these, can presently be achieved with neural network or equally elaborate techniques. However I'd like to summarize my understanding of a much more basic algorithm based on Independent Component Analysis, which is OK for simpler examples. Anton Antonov's answer proposed an implementation of this algorithm, presumably much more sophisticated than what follows, but I wanted to try by myself from scratch.

## The Problem

There are two independent sound sources. We could record them independently and stack the results in a vector $X_0$. But instead, we recorded them while both sources were playing, resulting in a measurement

$$X=AX_0$$ where $A$ is a $2\times 2$ matrix called mixing matrix. The objective is to estimate $X_0$ from the measurements $X$ only.

Source generation

For the first example, let's generate two sines with different frequencies and amplitudes:

X0 = {1.3*Table[Sin[2345.*t + 0.1], {t, Subdivide[2, 10000]}],
Table[Sin[2000.*t], {t, Subdivide[2, 10000]}]};
fs = 5000;
Audio[#, SampleRate -> fs] & /@ X0


We can check that the source signals are far from Gaussian, and that will help.

GraphicsRow[Histogram /@ X0]


We can also check that the covariance is almost diagonal, meaning the source signals are decorrelated:

Covariance[Transpose@X0]
(* {{0.845071, -0.000948025}, {-0.000948025, 0.499964}} *)


As a remainder, covariance the matrix whose coefficient $i,j$ is

$$\dfrac{1}{n-1} \sum_{k=1}^n (x_i - \overline{x}_i)(x_j-\overline{x}_j)$$

Mixing source signals

A = {{.3, .7}, {1.6, 1}};
X = A.X0;
audios = Audio[#, SampleRate -> fs] & /@ X


## Resolution

Whitening (or sphering) While the source were uncorrelated, we can check the measurements are not:

Covariance[Transpose@X]
(* {{0.320641, 0.754263}, {0.754263, 2.66031}} *)


So the first step is to uncorrelate the data. This is done by left-multiplying $X$ with the inverse of the square root of the correlation matrix:

Y = MatrixPower[Covariance[Transpose@X], -.5].(X - Mean /@ X)


This operation has a geometrical interpretation if we display $X_0$, $X$, or $Y$ as a list of points:

GraphicsRow[ListPlot[Transpose@#, PlotRange -> Full, AspectRatio -> Automatic,
ImageSize -> Small] & /@ {X0, X, Y}]


The matrix $A$ transformed the rectangle (of dimensions the sine amplitudes) as a diamond, then the uncorrelation transformed the diamond into a (rotated) square.

Rotating Now we need to rotate $Y$ in order to recover $X_0$, or more precisely something in the form diagonal matrix times $X_0$ (indeed, the present algorithm is not unique, there is no control on the amplitude of the estimated sources).

In the present case, the angle could be found easily, but in real life example $Y$ does not correspond to a rectangle. So the idea, and that's the key idea of this method, is to find the angle that correspond to "least Gaussian" signals. It relies on the Central Limit theorem: a linear combination of two independent random variables is more Gaussian than the random variables themselves.

This can be seen below: the linear combinations $X$ and $Y$ are closer to Gaussians than $X0$:

GraphicsGrid[Histogram /@ # & /@ {X0, X, Y}]


So now the question is how to measure "Gaussianity"? There are several techniques, some based on negentropy, some based on kurtosis. Kurtosis is the fourth moment, and is one measure of Gaussianity. A normal distribution has a kurtosis of 3, so we can define a cost function as:

J[theta_] = Kurtosis /@ (RotationMatrix[theta].Y) - 3 // Simplify


Then we maximize the non-Gaussianity to obtain $Z$, our estimation of $X_0$:

thetaOpt = theta /. Last@FindMaximum[Norm@J[theta], {theta, 1.2}]
Z = RotationMatrix[thetaOpt].Y;


The estimation of $A$ is:

Atilde = Inverse[RotationMatrix[thetaOpt].MatrixPower[Covariance[Transpose@X], -.5]


To check the results:

 Audio[#, SampleRate -> fs] & /@ # & /@ {X0, X, Z}


It's not perfect, but not too bad!

Other examples

You can try by mixing the instruments, as in Anton Antunov's answer:

tl = 3;
fs = 44000;
audio = Audio[#, SampleRate -> fs] & /@ {SoundNote["C", {0, tl},
"Piano", SoundVolume -> 1],
SoundNote["C", {0, tl}, "SynthVoice", SoundVolume -> 1]};
X0 = (AudioData /@ audio)[[All, 1]];
A = {{1.2, .6}, {.4, 1.1}};
X = A.X0;
Audio[#, SampleRate -> fs] & /@ X


The results is not too bad either, see the histograms (note that the mix, in the second line, is pretty close to a Gaussian!):

GraphicsGrid[Histogram[#, 100] & /@ # & /@ {X0, X, Z}]


You can also try this cost function:

J2[theta_] = With[{y = RotationMatrix[theta].Y},
1/12*Mean[#^3]^2 + 1/48*Kurtosis[#]^2 & /@ y] // Simplify


Unfortunately, this code yields very poor results on a more complex examples such as section 1.

References