Exploration of Digital Filter Dynamics

1. Introduction

This discussion will explore advanced signal processing techniques, and extend some of those techniques to two-dimensional image-processing applications. Many of examples presented can be reproduced in MATLAB, including two-dimensional Fourier transforms, deconvolution, and bilateral filtering.

2. Linear Filters

Given a signal x, and some impulse response h, a new signal can be defined as a moving average (MA) by a summation.

y_t = S_k h_t-k x_k

The summation symbol, S, will use an index k covering all integers as a summation index. This impulse response may have infinitely many elements in the case of IIR filters. Although in practice, filters will work on time-limited signals using a time-limited response, so that x and h arrays in a computational model will be finite in extent. When implementing stable filters, values may be taken to be zero when results are within an acceptable range. The resulting sequence is known as the convolution of x with h, which can be written y = h*x.

These finite impulse response (FIR) filters are linear, and behave in accordance with the principles of superposition. This allows the effects of individual filters to be combined easily by convolution. The effects of linear filters are independent of the particular order in which two filters are applied, and that, under certain conditions, a stable inverse filter can often undo the filtering.

3. Applying linear filters to two dimensions

In two dimensions, a signal and the impulse response, in the linear shift invariant sense, can be extended into the two-dimensional spatial domain by modifying the definition of convolution to include two-dimensional arrays.

y(n₁,n₂) = S_k1 S_k2 h(n₁-k₁,n₂-k₂) x(k₁,k₂)

In some cases, a filter may be factored into two filters, each dependent on only one spatial dimension. Such a filter is called separable. Separability permits techniques from the 1-D domain, such as factoring into couplets, to applied quite readily. Separable filters may be written in a fashion which shows these constraints on functional dependence.

y(n₁,n₂) = S_k1 S_k2 h₁(n₁-k₁) h₂(n₂-k₂) x(k₁,k₂)

In addition, some impulse responses may define rotationally-symmetric filters. Particularly for image-processing, such filters can serve to minimize effects caused by the coordinate system in which graphics pixels are arranged, or direction in which a signal is oriented. For instance, a linear smoothing filter might act independently of relative angle, and be dependent solely on the square of the Euclidean distance:

h(n₁-k₁,n₂-k₂) = h_r ((n₁-k₁)² +(n₁-k₁)²)

As an example, one might wish to design a zero-phase filter based on the function h(n₁,n₂) = 1 /(1 + n₁² + n₂²), where the denominator of a pixel weight is one larger than the square of the distance from the pixel being evaluated. In the Fourier analysis, this causes frequency response to fall off without strong dependence on direction of sinusoidal input. When graphed, the lines of equal response are directed along concentric circles. To keep the response finite, we could also reject any pixels more distant than the square root of 5, yielding an approximately rotationally-symmetric smoothing filter:

Impulse Response Values h(n1,n2)					Frequency Response |H( w1,w2)|
0	1/6	1/5	1/6	0
1/6	1/3	1/2	1/3	1/6
1/5	1/2	1	1/2	1/5
1/6	1/3	1/2	1/3	1/6
0	1/6	1/5	1/6	0

Although not shown in the graph, the frequency response will climb back to high levels as frequencies exceed the Nyquist frequency. When evaluating the frequency response, each frequency component may independently approach multiples of 2p when considering aliasing effects. Therefore, components are necessarily aliased with strong dependence on angle of the wavefront, not solely on spatial frequency of a sinusoid. For instance, along a diagonal, aliased frequencies are spaced farther apart, having frequency response periodicity 2p multiplied by the square root of 2. In any general sense, we can never truly escape functional dependence on direction. Therefore, the notion of rotationally symmetric filtering is only valid at reasonably low spatial frequencies where aliasing is not encountered.

The Z transform and the Fourier transform may be extended into two dimensions.

Two-dimensional Z transform:     X(Z₁,Z₂))= S_n1 S_n2 x(n₁,n₂) Z₁^n₁ Z₂^n₂

Fourier transform:     X(w₁,w₂) = S_n1 S_n2 x(n₁,n₂)e ^{-iw₁n₁ -iw₂n₂}

In the case where x is a sampled sinusoidal signal directed at some angle q, such as the signal x(n₁,n₂) = e ^{iw₁n₁ +iw₂n₂}, we see that component frequencies w₁ and w₂ are dependent on both the angle q and on the spatial frequency w. When evaluating frequency response of the transfer function using the Z transform, we can use the component frequencies on the unit circle, or their trigonometric representation in terms of w and q for each dimension.

