The Fourier Transform Part 1 (Biomedical Image Analysis)

French mathematician Joseph Fourier found that any 2^–periodic signal f (t) can be represented by an infinite sum of sine and cosine terms according to

 

tmp2F-89_thumb

When a signal f (t) is given, the coefficients ak and bk can be determined using Fourier analysis of the signal f (t):

tmp2F-90_thumb

where ak and bk are called the Fourier coefficients. The term 2^-periodic indicates that the signal repeats after integer multiples of 2^, that is, f (t) = f (t + 2^n), where n is an integer number. The periodicity is a consequence of the sine and cosine terms in Equation (3.1) and will become a fundamental assumption in the Fourier analysis of signals and images. The cosine and sine functions in Equation (3.1) form an orthogonal basis, since their inner product is zero:


tmp2F-91_thumb

Equation (3.3) needs to be evaluated for k = 0, 1,2, …, to to obtain the infinite number of coefficients that are needed to exactly represent the signal f(t). The special case k = 0 yields b0 = 0, and a0 represents the mean value of the signal over one period. In practice, a limited number of coefficients is sufficient to represent f(t) with satisfactory accuracy. The coefficients ak and bk identify the contribution of the kth harmonic oscillation to the overall signal. The kth harmonic signifies an oscillation with k times the frequency of the base oscillation (i.e., the oscillation with a period of 2^). For this reason, Equation (3.3) is also called harmonic analysis.

Consider the piecewise linear, periodic oscillations in Figure 3.1A (square wave) and Figure 3.1C (triangular wave). The triangular wave closely approximates a sine wave (gray curve in Figure 3.1A). The sine wave has only one nonzero coefficient, b1 = 1, as can be seen in Equation (3.1). Since both the square wave and the triangular wave can be described analytically, the integrals in Equation (3.3) can be solved to obtain the Fourier coefficients ak and bk. The average of the signal over one period is zero, and a0 = 0 follows immediately. In fact, all ak are zero because of the function’s symmetries. Furthermore, solving the second integral in Equation (3.3) yields bk = 0 for all even values of k and nonzero coefficients bk, as described for a square wave and a triangular wave, respectively:

tmp2F-93_thumb

Fourier analysis of two periodic functions. A square function (A) and a triangular function (C) are shown next to a sine wave (gray curve in A and C). Solving the integrals in Equation (3.3) yields ak = 0 for all k in both cases as a consequence of the symmetry. The bk are plotted in B for the square wave and in D for the triangular function. It can be seen that the Fourier coefficients of the triangular function drop off more rapidly for higher k than those of the square wave, because the triangular function is more closely related to the sine wave than is the square wave.

FIGURE 3.1 Fourier analysis of two periodic functions. A square function (A) and a triangular function (C) are shown next to a sine wave (gray curve in A and C). Solving the integrals in Equation (3.3) yields ak = 0 for all k in both cases as a consequence of the symmetry. The bk are plotted in B for the square wave and in D for the triangular function. It can be seen that the Fourier coefficients of the triangular function drop off more rapidly for higher k than those of the square wave, because the triangular function is more closely related to the sine wave than is the square wave.

Reconstruction of a signal from the Fourier coefficients is demonstrated in Figure 3.2. The individual harmonic oscillations for k = 1, 3, and 5 are shown in Figure 3.2A and C. According to Equation (3.1), the individual oscillations have integer multiples of the frequency and an amplitude decreasing with the reciprocal of the frequency. The summation of the oscillations (Figure 3.2B) shows convergence toward the final shape of the square-wave function. The harmonic coefficients of the triangle wave drop off more rapidly [according to Equation (3.5) with k2], and convergence is much faster, as can be seen in Figure 3.2D.

Fourier synthesis of the square-wave and triangular functions. The individual harmonic contributions (A and C, respectively), when added up, converge toward the original function as shown in B and D, respectively. For k ^<x>, the functions in Figure 3.1A and C emerge. Because the harmonic coefficients drop off toward higher frequencies faster in the triangular wave than in the square wave, the triangular function is better approximated with a lower number of coefficients.

