Audio Source Separation with Sparsity¶
Separate several audio sources from a smaller number of mixtures using sparse time-frequency structure. Examine spectrograms, estimate mixing directions, and listen to reconstructed signals to understand both the promise and limitations of masking.
Run this tour¶
Run the cells in order with a Python 3 kernel. The first cell locates the companion data and toolbox and installs missing dependencies when needed. All worked examples include their implementation directly in this notebook. Random seeds make comparisons reproducible; you can change them to explore other samples.
# Locate the companion toolbox locally, or fetch it for a standalone/Colab copy.
from pathlib import Path
import importlib.util
import os
import subprocess
import sys
working = Path.cwd()
candidates = [working, working / "python", working.parent / "python"]
python_dir = next((p for p in candidates if (p / "nt_toolbox").is_dir()), None)
if python_dir is None:
checkout = working / "numerical-tours-support"
if not checkout.exists():
subprocess.run(
[
"git",
"clone",
"--depth",
"1",
"--branch",
"master",
"https://github.com/gpeyre/numerical-tours.git",
str(checkout),
],
check=True,
)
python_dir = checkout / "python"
os.chdir(python_dir)
if str(python_dir) not in sys.path:
sys.path.insert(0, str(python_dir))
requirements = python_dir / "requirements.txt"
if any(
importlib.util.find_spec(name) is None
for name in [
"numpy",
"scipy",
"matplotlib",
"skimage",
"sklearn",
"pywt",
"ipywidgets",
"cvxpy",
"skfmm",
"autograd",
"progressbar",
"celer",
]
):
subprocess.run(
[sys.executable, "-m", "pip", "install", "-r", str(requirements)], check=True
)
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
plt.rcParams.update(
{
"figure.figsize": (8, 4),
"figure.dpi": 100,
"axes.spines.top": False,
"axes.spines.right": False,
"font.size": 11,
"image.cmap": "gray",
}
)
%matplotlib inline
$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This numerical tour explore local Fourier analysis of sounds, and its application to source separation from stereo measurements.
import numpy as np
import scipy as scp
import pylab as pyl
import matplotlib.pyplot as plt
from nt_toolbox.general import np, plt, pylab
from nt_toolbox.signal import imageplot, np, plt, pylab
import warnings
%matplotlib inline
Sound Mixing¶
We load 3 sounds and simulate a stero recording by performing a linear blending of the sounds.
Sound loading.
from nt_toolbox.load_sound import load_sound, np
n = 1024 * 16
s = 3 # number of sounds
p = 2 # number of micros
x = np.zeros([n, 3])
x[:, 0] = load_sound("nt_toolbox/data/bird.wav", n)
x[:, 1] = load_sound("nt_toolbox/data/female.wav", n)
x[:, 2] = load_sound("nt_toolbox/data/male.wav", n)
Normalize the energy of the signals.
x = x / np.tile(np.std(x, 0), (n, 1))
We mix the sound using a $2\mathrm{x}3$ transformation matrix. Here the direction are well-spaced, but you can try with more complicated mixing matrices.
Compute the mixing matrix
theta = np.linspace(0, np.pi, s + 1)[:-1]
theta[0] = 0.2
M = np.vstack((np.cos(theta), np.sin(theta)))
Compute the mixed sources.
y = np.dot(x, np.transpose(M))
Display of the sounds and their mix.
plt.figure(figsize=(10, 10))
for i in range(s):
plt.subplot(s, 1, i + 1)
plt.plot(x[:, i])
plt.xlim(0, n)
plt.title("Source #%i" % (i + 1))
Display of the micro output.
plt.figure(figsize=(10, 7))
for i in range(p):
plt.subplot(p, 1, i + 1)
plt.plot(y[:, i])
plt.xlim(0, n)
plt.title("Micro #%i" % (i + 1))
Local Fourier analysis of sound.¶
In order to perform the separation, one performs a local Fourier analysis of the sound. The hope is that the sources will be well-separated over the Fourier domain because the sources are sparse after a STFT.
First set up parameters for the STFT.
w = 128 # size of the window
q = w // 4 # overlap of the window
Compute the STFT of the sources.
from nt_toolbox.perform_stft import np, perform_stft, pyl
from nt_toolbox.plot_spectrogram import np, plot_spectrogram, plt
X = np.zeros([w, 4 * w + 1, s], dtype=complex)
Y = np.zeros([w, 4 * w + 1, p], dtype=complex)
for i in range(s):
X[:, :, i] = perform_stft(x[:, i], w, q, n)
plt.figure(figsize=(15, 10))
plot_spectrogram(X[:, :, i], "Source #%i" % (i + 1))
Worked example 1
Compute the STFT of the micros, and store them into a matrix |Y|.
for i in range(p):
Y[:, :, i] = perform_stft(y[:, i], w, q, n)
plt.figure(figsize=(15, 10))
plot_spectrogram(Y[:, :, i], "Source #%i" % (i + 1))
Estimation of Mixing Direction by Clustering¶
Since the sources are quite sparse over the Fourier plane, the directions are well estimated by looking as the direction emerging from a point clouds of the transformed coefficients.
First we compute the position of the point cloud.
mf = np.shape(Y)[0]
mt = np.shape(Y)[1]
P = np.reshape(Y, (mt * mf, p))
P = np.vstack((np.real(P), np.imag(P)))
Then we keep only the 5% points with largest energy.
Display some points in the original (spacial) domain.
Number of displayed points.
npts = 6000
Display the original points.
from numpy import random
sel = random.permutation(n)
sel = sel[:npts]
plt.figure(figsize=(7, 5))
plt.plot(y[sel, 0], y[sel, 1], ".", ms=3)
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.title("Time domain")
plt.show()
Worked example 2
Display some points of $P$ in the transformed (time/frequency) domain.
sel = random.permutation(n)
sel = sel[:npts]
plt.figure(figsize=(7, 5))
plt.plot(P[sel, 0], P[sel, 1], ".", ms=3)
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.title("Transformed domain")
plt.show()
We compute the angle associated to each point over the transformed domain. The histogram shows the main direction of mixing.
import math
nrow = np.shape(P)[0]
Theta = np.zeros(nrow)
for i in range(nrow):
Theta[i] = math.atan2(P[i, 1], P[i, 0]) % np.pi
Display histogram.
nbins = 100
t = np.linspace(np.pi / 200, np.pi, nbins)
hist = np.histogram(Theta, t)
h = hist[0] / np.sum(hist[0])
t = t[:-1]
plt.figure(figsize=(7, 5))
plt.bar(t, h, width=np.pi / nbins, color="darkblue", edgecolor="darkblue")
plt.xlim(0, np.pi)
plt.ylim(0, np.max(h))
plt.show()
Worked example 3
The histogram computed from the whole set of points are not peacked enough. To stabilize the detection of mixing direction, compute an histogram from a reduced set of point that have the largest amplitude. Compute the energy of each point. Extract only a small sub-set.
d = np.sum(P**2, 1)
rho = 0.1
v = np.sort(d)
I = np.argsort(d)[::-1]
# transformed points
I = I[np.arange(1, round(rho * len(I)) + 1)]
P1 = P[I, :]
# compute Theta
nrow = np.shape(P1)[0]
Theta = np.zeros(nrow)
for i in range(nrow):
Theta[i] = math.atan2(P1[i, 1], P1[i, 0]) % np.pi
nbins = 200
hist = np.histogram(Theta, nbins)
h = hist[0] / np.sum(hist[0])
t = hist[1][:-1]
plt.figure(figsize=(7, 5))
plt.bar(t, h, width=np.pi / nbins, color="darkblue", edgecolor="darkblue")
plt.xlim(0, np.pi)
plt.ylim(0, np.max(h))
plt.show()
Worked example 4
Detect the direction $M_1$ approximating the true direction $M$ by looking at the local maxima of the histogram. First detect the set of local maxima, and then keep only the three largest. Sort in descending order.
s1 = np.hstack((np.arange(2, nbins + 1), np.array([nbins - 1]))) - 1
s2 = np.hstack((np.array([2]), (np.arange(1, nbins)))) - 1
I = np.where((h[s1] < h) & (h[s2] < h))
v = np.sort(h[I])
u = np.argsort(h[I])[::-1]
theta1 = t[I[0][u[:3]]]
M1 = np.vstack((np.cos(theta1), np.sin(theta1)))
print("--- M ---")
print(M)
print("--- M1 ---")
print(M1)
Separation of the Sources using Clustering¶
Once the mixing direction are known, one can project the sources on the direction.
We compute the projection of the coefficients Y on each estimated direction.
A = np.reshape(Y, (mt * mf, p))
Compute the projection of the coefficients on the directions.
C = abs(np.dot(np.transpose(M1), np.transpose(A)))
At each point $x$, the index $I(x)$ is the direction which creates the largest projection.
$I$ is the index of the closest source.
tmp = np.max(C, 0)
I = np.argmax(C, 0)
I = np.reshape(I, (mf, mt)) + 1
An additional denoising is achieved by removing small coefficients.
T = 0.05
D = np.sqrt(np.sum(abs(Y) ** 2, 2))
I = I * (D > T)
We can display the segmentation of the time frequency plane.
plt.figure(figsize=(15, 10))
imageplot(I[: mf // 2, :])
plt.imshow(I[: mf // 2, :], cmap=plt.get_cmap("jet"), interpolation="nearest")
plt.show()
The recovered coefficients are obtained by projection.
Proj = np.dot(np.transpose(M1), np.transpose(A))
Xr = np.zeros([w, 4 * w + 1, s], dtype=complex)
for i in range(s):
Xr[:, :, i] = np.reshape(Proj[i, :], (mf, mt)) * (I == i)
The estimated signals are obtained by inverting the STFT.
xr = np.zeros([n, s])
for i in range(s):
xr[:, i] = perform_stft(Xr[:, :, i], w, q, n)
One can display the recovered signals.
plt.figure(figsize=(10, 10))
for i in range(s):
plt.subplot(s, 1, i + 1)
plt.plot(xr[:, i])
plt.xlim(0, n)
plt.title("Estimated source #%i" % (i + 1))
One can listen to the recovered sources.
i = 1
from IPython.display import Audio
Audio(x[:, i], rate=15000)
Audio(xr[:, i], rate=15000)
References and further reading¶
Özgür Yılmaz and Scott Rickard. Blind Separation of Speech Mixtures via Time-Frequency Masking. 2004, IEEE Transactions on Signal Processing 52(7), 1830–1847. Separating audio sources through sparse time-frequency structure.
Aapo Hyvärinen and Erkki Oja. Independent Component Analysis: Algorithms and Applications. 2000, Neural Networks 13(4–5), 411–430. Statistical assumptions and identifiability in blind source separation.
Stéphane Mallat. A Wavelet Tour of Signal Processing: The Sparse Way. 2009, 3rd ed., Academic Press. Multiresolution analysis, sparse approximation, and wavelet algorithms.
Gabriel Peyré. Advanced Signal, Image and Surface Processing. 2010, course notes. A mathematical companion to the Numerical Tours.
Pauli Virtanen et al.. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. 2020, Nature Methods 17, 261–272. The numerical routines used for transforms, interpolation, and optimization.