5.2 Deconvolution with inverse and weiner filters

Problem statement

The first task is to deconvolve the reference pulse from all the measured waveforms to obtain the actual projections of our specimen. In general the deconvolution problem can be modeled as in the following figure.

In the above figure, p(t) is the projection of the object slice that is to be used in the object reconstruction. The reference pulse is denoted by r(t) and interpreted as a low- pass degradation filter that convolves with the projections. The result of this convolution plus an additive noise n(t) is the measured waveforms s(t). In this case, both r(t) and s(t) are known and the goal is to obtain p(t) by deconvolving r(t) from s(t) and denoising. Mathematically the above block diagram can be written as:

Inverse filter

Most of deconvolution schemes are implemented in the frequency domain. This makes sense since convolution in the time domain is mapped as multiplication in the frequency domain. Inverse filtering is the simplest and most naïve method for deconvolution. No denoising takes place in this case.

Equation 1 can be written in the frequency domain as:

The main drawback of the method is that R(ω) is usually a low-pass filter and therefore 1/R(ω) is a high-pass filter which attains large values as the frequency increases. Thus, equation (3) becomes numerically unstable for small R(ω) values and greatly amplifies the high- frequency noise contribution. This makes inverse filtering very sensitive to even the small amounts of high frequency noise which exists in the measured waveforms.

One method of decreasing the noise sensitivity inherent when inverse filtering is to bound the frequency response 1/R(ω) to some threshold γ as follows:

The code that implements inverse filtering is inverseFilter.m

Wiener filter

One of the most widely used restoration techniques is the Wiener filter . Contrary to the inverse filtering this method also attempts to diminish noise while restoring the original signal. It executes an optimal balance between inverse filtering and noise smoothing in the mean square error sense. Assuming white Gaussian noise, the Wiener filter expressed in the Fourier domain is written as:

In this case the projections are unknown and an estimate of Spp is needed. This is provided by noting that Sss(ω)= Spp(ω)•|R(ω)|2, solving for Spp(ω) and substituting into equations (5) to get:

Interestingly, the last equation can be interpreted as two different filters in cascade in the frequency domain. The first is denoising, while the second is exactly the inverse filter considered in the previous section.

The code that implements Wiener filtering is wienerFilter.m .

Driver for filters