FIGURE 3.2 Fourier synthesis of the square-wave and triangular functions. The individual harmonic contributions (A and C, respectively), when added up, converge toward the original function as shown in B and D, respectively. For k ^<x>, the functions in Figure 3.1A and C emerge. Because the harmonic coefficients drop off toward higher frequencies faster in the triangular wave than in the square wave, the triangular function is better approximated with a lower number of coefficients.

The examples above show the transform of a continuous function into discrete Fourier coefficients. To apply the Fourier transform on discretely sampled signals and images, a discrete Fourier transform needs to be defined where the integration is replaced by a summation over the discrete samples of the signal. The one-dimensional discrete Fourier transform and the inverse transform are defined, respectively, by

tmp2F-95_thumb

where the Fr are the complex Fourier coefficients of the signal fk = f (tk) that is sampled at discrete time pointstmp2F-96_thumb. Equation (3.6) makes use of the Euler relationship ( j is the imaginary unit)

tmp2F-97_thumb

and is therefore closely related to the discrete summation in Equation (3.1), with its orthogonal sine and cosine base. To obtain all possible Fourier coefficients, the summation in Equation (3.6) needs to be evaluated for r = 0,1, …, N-1, where N is the number of data points in the time series.

Implicitly, periodic behavior of the time series fk is assumed; that is, fk = fk+mN, where m may be any integer number. Due to the periodicity of the complex exponential, the same periodicity applies to the Fourier coefficients Fr: namely, Fr = Fr+mN. In fact, the Fourier coefficients are commonly interpreted such that the indices r = 0, 1,2,…, N/2 – 1 are considered coefficients for positive frequencies, and the indices r = N – 1, N – 2, … , N/2 are considered negative frequencies and are shifted by subtracting N from r. The coefficient F0, similar to a0 in Equation (3.1), is the average value of the signal. F0 is always real-valued. The indices r of the Fourier transform can be calibrated to relate to actual frequencies. When the signal has the period T = N A t, the base frequency is 1/T. Therefore, the frequency of the coefficient Fr for r = ±1 is ±1/T, and the frequency of each coefficient Fr is r/T. The highest possible frequency that the discrete Fourier transform provides is N/2T or 1/2A t, in agreement with the Nyquist sampling theorem.

The same considerations apply to images. One horizontal line of pixels provides discrete values fk, sampled with a constant spatial sampling interval A x. Analogous to the time series, the indices r of the Fourier transform represent the spatial frequency ur = r/(MA x). By convention, frequencies in the horizontal direction are given the index u, and frequencies in the vertical direction are given the index v. The two-dimensional Fourier transform of a two-dimensional image is computed by first computing the horizontal one-dimensional Fourier transform line by line, followed by the column-by-column Fourier transform of the horizontal Fourier transforms:

tmp2F-98_thumb

The second summation term represents the horizontal Fourier transform, and the first summation represents the vertical Fourier transform. Since the first exponential term is constant with respect to the second summation, the two exponential terms can be joined, leading to the most commonly used definition for the two-dimensional discrete Fourier transform:

tmp2F-99_thumb

In Equations (3.8) and (3.9), I(x,y) are the image values, M and N are the image size (number of pixels) in the horizontal and vertical directions, respectively, and F(u,v) is the Fourier transform. Notably, the Fourier transform can be seen as an image with the same dimensions as the input image I, albeit with complex pixel values. Like the one-dimensional Fourier transform, those values of u that are larger than M/2 and those values of v that are larger than N/2 are converted to negative frequencies by subtracting M and N, respectively. The result is an image F(u,v) with the origin (u = 0, v = 0) in the center, as illustrated in Figure 3.3. The value at the origin, F(0,0), contains the image average as can be seen from Equation (3.9). A single row of pixels parallel to the x-axis would represent a sine wave with 32 full periods over M = 512 pixels or one period over 16 pixels. A one-dimensional Fourier transform of the single row would have three peaks: a peak in the center at u = 0 (the average value) and two symmetrical peaks 32 pixels right and left of the origin.

