Description of Source Separation Algorithm
The approach that I would take is based on a Semi-supervised Non-negative Matrix Factorization. This would work assuming that:
- The audio file that you are playing back is the same every time.
- You are able to capture the played back recording on its own and, without background noise (i.e. man talking).
- The algorithm will work best in the same room/environment and with the same setup of source-recorder.
In that case, all you have to do is:
- Record the playback signal $x_p$.
- Calculate the amplitude spectrogram $|X_p|$ of it - you can ignore the phase. It will have a size of $N\times T$, where $N$ is the number of frequency bins and $T$ is the number of frames.
- Decompose the $X_p$ using NMF of some rank $R_p$. This will give you dictionary $W_p$ and weights $H_p$. The dictionary will be of size $N\times R$, and weights will have size $R\times T$.
Now you are able to reconstruct the original spectrogram simply by doing the multiplication $\tilde{X_p} = W_p H_p$. The basis matrix contains the atoms, which are basically a frequency patterns that describe nature of your signal. For the final part:
- Record the signal $x$ that contains both the playback signal $x_p$ and unknown speech signal $x_s$.
- Extract the spectrogram $|X|$ similarily to the previous step. This time please store the phase $\angle X$ - it will be required for the reconstruction.
- Assuming the rank of $R_s$ for speech, create the new, concatenated basis and weight matrices, where $H_p$ is the same as calculated before and $W_p$, $W_s$, $H_s$ are random numbers:
$$ H = [H_p, H_s]$$
$$ W = [W_p; W_s]$$
The sizes are:
- $H_s$: $N\times R_s$
- $W_s$: $R_s \times T$
- $H$: $N\times (R_p + R_s)$
- $W$: $(R_p + R_s) \times T$
Now:
- Apply the partial updates to learn the $H_p$, $H_s$ and $W_s$. The basis $W_p$ stays fixed!
- After the decomposition is done, create the masks for playback and speech:
$$M_p = \frac{W_pH_p}{WH}$$
$$M_s = \frac{W_sH_s}{WH}$$
- Multiply element-wise the spectrogram $X$ by the mask $M_p$ to get the playback part $\tilde{X}_p=X\circ M_p$ and do the same to obtain $\tilde{X}_s=X\circ M_s$.
- Using the stored phase $\angle X$, perform ISTFT (inverse spectrogram) on both $\tilde{X}_p$ and $\tilde{X}_s$, which will return the time domain signals $\tilde{x}_p$ and $\tilde{x}_s$.
In a nutshell, how the algorithm works:
- You learn a dictionary with contains patterns specific to your signal $X_p$ via NMF decomposition.
- You create a new, bigger dictionary, which contains both pre-learned atoms that are fixed, as well as new atoms, which will be learned through the decomposition - these should effectively capture the nature of anything that is not $X_p$.
- During that decomposition step, weights in matrix $H$ are learned from the scratch.
- Once the whole factorization is finished, you can get the partial representation for both speech and playback and create the masks (through Wiener filtering).
Regarding factorization, it can be done either on a whole file or processed per blocks (for example 10 seconds). I would suggest to start with the whole file, play with the number of atoms $R_s$ and see how it works.
Semi-supervised NMF
All of that might sound very complicated, but in fact NMF is extremely easy to implement. Additionally there are many packages out there. However, for sake of completeness I am going to provide the full and partial update rules and roughly explain the algorithm.
Since there are many flavours of NMF with various divergences, I will focus on Kullback-Leibler (KL) and Itakura-Saito (IS) NMF. The latter one is in theory a better fit for spectrogram decomposition.
Core algorithm
Training
Inputs:
- $X_p$: magnitude spectrogram of the training signal, size $N\times T$
- $R_p$: rank of the decomposition (number of atoms to use), should be less than $T$. This is the parameter to be optimized.
- $\alpha$: sparsity penalty, should be in range $0 - 1$. $0$ means that there is no penalty - a good starting point. This parameters enforces sparse weight activations.
- $I$: maximum number of iterations for factorization. This is the stopping criterion.
Outputs:
- $W_p$: basis matrix, which contains the atoms describing the training pattern, size $N \times R_p$
Algorithm:
- Check that all values of the $X_p$ are positive.
- Add a small positive random numbers to $X_p$ to eliminate any zeros.
- Initialize $H_p$ and $W_p$ with positive random values (uniform distribution between $0$ and $1$).
- Add small random offsets to $H_p$ and $W_p$.
- Normalize columns of $W_p$ by their sum.
- Update $H_p$ with sparsity $\alpha$.
- Update $W_p$
- Normalize columns of $W_p$ by their sum.
- Repeat steps 6-8 for number of iterations $I$.
- Return basis $W_p$.
Factorization
Inputs:
- $X_s$: magnitude spectrogram of the test signal, size $N\times T$
- $R_p$: number of atoms for input pattern
- $R_s$: number of atoms for anything else
- $W_p$: basis matrix with trained atoms for the input pattern
- $\alpha$: sparsity penalty, should be in range $0 - 1$. $0$ means that there is no penalty - a good starting point. This parameters enforces sparse weight activations.
- $I$: maximum number of iterations for factorization. This is the stopping criterion.
Outputs:
- $W$: basis matrix containing both atoms for $W_p$ target pattern and other signal $W_s$
- $H$: weight matrix containing both activations for $H_p$ and other signal $H_s$
Algorithm:
- Check that all values of $X_s$ are positive
- Add a small random offset to $X_s$
- Create the basis matrix $W$ of size $N \times (R_p + R_s)$ and weight matrix $H$ of size $(R_p + R_s) \times T$.
- Initialize the $H$ and $W$ with random numbers (uniform distribution).
- Add a small random offset to $H$ and $W$.
- Set the first $R_p$ columns of $W$ to $W_p$.
- Normalize columns of $W_p$ by their sum.
- Partially update $W_p$
- Partially update $W_s$
- Partially update $H_s$
- Normalize the columns of $W$ by their sum
- Repeat steps 8-11 for $I$ iterations
Update rules
The $\circ$ means Hadamard product (element-wise multiplication).
$\mathcal{I}$ is the identity matrix (filled with ones) having the same size as $X$.
Division in terms of matrices is an element-wise division.
Raising to the power is element-wise opration.
Partial updates for $s$ are the same as for $p$.
Kullback-Leibler NMF
$H$ update:
$$H \leftarrow H \circ \dfrac{ H^T \dfrac{X} {WH} + \alpha}{H^T\mathcal{I}}$$
$W$ update:
$$W \leftarrow W \circ \dfrac{ \dfrac{X} {WH} W^T}{\mathcal{I}W^T}$$
Partial $H$ update:
$$ H_p \leftarrow H_p \circ \dfrac{W_p^T \dfrac{X}{WH} + \alpha}{W_p^T\mathcal{I}} $$
Partial $W$ update:
$$ W_p \leftarrow W_p \circ \dfrac{\dfrac{X}{WH} H_p^T }{\mathcal{I} H_p^T} $$
Itakura-Saito NMF
$H$ update:
$$ H \leftarrow H \circ \dfrac{W^T \dfrac{X}{[WH]^2} + \alpha}{W^T\dfrac{1}{WH}} $$
$W$ update:
$$ W \leftarrow W \circ \dfrac{\dfrac{X}{[WH]^2} H^T }{\dfrac{1}{WH} H^T} $$
Partial $H$ update:
$$ H_p \leftarrow H_p \circ \dfrac{W_p^T \dfrac{X}{[WH]^2} + \alpha}{W_p^T\dfrac{1}{WH}} $$
Partial $W$ update:
$$ W_p \leftarrow W_p \circ \dfrac{\dfrac{X}{[WH]^2} H_p^T }{\dfrac{1}{WH} H_p^T} $$
Cost functions
These are useful if you want to debug your code or maybe include yet another stopping criterion based on a minimum change between previous and current value of the cost function. The sum $\Sigma$ is a sum of all elements.
Kullback-Leibler
$$D_{KL}(X, WH)=\Sigma\left[ X\log \left(\dfrac{X}{WH}\right)-X+WH\right]$$
Itakura-Saito
$$D_{IS}(X, WH)=\Sigma \left[ \dfrac{X}{WH} -\log \left( \dfrac{X}{WH}\right) -1 \right]$$
