# Remove a known wav file from recorded file

I looked at many similar questions but still cannot find an answer to my problem.

I'm playing an audio file N, while recording a man talking (b). The problem is that the microphone gets the played sound N together with the man talking, b. so actually I have y=N*+b in my recording, where N* is the played sound with a short time delay and a bit noise (due to the recording and the phone speaker).

My goal is to totally remove the N* from the recorded file. It seems like an easy problem cause I have the original played file N, but I find it more complicated.

What I tried:

1. finding the delay between the N and y, invert N and add to y. Seems like it's not working and the result is even worse probably because N* has different amplitudes and a phase shift.
2. Tried the same thing but subtracting in the frequency domain between y and N. Still no good.

Any suggestions?

• Well, this would work if you actually detected the delay including the phase shift (unless effects like time-variant channels and sampling rate mismatches come into play), and also the amplitude change. So, how do you go about "finding the delay between N and y", exactly? – Marcus Müller Nov 25 '17 at 18:27
• I use the matlab commands "finddelay" and "alignsignals" to align the two wav files. It finds the delay and aligns them correctly. But as I said, probably the recorded N* has more differences, not only a time shift, causing the subtraction (after the alignment) to be useless. – user2630165 Nov 25 '17 at 18:54
• huh, never heard of ´finddelay´ nor ´alignsignals´. Does finddelay give you a real number, or just an integer sample offset? As you said, if your delay estimate doesn't include a phase shift, this is not going to work. The usual way to go about this is working with the position and phase and magnitude of the maximum of the crosscorrelation function; the position gives you the integer samples delay, the phase shift can then be used e.g. with a MMSE resampler or just a complex rotator (that won't help much with your real-valued signals), and the magnitude to normalize the signals. – Marcus Müller Nov 25 '17 at 18:57
• Finddelay simply gives a sample offset and then I can align the samples (only by time). Do you have any idea why the recorded audio has a phase shift so trivial samples subtraction cannot work? Is there another way to solve this problem ? – user2630165 Nov 25 '17 at 19:19
• I've explained the solution. And as you noticed, the delay that N experiences is not inherently an integer number of samples – Marcus Müller Nov 25 '17 at 19:28

Assuming that there are no non-linear effects between the sound N source and the microphone (such as AGC), you might try to estimate the impulse response of the channel between the source for sound N and the microphone data. The better estimate you can make of that impulse response, the more of sound N that you can remove. It’s unlikely that the negative of the impulse response is just the data inverted. That's due to the fact that the amplifier, speaker, microphone, and room acoustics (etc.) probably have frequency dependent phase delays and/or multiple reverb/echoes.

Channel estimation in noise is a non-trivial problem.

This is not just a channel estimation problem. Because the known reference signal and the recording have not been clock-synchronised you will have a small amount of temporal drift between the clocks. This makes the relationship between recording and reference non-trivial. If you can control the recording setup you should therefore synchronise the playback clock with the recording clock. If you cannot, then you need to first estimate the clock drift and compensate for it using some form of time-warping.

With synchronised clocks you do in fact have a channel estimation problem, but the channel will in general time-variant. That means the channel parameters cannot be estimated globally but need to be found locally. A good approach for such problems is a matching pursuit with a dictionary of the family of time-shifted copies of the reference recording. It is straight forward to show that such a matching pursuit would with each step reduce the error until you hit the level of incoherence (where all inner products of the family are more or less of the same magnitude) where you should stop. Taking into account the generally time-varying nature of the channel, you should not follow the matching pursuit globally but for shorter sections of the reference signal. The length of the section needs to be traded off against the channel variation. Shorter segments will adapt more quickly to channel changes, longer segments will reject noise better. The segments should be chosen so that they overlap and add up to the reference signal. Overlapping triangular windows are suitable candidates for this.

# 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:

1. Check that all values of the $X_p$ are positive.
2. Add a small positive random numbers to $X_p$ to eliminate any zeros.
3. Initialize $H_p$ and $W_p$ with positive random values (uniform distribution between $0$ and $1$).
4. Add small random offsets to $H_p$ and $W_p$.
5. Normalize columns of $W_p$ by their sum.
6. Update $H_p$ with sparsity $\alpha$.
7. Update $W_p$
8. Normalize columns of $W_p$ by their sum.
9. Repeat steps 6-8 for number of iterations $I$.
10. 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:

1. Check that all values of $X_s$ are positive
2. Add a small random offset to $X_s$
3. 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$.
4. Initialize the $H$ and $W$ with random numbers (uniform distribution).
5. Add a small random offset to $H$ and $W$.
6. Set the first $R_p$ columns of $W$ to $W_p$.
7. Normalize columns of $W_p$ by their sum.
8. Partially update $W_p$
9. Partially update $W_s$
10. Partially update $H_s$
11. Normalize the columns of $W$ by their sum
12. 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]$$

• May I ask in technical question? For sparsity at KL-NMF, is it common to add $\alpha$ ratio directly? I used to see $\alpha$ as a function and it is added at partial H update only. – Jan Mar 30 '18 at 6:10

What you describe as N* is actually a little more complicated than just a delayed N + noise. N* also contains the result of the convolution that happens along the signal path, if you consider the path as an LTI system. The impulse response with which N gets convolved is the result of your D/A converter that is connected to your speakers which are operating in a room and which are picked up by a microphone that is connected to an A/D converter. All of these elements contribute to the impulse response and therefore to N*.

Unfortunately, it's also nearly impossible to obtain the exact impulse response N was convolved with in your recording. It is affected by too many factors for there to be a chance of actually replicating it. Otherwise it would be trivial—just do another recording but without someone talking, invert this new recording and subtract it from the first recording. That however most likely won't work since it's not feasible to recreate the exact situation in which the first recording happened.

So to answer your question: There is no surefire way of doing what you're set out to do. If you have a lot of control over the circumstances in which the recording takes place (anechoic room, reliable hardware with known characteristics), there might be a chance. For "live" situations, there is little hope.

• Thanks for the answer. Its very interesting. But if I understand correctly, you're actually saying it will be a waste of time to try out Marcus Müller suggestion? – user2630165 Nov 25 '17 at 20:15
• @user2630165 If you're working solely with N when trying to fully remove N* from the recording, then no, that definitely won't work (even when disregarding the aspect of noise that is also contained in N*). As stated above, the situation is less dire the more control you have over the recording siuation, but realistically you shouldn't expect being able to isolate the speech signal from the recording. – S. R. Nov 25 '17 at 20:24

So I actually just spent the last year working on a website that does just this, www.tuneout.io. I just launched a free Beta version. Curious to hear if it works for you.

Since SO doesn't like here's-a-website-that-exactly-solves-your-problem answers, I'll address some points brought up by other commenters. The problem is mainly the channel distortion on the audio (room acoustics, speaker distortion, etc.). Phase inversion only really works if the sound is very clean, with little to no reverberation, the speaker and mic have flat frequency responses, and the volume is constant. Doing the subtraction in the frequency domain is actually mathematically the same as subtraction in the time domain (assuming no windowing). So it's expected you'll get the same result. Another issue is that the song in the recording will often be a different length than your audio file of the song. This is because real-life devices all have slightly different playback speeds due to slight hardware differences. So, even if you find the exact delay using finddelay(), if there is "time drift" in the file, the delay will also drift, quickly rendering phase inversion useless. Add to this you could be starting with a file of slightly different length due to compression (mp3s can be several seconds longer than a wav file sometimes).

My algorithm uses a special blend of spectral subtraction and machine learning to address all of these problems in a robust (hopefully) fully automated way. Hope it works for you!