Seeing through Chaos: Noise Filtering

Mind is a powerful thing. It is both the cause of chaos in our lives and the source of our redemption from it.

 This chaos which religion describes as an "Act of God", is known to Geeks as "Noise". Obviously Rockstars are oblivious to the fact that they are actually pleading for "Chaos" when they ask their audience to "Make some Noise \m/".

Today let us attempt to see through this veil of randomness and uncover the Signals hidden behind. By the way, a very prosperous and a happy new year to all my readers, if that makes you feel any better :p. 

Noise removal using Wiener Filter

A noisy image is an eyesore, yuck! So do we fold our hands and stare at it helplessly. No, we protest....vehemently. We oppose it, suppress it and humiliate it....yeah!! You can't disagree once you've seen this:

Left: Original Image                                                    Right: Noisy Image 

The noise that you are seeing here is called 'Additive White Gaussian Noise'. It is additive because the noise signal (yeah! noise is also a signal, a pathetic one at that but a signal indeed) is algebraically added to the pixel intensities. It is 'Gaussian' because the magnitude of noise at any particular pixel is Gaussian distributed random variable. It is 'white' because just as white light is composed of all possible wavelengths of light, white noise contains all possible frequency components i.e. its Power Spectral Density (PSD) is flat. But what the hell is 'Power Spectral Density' right? Let us dig deeper here because wiener filter stands on the pedestal of PSD and its estimation for the signals in question is absolutely essential.

Power Spectral Density

The idea of PSD resonates with that of Fourier Transforms, in fact by definition for continuous time signals:

\begin{align}PSD(\omega)=E[|F(\omega)|^2] \end{align}

where \(F(\omega)\) is the Fourier Transform of signal \(f(t)\) and \(E[.]\) denotes the expected value. Basically as the name suggests it is a measure of the power contributed by the different frequency components to the original signal. The task of finding the PSD for any deterministic signal such as an image is simple. Because images are discrete signals (at least on computer where the computation takes place), we shall use Discrete Fourier Transform (DFT) instead of the continuous one. For deterministic images, such as the original cameraman image above:

\begin{align}PSD[m,n] = \frac{1}{MN}|F[m,n]|^2 \end{align}

Let me tell you a secret - Continuous-time Fourier Transform (CFT), Discrete-Time Fourier Transform (DTFT) and Discrete Fourier Transform (DFT) used to give me a headache until recently when the concepts began to gel in. So, if you are passing through the same bad weather, may God bless you with one of my future posts :p. Here's a quick hand-wavy recap:

When a continuous time function is sampled, then in frequency domain its transform replicates after intervals of \(\frac{2 \pi}{T}\). Further the transform undergoes both frequency scaling (by a factor of T) and amplitude scaling (by a factor of  1/T) where T is the time period. However for discrete computations in frequency domain in computers we need to sample even the DTFT and because of the periodicity property only the section between \([0, 2 \pi]\) is used. So in this interval for a 1-D signal, for example, the relation between DTFT (\(F(\omega)\)) and DFT \(F[k]\) is given by:

\begin{equation}F[k]=F(\omega)|_{\omega=\frac{2 \pi k}{N}}\end{equation}

As far as the noise is concerned, since we are dealing with AWGN, its PSD is simply the noise variance at all frequencies. What could be simpler than that right?

The Wiener Filter

The idea of Wiener Filter is intuitive. We want to suppress those DFT coefficients where the Signal-to-Noise ratio is small while leaving those coefficients untouched where the Signal-to-Noise ratio is large. Wiener Filter is actually an optimal filter in the sense that it minimizes the mean squared error between the noisy signal and the underlying signal. A rigorous derivation is uploaded here.  But the final result i.e. the filter in frequency domain is give by:

\begin{equation}H(\omega)=\frac{S_u}{S_u+S_\eta}\end{equation}

where \(S_u\) and \(S_\eta\) are the PSD of the original image and noise respectively. This was what the professor who taught me image processing, but who shall remain unnamed here, told me. And I was like: This dude must be crazy! How is it possible to find the PSD of the original image when we are given just the noisy image. But it turns out that a reasonable estimate can be obtained from the noisy image as well. Now I don't know if this is a standard practice but my approach was this:

The trick here is to first perform a Gaussian smoothing operation on the noisy image. This acts as a low pass filter and filters out the high frequency noise components. The PSD of the resulting image would then serve as a not-too-bad estimate of the original image. Of course, now there is a trade-off because Gaussian smoothing also smudges the image. More the smudge, more is the noise that would be filtered out and vice versa. Take your time to appreciate the results below and you can find the entire MATLAB code here.

Noisy image after Gaussian smoothing operation

           Left: Noisy Image                       Right: Wiener Filtered Image

Wiener Filtering using the PSD of original image

There was some gap between my last blog and this one. The reason being I promised you, my readers, that I would tell you how to filter out noise from a signal. I did not expect that I would run into practical issues which we seldom think about while trying to understand the theory. The point is I learnt a lot myself while preparing for this particular blog. Hope you enjoyed it and found it equally informative.