Z₁ = e ^-iw₁ = e ^{-iw cos q}

Z₂ = e ^-iw₂ = e ^{-iw sin q}

This technique permits frequency response to be expressed in polar coordinates as H(w, q), in addition to the usual component expression H(w₁, w₂).

In examining the dynamics of certain filters, it is worthwhile to consider what ordinarily happens to sharp signals as high-frequency components are suppressed by symmetric averaging, or equivalently multiplication in the frequency domain. The effects of limiting the spectral bandwidth can be generalized to the spatial domain, and interpreted as a smoothing of the sharp edges often present in photographic images.

In the following graph, a boxcar function, or a short pulse of ones, is shown. Its Fourier transform is related to the sinc function, and the magnitude spectrum is shown with sinc-like ripples. Using only the low frequency terms from the Fourier transform, an approximation of the function can be generated. The result is roughly similar to the concept of linear smoothing in the time domain, as high frequency signal components are suppressed. Approximations of the boxcar function are shown with 2, 4, and 8-term approximations.

Similarly, a Gaussian function can be examined in this manner. One notable difference is that the magnitude spectrum is also a Gaussian function.

In either case, there is a continuing trade-off between signal preservation, and attenuation of potential high frequency noise. The area under the curves is constant, so each approximation is partially a matter of how "spread out" the pulse has become as it lacks components of higher frequency.

In two-dimensions, we might consider the boxcar function rotated about a central axis. The resulting signal appears as a disc, with a sharp circular boundary suggesting sinc-like behavior in the frequency domain. The importance of examining the magnitude with a logarithmic scale is apparent in this example, where variations in high-frequency response become nearly invisible on the linear scale. In the spatial domain, this circular window might represent a telescope mirror, which only admits light striking the disc-shaped area.

Circular Window

Magnitude Spectrum (linear scale)

Magnitude Spectrum (log scale)

Some of these image-processing techniques can be generalized to three-dimensions, allowing a sequence of images to be processed using information about the optics. For instance, images might be sampled at varying focal planes of a microscope slide. The following images were collected from onion cells, stained with methylene blue.

Onion cells, image A:

Onion cells, image B:

After deconvolution, some the blur imaged from adjacent focal planes is removed, although efficient use of this process requires knowledge of the point spread function (psf) of the optics involved.

Processed Image, image A:

Processed Image, image B:

Currently, this technique is becoming more widespread in biological microscopes, which can also use deconvolution to visualize cells in three-dimensions, in addition to ordinary two-dimensional images.

4. Nonlinear Filtering - the Median Filter

As an alternative to using linear filters, one commonly employed technique is median filtering. This operation takes a given signal and evaluates a local median at each point using neighboring values as a basis. This nonlinear technique can remove single-point noise caused by experimental errors, or other sampling errors that occur at single points due to bad pixels or other systemic causes. In image-processing, such noise is often called salt-and-pepper noise, which describes the granular appearance of individual points of noise in the signal.

Landscape Image with Noise:

Median Filtered:

A median filter can effectively suppress this type of noise. But this type of filtering is not ideal for suppressing white noise, nor is it reversible. One advantage of the median over linear filters is the elimination of isolated points of noise.

5. Bilateral Filtering

A recent development is the bilateral filter, a technique proposed by Tomasi and Manduchi in 1998. This technique preserves edges by mixing a moving average technique with a nonlinear system of weights. Loosely speaking, it relies on an assumption that any noise is more uniformly distributed, and that signals have distinct edges or steps that can be detected by examining local differences. As with a linear smoothing filter, each neighboring value is weighted on its proximity in space or time (a domain weight). But then, a second weighting factor gives some measurement of local difference (a range weight). Often, these two weights are expressed as a pair of Gaussian distributions in a summation. In the time domain, this is expressed as:

S_k x_k W_D W_R
where W_D=e^{-(t-k)2 / s_D2} and W_R=e^{-(x_t-x_k)2 / s_R2}

As with linear smoothing filters, the effect of more distant values will be diminished, resulting in stability and a low-pass type of behavior. A mean value on the time index, s_D, determines the level of attenuation of higher frequencies. But also, another mean value, s_R, weighs the effect of relative radiosity, or, most specific to time-domain data, a local change amplitude. Large changes in amplitude will be assumed to be signal-driven, and given less weight in the smoothing sum. Such a measurement could be interpreted as allowing signals, which sometimes have discrete steps or edges, to be less attenuated than the more randomly occurring noise at corresponding frequencies.