Two-dimensional Fourier transform. The spatial-domain image (A) contains periodic sinusoidal intensity variations; the image size is 512 x 512 pixels. Along the y-axis there are eight full waves, and along the x-axis, 32. Its Fourier transform (B, inverted with dark dots representing peak values for better visualization) shows five peaks. The central peak represents the image average. The two peak pairs along the u and v axes are the positive and negative frequency of the sinusoidal intensity variations along the x and y axes, respectively, of the original image. The spatial frequency along the x-axis is higher than that along the y-axis; correspondingly, the peaks on the u-axis are at higher u values (farther away from the origin) than are the peaks on the v-axis.

FIGURE 3.3 Two-dimensional Fourier transform. The spatial-domain image (A) contains periodic sinusoidal intensity variations; the image size is 512 x 512 pixels. Along the y-axis there are eight full waves, and along the x-axis, 32. Its Fourier transform (B, inverted with dark dots representing peak values for better visualization) shows five peaks. The central peak represents the image average. The two peak pairs along the u and v axes are the positive and negative frequency of the sinusoidal intensity variations along the x and y axes, respectively, of the original image. The spatial frequency along the x-axis is higher than that along the y-axis; correspondingly, the peaks on the u-axis are at higher u values (farther away from the origin) than are the peaks on the v-axis.

If the pixel size is known, the peak location translates into a spatial frequency. If we assume that this image has pixels of size 0.2 x 0.2 mm, one period along the x-axis is 3.2 mm. The highest spatial frequency (at u = 256) is 2.5 mm-1 and the spatial frequency at u = 32 is 2.5 mm-1 32/256, or 0.3125 mm-1, the reciprocal of 3.2 mm. Correspondingly, each column in Figure 3.3A shows sinusoidal oscillations of 64 pixels, or 12.8 mm along the y-axis. The peaks on the v-axis can be expected at v = ±8, which corresponds to 2.5 mm-1 8/256 or 0.078125 mm-1, the reciprocal of 12.8 mm. The farther away a pixel in the Fourier transform is from the origin, the higher its associated spatial frequency. In fact, all points with the same distance from the origin (i.e., on concentric circles) have the same spatial frequency. Therefore, any point (u,v) in the frequency-domain plane has the spatial frequency w:

tmp2F-101_thumb

Since the Fourier transform is generally complex-valued, the visual representation is difficult. Most frequently, the magnitude of the frequency components IF(u,v)l is displayed and the phase information omitted. Furthermore, the dynamic range of the values of F(u,v) is huge, and often contrast is improved by displaying IF(u,v)l on a logarithmic scale, remapped to the displayable range from 0 to 255. Many times, the power spectrum, IF(u,v)I,2 is displayed. The power spectrum is a concept that originates in electrical engineering and relates to the distribution of the power of a signal over the frequency.

Demonstration of Fourier-based filtering. After removal of the peaks on the u-axis in Figure 3.3B and inverse Fourier transform, only the oscillations along the y-axis remain.

FIGURE 3.4 Demonstration of Fourier-based filtering. After removal of the peaks on the u-axis in Figure 3.3B and inverse Fourier transform, only the oscillations along the y-axis remain.

The basic idea behind Fourier-domain filters can be derived from the following experiment: If the two outer peaks on the u-axis in Figure 3.3B are removed, only the oscillation along the y-axis should remain. The experimental proof is given in Figure 3.4. The image in Figure 3.3A was Fourier-transformed, and the two symmetrical peaks on the u-axis were removed manually. After performing the inverse Fourier transform, only the oscillations along the y-axis remain. Filtering methods in the frequency domain and their applications are discussed in Section 3.2.

Next post:

Previous post: