5.1 Ideal Filters and the Ringing Problem

The Ringing Problem. We have seen that the inverse discrete Fourier transform (IDFT) of an ideal lowpass filter with cutoff frequency d₀ and total pass bandwidth is the the sinc function

h(d) = Asinc(2d₀ d) = Asinc(2d₀ d)/(2d₀ d) (5.1)

where d is the distance to the origin in the spatial domain. The DFT-IDFT pair is

h(d) = Asinc(2d₀ d) [A/(2d₀ )]R(d_f /(2d₀ )) = H(d_f ) (5.2)

where the rectangular function is defined by R(d_f /(2d₀)) = 1 for |d_f | d₀, R(d_f /(2d₀)) = 0 for |d_f | > d₀, and d_f is the distance to the origin in the frequency domain.

Figure 5.1 shows the right half of the sinc function h(d_f). This function continues to oscillate with both positive and negative values as |d| . In order to use it as a convolution mask of finite size, we must truncate it to the function h^~(d) that is equal to 0 for d larger than some reasonably small value. But when we truncate h(d) to h^~(d), we change the filter H(d) to H^~(d). Figure 5.2 presents the new filter H^~(d) due to the truncation, which oscillates in a fashion known as the Gibbs phenomenon. The oscillation of the magnitude at the various frequencies is called ringing and it has undesirable effects on the filtered G(u,v) and therefore on the processed image g(m,n).

Figure 5.1. A sinc h(d) for convolution. Figure 5.2. The Gibbs phenomenon.

Upon putting A = 2d₀ in Equation (5.1), h(d) is a sinc function with height 2d₀ such that

(A/2d₀ )R(d_f /(2d₀)) = 1, |d_f | d₀ (5.3)

The positive half of the corresponding impulse response function is shown in Figure 5.1, where the discrete variable m is used in place of d. The first zeros of sinc(m) are at m = ±1/(2d₀). The Fourier transform pair therefore becomes

(2d₀ )sinc(2d₀ d) = (2d₀ )sin(2d₀ d)/(2d₀ d) H(d_f /2d₀ ) = 1, |d_f | d₀ ; else 0 (5.4)

Lowpass Sinc Convolution Masks. To obtain a convolution mask for the 2-dimensional lowpass filter for lowpass filtering in the spatial domain, we obtain the impulse response function of spatial distance d to be

h(d) = (2d₀ )sin(2d₀ d)/(2d₀ d) (5.5)

We obtain a 3x3 convolution mask from this as the array

h(d₂) h(d₁) h(d₂)

h(d₁) h(d₀) h(d₁) (5.6)

h(d₂) h(d₁) h(d₂)

We could use the city-block distance (norm), which sums the number of unit steps in the x and y directions, so that d₁ = 1, d₂ = 2. If we use the Euclidean norm, then d₁ = 1 and d₂ = 2. The maximum norm [6] is our choice here, for which d₁ = 1 and d₂ = 1 (the definition is (x,y) - (x₀,y₀)_max = max{|x-x₀|, |y-y₀|}).

The maximum norm allows all mask entries adjacent via a side or corner to the center value to have equal values, so that in the processing, the pixels adjacent to the center of a pixel neighborhood will have the same weight. The pxp mask determined by the sinc function has concentric squares of equal values around the center (the sinc function is centered on the center of the mask to obtain the values h(d_k)). Upon standardizing d to d/D, for D = max{M-1,N-1}, Equation (9) becomes

h(d) = (2d₀/D)sin(2d₀d/D²)/(2d₀d/D²) (5.7)

For example, if we have a 256x256 image and use the maximum norm, a 3x3 mask, and design a cutoff frequency of d₀ = 224 (to filter off the highest frequencies as noise), then

h(d) = (2x224/255)sin(2d₀ d/255)/(2d₀ d/255) 1.76sin(1.76d)/(1.76d) (5.8)

Because d = 0 or d = 1 for all entries in the 3x3 mask with maximum norm, we obtain the moderate sharpening mask (suggestive of a Laplacian mask [2])

