The Wiener-like filtering in the wavelet space

Let us consider a measured wavelet coefficient w_i at the scale i. We assume that its value, at a given scale and a given position, results from a noisy process, with a Gaussian distribution with a mathematical expectation W_i, and a standard deviation B_i:

$$P(w_i/W_i) = \frac{1}{\sqrt{2\pi}B_i} e^{- \frac{(w_i-W_i)^2} {2B_i^2}}$$ (14.76)

Now, we assume that the set of expected coefficients W_i for a given scale also follows a Gaussian distribution, with a null mean and a standard deviation S_i:

$$P(W_i) = \frac{1}{\sqrt{2\pi}S_i}e^{-\frac{W_i^2}{2S_i^2}}$$ (14.77)

The null mean value results from the wavelet property:

$$\int_{-\infty}^{+\infty} \psi^*(x) dx = 0$$ (14.78)

We want to get an estimate of W_i knowing w_i. Bayes' theorem gives:

 

$$P(W_i/w_i) = \frac{P(W_i)P(w_i/W_i)}{P(w_i)}$$ (14.79)

We get:

$$P(W_i/w_i) = \frac{1}{\sqrt{2\pi}\beta_i}e^{-\frac{(W_i-\alpha_i w_i)^2}{2\beta_i^2}}$$ (14.80)

where:

$$\alpha_i = \frac{S_i^2}{S_i^2+B_i^2}$$ (14.81)

the probability P(W_i/w_i) follows a Gaussian distribution with a mean:

$$m = \alpha_i w_i$$ (14.82)

and a variance:

$$\beta_i^2 = \frac{S_i^2B_i^2}{S_i^2 + B_i^2}$$ (14.83)

The mathematical expectation of W_i is $\alpha_i w_i$.

With a simple multiplication of the coefficients by the constant $\alpha_i$, we get a linear filter. The algorithm is:

1.

Compute the wavelet transform of the data. We get w_i.

2.

Estimate the standard deviation of the noise B₀ of the first plane from the histogram of w₀. As we process oversampled images, the values of the wavelet image corresponding to the first scale (w₀) are due mainly to the noise. The histogram shows a Gaussian peak around 0. We compute the standard deviation of this Gaussian function, with a $3\sigma$ clipping, rejecting pixels where the signal could be significant;

3.

Set i to 0.

4.

Estimate the standard deviation of the noise B_i from B₀. This is done from the study of the variation of the noise between two scales, with an hypothesis of a white gaussian noise;

5.

S_i² = s_i² - B_i² where s_i² is the variance of w_i.

6.

$\alpha_i = \frac{S_i^2}{S_i^2+B_i^2}$.

7.

$W_i = \alpha_i w_i$.

8.

i = i + 1 and go to 4.

9.

Reconstruct the picture from W_i.