Usually the sum must be divided by a normalizing factor. In the following example, a sampled 30 Hz sinusoidal function is modulated with 7 Hz square wave to create pulses, or abrupt changes in the signal, over a one second interval at 1 kHz. The result is then combined with white noise.

A 5-point linear filter h=(1 2 4 2 1) is applied in an attempt to smooth the signal, along with a 5-point median and a 5-point bilateral filter. The choice of mean values for the bilateral filter were s_D =1.5 ms, with a choice of s_R at half-amplitude. A sum-squared computation of error is performed to compare with the original modulated sinusoid before any noise was added.

In the linear filter, abrupt steps in the signal are attenuated along with the noise, creating clipped peaks and reduced sharpness of any discontinuous steps. The median filter also clips the signal in some areas, and does not suppress the white noise as well as linear filtering. The bilateral filter partially avoids these effects, with smoothing in the continuous parts of the signal being close to smoothing of the linear filter. Looking at the magnitude spectrum on logarithmic scale, we can view each signal using a Hamming window, and consider the magnitude spectra:

The linear filter more effectively suppresses high frequencies, which is evident by a decline in the magnitude spectrum. Although, when comparing with the spectrum of the original signal prior to any noise being added, the linear filter also suppresses valid signal components as well.

When applied to two-dimensions, including images with color components, the preservation of edges is evident. Since signal components with distinct edges are often visualized as sinc functions in the frequency domain, it can be expected that high-frequency components approaching the Nyquist frequency, will be preserved in the signal. However, the same frequencies will be attenuated in the noise.

Landscape Image with Noise:

Bilaterally Filtered Image:

6. Cepstrum Analysis

In the magnitude spectrum of the one-dimensional example, the original signal shows some degree of regularity in the magnitude response. The linear filter diminishes the high frequencies, but does not restore the shape of the spectrum as a whole. A technique which analyzes the spectrum is called the cepstrum (pronounced "kepstrum"), which derives its name from the word spectrum. This technique takes the logarithm of the Fourier transform, with or without the complex phase information, and then evaluates the Fourier transform of that result.

cepstrum(x) = IDFT(LOG(DFT(x)))

When using complex values, careful application of the logarithm must be applied to unwrap phase angles. Signal echoes or other repetitions of signal can be detected. This technique has applications in sonar, speech processing, and other fields. Although the cepstrum transforms the time domain data back into the time domain, we are now looking at a very different signal whose seconds are not points in time, but some measure of spectral periodicity. Looking at the cepstrum of the previous examples, we can see signal patterns that the bilinear filter partially recovers.

In two-dimensions, this technique can sometimes detect blurring, which suggests a linear deconvolution process for image restoration. An image of a wind generator is shown, with a blurring detected by 2D-cepstrum using a Gaussian window, and subsequent linear operations are used to correct the blur.

Blurred Image:

Unblurred by Deconvolution:

7. Conclusions

Each of the filters considers has characteristics which can be analyzed and considered when selecting filters for a particular application. The following summarizes some of this information.

Filter characteristics:

Filter	1-D Formula	Linear	Shift Invariant	Preserves Edges	Suppresses Single Noise	Inverse Possible
Moving average	yt = Sk hk-txk	Yes	Yes	No	No	Yes
Median Filter	y = median(xt-k,...,xt+k)	No	Yes	Yes	Yes	No
Bilateral Filter	yt = Sk xk WDWR / Sk WDWR	No	Yes	Yes	No	No

More work might be done in discovering cepstral properties of the sinc function, and how range weights are mathematically suited for enhancing cepstral components in edge-preserving smoothing filters.

References

[1] J. Karl, An Introduction to Digital Signal Processing, Academic Press 1989.

[2] J. Kovacevic, z-Transform, AT&T Bell Laboratories 1996.

[3] M. Elad, Analysis of the Bilateral Filter, Stanford University.

[4] W. Jiang, M. Baker, Q. Wu, C. Baja, and W. Chiu, Applications of a bilateral denoising filter in biological electron microscopy, Journal of Structural Biology, 2003.

[5] J. Long, Digital Image Processing: Homomorphic Image Enhancement, Image Restoration, and Phase-Only Image Processing, University of Wisconsin Oshkosh, 1994.

[6] C. Tomasi and J. Manduchi, Bilateral Filtering for Gray and Color Images, Proceedings of the 1998 IEEE International Conference on Computer Vision, Bombay, India.