h(1) h(1) h(1) -0.225079 -0.225079 -0.225079

(r) h(1) h(0) h(1) = (r) -0.225079 1.750000 -0.225079 (5.9)

h(1) h(1) h(1) -0.225079 -0.225079 -0.225079

For a cutoff frequency of d₀ = 64, we have the lowpass impulse response function

h(d) = (2x64/255)sin((2x64/255)d)/((2x64/255)d) = 0.5sin(d/2)/(d/2) (5.10)

that has the averaging mask form

h(1) h(1) h(1) 0.318310 0.318310 0.318310

h(1) h(0) h(1) = 0.318310 0.500000 0.318310 (5.11)

h(1) h(1) h(1) 0.318310 0.318310 0.318310

This is a smoothing (averaging) filter. It passes only the lower 1/4 of the frequencies and so we expect that it blurs the image.

Note that the entries do not sum to unity, so the processed pixel values will be biased toward the darker or brighter end of the gray scale. To prevent this, we sum up the entries h(d_k) to obtain a standardizing divisor S and the divide all entries by S, i.e., the mask we actually use would be {h(d_k)/S}. We do the above process for all pxq masks (which may be 5x5, 7x7, 9x9 or larger, or 5x7, etc.). Because the mask values may oscillate in sign as a function of distance from the center, we sum the positive and negative mask values separately and add a positive number to the central value, if necessary, to keep the total sum positive and approximately unity.

Figure 5.3 A deleted sinc function

Highpass Sinc Convolution Masks. Let H_low(u,v) be the ideal lowpass filter in the frequency domain with IDFT h(d). Then H_high(u,v) 1 - H_low(u,v) h_high(d) = (d) - Asinc(2d₀d) = (d) - 2d₀sinc(2d₀d) is the transform pair, where (d) is the delta impulse (infinite spike) function that is the inverse Fourier transform of the constant unity filter R(d_f) = 1, |d_f| 0. We call h_high(d) the deleted sinc function, shown in Figure 5. A similar function appears in [1], although sincs are somewhat discounted there because of ringing. Using a boost gain of ß (ß > 1) and an attenuation of (0 < < 1), we design a highpass filter that partially passes low frequencies, per

h_high(d) = ß(d) - sinc(d) = ß(d) -

(2d₀/D)sin(2d₀d/D)/(2d₀d/D) (5.12)

We sometimes use ß = 1.4 and = 0.4.

Notch and Bandpass Convolution Masks. Notch and bandpass filters are opposites in that

H_notch(d_f) = 1 - H_band(d_f) (5.13)

The notch filter can be constructed by summing a lowpass and a highpass filter that do not overlap, or

H_notch(d_f) = H_low1(d) + H_high2(d) (5.14)

for d_0,low1 < d_0,high2 (lowpass cutoff is less than the highpass cutoff).

This can also be put in the form

H_notch(d_f) = H_low1(d_f) + {1 - H_high2(d_f)} = 1 + H_low1(d_f) - H_high2(d_f) (5.15)

h_notch(d_f) = ß(d) + ₁(2d_0,low1/D)sin(2d_0,low1d/D)/(2d_0,low1d/D) -

₂(2d_0,high2/D)sin(2d_0,high2d/D)/(2d_0,high2d/D) (5.16)

A bandpass filter and its impulse response can be constructed via

H_band(d_f) = 1 - H_notch(d_f) = H_low2(d_f) - H_low1(d_f) (5.17)

h_band(d) = (2d_0,low2/D)sin(2d_0,low2d/D)/(2d_0,low2d/D) -

(2d_0,low1/D)sin(2d_0,low1d/D)/(2d_0,low1d/D) (5.18)

where d_0,low1 and d_0,low2 are respective cutoffs for lowpass and highpass filters.

5.2 Gaussian Windowed Sinc Processing

Figure 4 shows the original "Shuttle" image of 320x240 pixels taken from a photograph. It can be seen that although it is a good picture, it is not highly focused. Figure 5 shows the focusing effect of a wideband lowpass filter.

