Multivariate Wiener Filter Derivation

The Wiener filter is the optimal filter to use if you want to remove Gaussian noise from a Gaussian signal. Most useful applications though involve multiple (possibly correlated) signals with additive noise, so a multivariate Wiener filter is required. In this post I'll derive the multivariate Wiener filter as a Maximum A-Posteriori (MAP) estimator of the clean signal given the noisy signal. Because Gaussians are symmetric, the MAP estimator will be identical to the MMSE estimator, with the advantage that the MAP estimator is much easier to compute.

So our assumption is we have a noisy signal \(y(n)\) which is the result of adding gaussian noise \(d(n)\) to our clean signal \(x(n)\). Remember each of these are vectors, so \(y(n)\), \(d(n)\) and \(x(n)\) are vectors with \(k\) dimensions.

\[ y(n) = x(n) + d(n)\]

What we want to do is find the most likely value of \(x(n)\) given \(y(n)\). We never get to know \(d(n)\), but we do have an estimate of the covariance of \(d(n)\) which we will call \(\Sigma_d\). We assumed initially that \(x(n)\) and \(d(n)\) were distributed according to the multivariate Gaussian distribution, but because the pdf takes up so much space, I'll be using the following notation:

\[ a \sim \mathcal{N}(\boldsymbol\mu_a,\boldsymbol\Sigma_a) = (2\pi)^{-\frac{k}{2}}|\boldsymbol\Sigma_a|^{-\frac{1}{2}}\, e^{ -\frac{1}{2}(\mathbf{a}-\boldsymbol\mu_a)'\boldsymbol\Sigma_a^{-1}(\mathbf{a}-\boldsymbol\mu_a) } \]

This means \(a\) is distributed according to a multivariate Gaussian with mean \(\boldsymbol\mu_a\) and covariance \(\boldsymbol\Sigma_a\). So we have \(d(n) \sim \mathcal{N}(0,\boldsymbol\Sigma_d)\) and \(x(n) \sim \mathcal{N}(\boldsymbol\mu_x,\boldsymbol\Sigma_x)\).

From our intial assumtion that \(y(n) = x(n) + d(n)\), we can also determine the pdf of \(y(n)\) if we know \(x(n)\). This is simply \(p(y | x) = \mathcal{N}(x,\boldsymbol\Sigma_d)\), i.e. we know \(y\) will be centered around \(x\) with a variance due only to the noise. We now have the distribution of the noisy given the clean, unfortunately we won't know the clean signal. What we want is the distribution of the clean signal given the noisy signal. From here we can use Bayes Rule:

\[p(x | y) = \dfrac{p(y|x)p(x)}{p(y)} \]

This will give us the pdf of the clean signal given the noisy, then we just have to find the maximum of the distribution to find our estimate of \(x\). Note that \(p(y)\) in the above equation does not effect the maximisation (because it does not depend on \(x\)), so it will be left out. From here we write out everything we have, using the values for \(p(x)\) and \(p(y|x)\) we worked out above, hopefully fitting everything on the page:

\[ p(x | y) = p(y|x)p(x)\] \[=(2\pi)^{-\frac{k}{2}}|\boldsymbol\Sigma_d|^{-\frac{1}{2}}\, e^{ -\frac{1}{2}(y - x)'\boldsymbol\Sigma_d^{-1}(y-x)} \times (2\pi)^{-\frac{k}{2}}|\boldsymbol\Sigma_x|^{-\frac{1}{2}}\, e^{ -\frac{1}{2}(x-\boldsymbol\mu_x)'\boldsymbol\Sigma_x^{-1}(x-\boldsymbol\mu_x)}\]

We want to maximise the equation above, unfortunately it is rather complicated. To get around this we take the log of the whole thing, since maximising a function and maximising its log will give the same point. The result is much easier to differentiate:

\[ \log(p(x | y)) = \log(C) - \frac{1}{2}(y - x)'\boldsymbol\Sigma_d^{-1}(y-x)- \frac{1}{2}(x-\boldsymbol\mu_x)'\boldsymbol\Sigma_x^{-1}(x-\boldsymbol\mu_x) \]

where \(C\) consists of all the bits that don't depend on \(x\).

\[C = (2\pi)^{-\frac{k}{2}}|\boldsymbol\Sigma_d|^{-\frac{1}{2}} (2\pi)^{-\frac{k}{2}}|\boldsymbol\Sigma_x|^{-\frac{1}{2}}\]

Now, to maximise \(\log(p(x|y))\), we differentiate it. This involves some vector calculus, the rule we use can be seen in The Matrix Cookbook (page 11, eqn 85 and 86).

\[ \dfrac{d\log(p(x | y))}{dx} = -\boldsymbol\Sigma_d^{-1}(y-x) + \boldsymbol\Sigma_x^{-1}(x-\boldsymbol\mu_x) \]

To find the maximum we set the left hand side to zero, then solve for \(x\):

\[ 0 = -\boldsymbol\Sigma_d^{-1}(y-x) + \boldsymbol\Sigma_x^{-1}(x-\boldsymbol\mu_x) \]

With a tiny bit of algebra, we can manipulate the above equation to isolate \(x\),

\[\hat{x'} = \boldsymbol\Sigma_x(\boldsymbol\Sigma_x + \boldsymbol\Sigma_d)^{-1}(y - \boldsymbol\mu_x) + \boldsymbol\mu_x \]

and that is it. To get an estimate for the clean signal, we just use the equation above which takes the noisy signal. This is optimal (i.e. the best that is theoretically possible) if \(x\) and \(d\) are both Gaussian. It can still be O.K. if they are not, but you can easily repeat the above steps for different distributions. If the distributions aren't symmetric, then the MMSE estimator will usually perform better than the MAP estimator.

Expected Value of a Rayleigh Distribution

Computing the Mean

Rayleigh distributions pop up here and there, but the procedure outlined here for computing the expected value (i.e. the mean) of a distribution can be applied to any probability distribution function, it is just solving an integral. Sometimes these integrals can be pretty tough though, so it is best to use e.g. mathematica to solve them for you.

The Rayleigh Distribution is defined as follows:

\[f(x|\sigma) = \dfrac{x}{\sigma^2} \exp \left( -\dfrac{x^2}{2\sigma^2} \right) \]

The expected value of a distribution \(f(x)\) is defined as follows:

\[ \mathbb{E}[x] = \int x f(x) dx \]

By substituting our desired pdf into the above equation, and solving the integral we get the mean of the distribution:

\[ \mathbb{E}[x] = \int_0^\infty x \dfrac{x}{\sigma^2} \exp \left( -\dfrac{x^2}{2\sigma^2}\right) dx \]

We can now use the Wolfram integral calculator to solve it for us:

\[ \mathbb{E}[x] = \sigma \sqrt{\dfrac{\pi}{2}} \]

We can compare this to the value given on the Rayleigh Distribution page, so we know we got the right answer. Remember, this can be used for any distribution, not just the Rayleigh distribution.

Computing the Variance

Once we have the mean, we can also compute the variance using the following formula:

\[Var(x) = \mathbb{E}[x^2] - \mathbb{E}[x]^2 \]

Substituting in the value of \(\mathbb{E}[x]\) that we have already computed, we have the following equation:

\[Var(x) = \int_0^\infty x^2 \dfrac{x}{\sigma^2} \exp \left( -\dfrac{x^2}{2\sigma^2}\right) dx - \left( \sigma \sqrt{\dfrac{\pi}{2}} \right)^2 \]

Which simplifies very nicely after solving the integral with wolfram:

\[Var(x) = \sigma^2(2-\pi/2)\]

MMSE Estimators for Speech Enhancement

MMSE (Minimum Mean Square Error) estimators are commonly used in speech enhancement, with many different estimators applied to different speech representations. In its most basic form, an MMSE estimator is used to find the most likely clean signal given a noisy signal. I will start by explaining a time domain MMSE estimator (which won't be very useful), and then we can move on to better ones.

Some notation: \(x_n\) is our clean signal (\(n\) is the sample index), \(y_n\) is our noisy signal, and \(d_n\) will be our noise. We'll assume our noise is additive i.e. \[ y_n = x_n + d_n \]

We will now start applying some assumptions: we assume both \(x_n\) and \(d_n\) are gaussian variables, with different variances. As an assumption, this is not a very good one. Time domain speech is not gaussian, and most real noise is also not gaussian. One of the things this assumption has going for it though is that it results in a very simple estimator.

So I should now explain what an estimator actually is. Given a probability distribution function for a random variable, we want to pick the most likely value the random variable can take. There are two main ways of doing this: MMSE estimators, and MAP estimators. I'll explain MAP estimators in another post, they tend to be simpler to compute but they don't usually perform as well as MMSE estimators.

MMSE estimators are derived as follows: given a pdf \(p(x_n|y_n)\) (this is the probability of the clean signal given the noisy signal) we wish to find \(\hat{x}_n\) (this is called "x hat"),which is our estimate of the clean signal. To do this we find the expected value of the pdf \(p(x_n|y_n)\). We use the following formula for expectation:

\[\hat{x}_n = \int_{-\infty}^{\infty}x_n p(x_n|y_n)dx_n\]

Lets look at the pdfs we know (remember we assumed both \(x\) and \(d\) are gaussian distributed):

\[p(x_n) = \dfrac{1}{\sqrt{2\pi\sigma_x^2}} \exp \left( \dfrac{-x_n^2}{2\sigma_x^2} \right) \] \[p(d_n) = \dfrac{1}{\sqrt{2\pi\sigma_d^2}} \exp \left( \dfrac{-d_n^2}{2\sigma_d^2} \right) \]

Using the 2 pdfs above, we can derive the pdf for \(p(y_n|x_n)\) (the pdf of the noisy signal given the clean signal). Since we assumed \(y_n = x_n + d_n\), if we are given the value \(x_n\), we know \(y_n\) must be close to it with a variance of \(\sigma_d^2\) i.e. the noise variance. This results in:

\[p(y_n|x_n) = \dfrac{1}{\sqrt{2\pi\sigma_d^2}} \exp \left(\dfrac{ -(y_n-x_n)^2 }{ 2\sigma_d^2 } \right) \]

Despite how easy it is to calculate \(p(y_n|x_n)\), it is not the pdf we want. We can compute \(p(x_n|y_n)\) using Bayes rule:

\[p(x_n|y_n) = \dfrac{p(y_n|x_n)p(x_n)}{\int_{-\infty}^\infty p(y_n|x)p(x) dx}\]

The value we actually want is:

\[\hat{x}_n = \int x_n p(x_n|y_n) dx_n = \int \dfrac{x_n p(y_n|x_n)p(x_n)}{\int_{-\infty}^\infty p(y_n|x)p(x) dx} dx_n\]

Now we substitute the expressions for the pdfs into our equations to get a bit of a monster:

\[\hat{x}_n = \int{ \dfrac{x_n \dfrac{1}{\sqrt{2\pi\sigma_d^2}} \exp \left(\dfrac{ -(y_n-x_n)^2 }{ 2\sigma_d^2 } \right) \dfrac{1}{\sqrt{2\pi\sigma_x^2}} \exp \left( \dfrac{-x_n^2}{2\sigma_x^2} \right)}{\int_{-\infty}^\infty \dfrac{1}{\sqrt{2\pi\sigma_d^2}} \exp \left(\dfrac{ -(y_n-x)^2 }{ 2\sigma_d^2 } \right) \dfrac{1}{\sqrt{2\pi\sigma_x^2}} \exp \left( \dfrac{-x^2}{2\sigma_x^2} \right) dx} dx_n}\]

If you want to do it by hand you can simplify the above expression a bit and use a book of integrals to solve it, but the easier way is just to use mathematica to do it for you (or this page which helps with integrals). The result is the following very elegant solution:

\[\hat{x}_n = \dfrac{\sigma_x^2}{\sigma_x^2 + \sigma_d^2}y_n \]

If \(x_n\) and \(d_n\) were actually gaussian random variables, then the estimator we derived above is optimal i.e. no better estimator of \(y_n\) exists. In reality, we would also need to know \(\sigma_d^2\) and \(\sigma_x^2\), and these can be very tricky to estimate. In the next post I will use the same method but different probability distributions to derive other estimators.

In the next post I'll do a slightly better MMSE estimator that doesn't use gaussians. By using more complicated pdfs we can improve our estimator, but we will also get much more complicated integrals. In some cases the integrals won't be analytically solvable at all, so we'll have to be happy approximating them numerically. If this is not possible because of e.g. timing constraints, we'll have to settle on a MAP estimator.