A narrow (wide) Gaussian maps via the Fourier transform to a wide (narrow) Gaussian. Gaussians have a gentle roll-off. With a Gaussian window (envelope) imposed on the sinc, the lowpass filter is the convolution of a pulse with a Gaussian in the frequency domain. The resulting filter has the shape of a wide bell shaped curve. In the spatial domain, the window has the form

w(d) = exp[-d²/(2²)] (5.19)

Figure 6 shows the result of processing the original image with a lowpass sinc in a Gaussian window with mask size 3x3, d₀ = 80 and = 1.0. The resulting image is blurred. Figure 7 is the result of a lowpass 5x5 mask, d₀ = 300 and = 2.0 for the Gaussian window. The processed image is very sharp. Figure 8 used a lowpass 9x9 mask, d₀ = 300 and = 4.0. It is rather oversharpened. Figure 9 is the result of highpass Gaussian windowed sinc convolution with 3x3 mask, d₀ = 100 and = 1.0. Thus the passband is 100 to 320. It is lightly sharpened.

5.3 Gaussian Windowed Sincs

We sometimes need a window that has a more rapid roll-off to allow us to use smaller convolution masks without ringing. We could use Hamming or Hanning raised cosine windows [6, p. 436], but those are actually bell shaped curves very much like Gaussians. We choose the square of the Gaussian because it has an increasingly more rapid roll-off farther from the center. The effect is a more slender bell shaped curve that goes to zero more quickly. The filter in the frequency domain that corresponds to this is a convolution of a Gaussian with the bell shaped filter (the convolution of a Gaussian with the convolution of a Gaussian and a pulse).

Figure 10 presents the results of filtering the original image with a 5x5 squared Gaussian sinc mask, d₀ = 100 and = 2.0. Figure 11 shows the results when Figure 10 is passed a second time through the same filter. These results are rather striking.

Fig. 4. Original shuttle image. Fig. 5. 3x3 wideband lowpass, d₀ = 200.

Fig. 6. Gaussian sinc, 3x3, d₀=80. Fig. 7. Gaussian sinc, 5x5, d₀=300.

Fig. 8. Gaussian sinc, 9x9, d₀=100. Fig. 9. Gaussian sinc, 3x3, d₀=100, =1.0.

Fig. 10. Squared Gaussian windowed sinc, Fig. 11. Double Squared Gaussian windowed

5x5, d₀=100, =2.0. sinc, 5x5, d₀=100, =2.0.

5.4 Gaussian Filters

Gaussian filters can be implemented in the spatial domain because the IFT of a Gaussian filter in the frequency domain is a Gaussian convolution function in the spatial domain whose spread parameter (standard deviation ) is essentially the reciprocal of the spread parameter of the frequency filter. Gaussians are gentle curves with no sharp corners so that ringining/aliasing is not a problem. To form a lowpass Gaussian filter, we put a wide Gaussian centered on the origin in the frequency domain {(u,v)}. The spatial domain filter that corresponds to it is the narrow Gaussian convolution function that determines a convolution mask that is centered on each pixel in turn to convolve to obtain the new pixel value.

A high pass filter in the frequency domain is the constant unity minus the lowpass Gaussian. To boost the high frequencies, we usually take 2 times the constant to obtain the high pass filter

2 - exp[- d² / (2² )] (5.20)

The convolution function that results from taking the IFT of this is

2 - exp[-(d - d₀ )² / (2² )] (5.21)

5.5 Homomorphic Filtering

The homomorphic filter moderately attenuates the power level at the lower frequencies while it also moderately boosts the power at the high frequencies. The overall effect is the strong enhancement of an image. Here we demonstrate that this type of filter can also be fashioned from a constant and a Gaussian. Consider Figure 5.12. The Gaussian is subtracted from the larger one to yield a filter that suppresses the lower frequencies while boosting the higher ones.

The is accomplished by taking the constant > 1 and multiplying Gaussian by a constant < 1 to obtain

H(d) = - exp[-d)² / (2² )] (5.22)

Figure 5.12. A homomorphic filter.

5.7 Computer Programs

5.9 Exercises