分享

傅里叶转换公式及在图像处理中的应用

 pplqingshi 2012-04-13

Fourier Transform 

The Fourier transform is a generalization of the complex Fourier series in the limit as L->infty. Replace the discrete A_n with the continuous F(k)dk while letting n/L->k. Then change the sum to an integral, and the equations become

f(x)=int_(-infty)^inftyF(k)e^(2piikx)dk
(1)
F(k)=int_(-infty)^inftyf(x)e^(-2piikx)dx.
(2)

Here,

F(k)=F_x[f(x)](k)
(3)
=int_(-infty)^inftyf(x)e^(-2piikx)dx
(4)

is called the forward (-i) Fourier transform, and

f(x)=F_k^(-1)[F(k)](x)
(5)
=int_(-infty)^inftyF(k)e^(2piikx)dk
(6)

is called the inverse (+i) Fourier transform. The notation F_x[f(x)](k) is introduced in Trott (2004, p. xxxiv), and f^^(k) and f^_(x) are sometimes also used to denote the Fourier transform and inverse Fourier transform, respectively (Krantz 1999, p. 202).

Note that some authors (especially physicists) prefer to write the transform in terms of angular frequency omega=2pinu instead of the oscillation frequency nu. However, this destroys the symmetry, resulting in the transform pair

H(omega)=F[h(t)]
(7)
=int_(-infty)^inftyh(t)e^(-iomegat)dt
(8)
h(t)=F^(-1)[H(omega)]
(9)
=1/(2pi)int_(-infty)^inftyH(omega)e^(iomegat)domega.
(10)

To restore the symmetry of the transforms, the convention

g(y)=F[f(t)]
(11)
=1/(sqrt(2pi))int_(-infty)^inftyf(t)e^(-iyt)dt
(12)
f(t)=F^(-1)[g(y)]
(13)
=1/(sqrt(2pi))int_(-infty)^inftyg(y)e^(iyt)dy
(14)

is sometimes used (Mathews and Walker 1970, p. 102).

In general, the Fourier transform pair may be defined using two arbitrary constants a and b as

F(omega)=sqrt((|b|)/((2pi)^(1-a)))int_(-infty)^inftyf(t)e^(ibomegat)dt
(15)
f(t)=sqrt((|b|)/((2pi)^(1+a)))int_(-infty)^inftyF(omega)e^(-ibomegat)domega.
(16)

The Fourier transform F(k) of a function f(x) is implemented as FourierTransform[f, x, k], and different choices of a and b can be used by passing the optional FourierParameters-> {a, b} option. By default, Mathematica takes FourierParameters as (0,1). Unfortunately, a number of other conventions are in widespread use. For example, (0,1) is used in modern physics, (1,-1) is used in pure mathematics and systems engineering, (1,1) is used in probability theory for the computation of the characteristic function, (-1,1) is used in classical physics, and (0,-2pi) is used in signal processing. In this work, following Bracewell (1999, pp. 6-7), it is always assumed that a=0 and b=-2pi unless otherwise stated. This choice often results in greatly simplified transforms of common functions such as 1, cos(2pik_0x), etc.

Since any function can be split up into even and odd portions E(x) and O(x),

f(x)=1/2[f(x)+f(-x)]+1/2[f(x)-f(-x)]
(17)
=E(x)+O(x),
(18)

a Fourier transform can always be expressed in terms of the Fourier cosine transform and Fourier sine transform as

 F_x[f(x)](k)=int_(-infty)^inftyE(x)cos(2pikx)dx-iint_(-infty)^inftyO(x)sin(2pikx)dx.
(19)

A function f(x) has a forward and inverse Fourier transform such that

 f(x)={int_(-infty)^inftye^(2piikx)[int_(-infty)^inftyf(x)e^(-2piikx)dx]dk   for f(x) continuous at x; 1/2[f(x_+)+f(x_-)]   for f(x) discontinuous at x,
(20)

provided that

1. int_(-infty)^infty|f(x)|dx exists.

2. There are a finite number of discontinuities.

3. The function has bounded variation. A sufficient weaker condition is fulfillment of the Lipschitz condition

(Ramirez 1985, p. 29). The smoother a function (i.e., the larger the number of continuous derivatives), the more compact its Fourier transform.

The Fourier transform is linear, since if f(x) and g(x) have Fourier transforms F(k) and G(k), then

int[af(x)+bg(x)]e^(-2piikx)dx=aint_(-infty)^inftyf(x)e^(-2piikx)dx+bint_(-infty)^inftyg(x)e^(-2piikx)dx
(21)
=aF(k)+bG(k).
(22)

Therefore,

F[af(x)+bg(x)]=aF[f(x)]+bF[g(x)]
(23)
=aF(k)+bG(k).
(24)

The Fourier transform is also symmetric since F(k)=F_x[f(x)](k) implies F(-k)=F_x[f(-x)](k).

Let f*g denote the convolution, then the transforms of convolutions of functions have particularly nice transforms,

F[f*g]=F[f]F[g]
(25)
F[fg]=F[f]*F[g]
(26)
F^(-1)[F(f)F(g)]=f*g
(27)
F^(-1)[F(f)*F(g)]=fg.
(28)

The first of these is derived as follows:

F[f*g]=int_(-infty)^inftyint_(-infty)^inftye^(-2piikx)f(x^')g(x-x^')dx^'dx
(29)
=int_(-infty)^inftyint_(-infty)^infty[e^(-2piikx^')f(x^')dx^'][e^(-2piik(x-x^'))g(x-x^')dx]
(30)
=[int_(-infty)^inftye^(-2piikx^')f(x^')dx^'][int_(-infty)^inftye^(-2piikx^(''))g(x^(''))dx^('')]
(31)
=F[f]F[g],
(32)

where x^('')=x-x^'.

There is also a somewhat surprising and extremely important relationship between the autocorrelation and the Fourier transform known as the Wiener-Khinchin theorem. Let F_x[f(x)](k)=F(k), and f^_ denote the complex conjugate of f, then the Fourier transform of the absolute square of F(k) is given by

 F_k[|F(k)|^2](x)=int_(-infty)^inftyf^_(tau)f(tau+x)dtau.
(33)

The Fourier transform of a derivative f^'(x) of a function f(x) is simply related to the transform of the function f(x) itself. Consider

 F_x[f^'(x)](k)=int_(-infty)^inftyf^'(x)e^(-2piikx)dx.
(34)

Now use integration by parts

 intvdu=[uv]-intudv
(35)

with

du=f^'(x)dx
(36)
v=e^(-2piikx)
(37)

and

u=f(x)
(38)
dv=-2piike^(-2piikx)dx,
(39)

then

 F_x[f^'(x)](k)=[f(x)e^(-2piikx)]_(-infty)^infty-int_(-infty)^inftyf(x)(-2piike^(-2piikx)dx).
(40)

The first term consists of an oscillating function times f(x). But if the function is bounded so that

 lim_(x->+/-infty)f(x)=0
(41)

(as any physically significant signal must be), then the term vanishes, leaving

F_x[f^'(x)](k)=2piikint_(-infty)^inftyf(x)e^(-2piikx)dx
(42)
=2piikF_x[f(x)](k).
(43)

This process can be iterated for the nth derivative to yield

 F_x[f^((n))(x)](k)=(2piik)^nF_x[f(x)](k).
(44)

The important modulation theorem of Fourier transforms allows F_x[cos(2pik_0x)f(x)](k) to be expressed in terms of F_x[f(x)](k)=F(k) as follows,

F_x[cos(2pik_0x)f(x)](k)=int_(-infty)^inftyf(x)cos(2pik_0x)e^(-2piikx)dx
(45)
=1/2int_(-infty)^inftyf(x)e^(2piik_0x)e^(-2piikx)dx+1/2int_(-infty)^inftyf(x)e^(-2piik_0x)e^(-2piikx)dx
(46)
=1/2int_(-infty)^inftyf(x)e^(-2pii(k-k_0)x)dx+1/2int_(-infty)^inftyf(x)e^(-2pii(k+k_0)x)dx
(47)
=1/2[F(k-k_0)+F(k+k_0)].
(48)

Since the derivative of the Fourier transform is given by

 F^'(k)=d/(dk)F_x[f(x)](k)=int_(-infty)^infty(-2piix)f(x)e^(-2piikx)dx,
(49)

it follows that

 F^'(0)=-2piiint_(-infty)^inftyxf(x)dx.
(50)

Iterating gives the general formula

mu_n=int_(-infty)^inftyx^nf(x)dx
(51)
=(F^((n))(0))/((-2pii)^n).
(52)

The variance of a Fourier transform is

 sigma_f^2=<(xf-<xf>)^2>,
(53)

and it is true that

 sigma_(f+g)=sigma_f+sigma_g.
(54)

If f(x) has the Fourier transform F_x[f(x)](k)=F(k), then the Fourier transform has the shift property

int_(-infty)^inftyf(x-x_0)e^(-2piikx)dx=int_(-infty)^inftyf(x-x_0)e^(-2pii(x-x_0)k)e^(-2pii(kx_0))d(x-x_0)
(55)
=e^(-2piikx_0)F(k),
(56)

so f(x-x_0) has the Fourier transform

 F_x[f(x-x_0)](k)=e^(-2piikx_0)F(k).
(57)

If f(x) has a Fourier transform F_x[f(x)](k)=F(k), then the Fourier transform obeys a similarity theorem.

 int_(-infty)^inftyf(ax)e^(-2piikx)dx=1/(|a|)int_(-infty)^inftyf(ax)e^(-2pii(ax)(k/a))d(ax)=1/(|a|)F(k/a),
(58)

so f(ax) has the Fourier transform

 F_x[f(ax)](k)=|a|^(-1)F(k/a).
(59)

The "equivalent width" of a Fourier transform is

w_e=(int_(-infty)^inftyf(x)dx)/(f(0))
(60)
=(F(0))/(int_(-infty)^inftyF(k)dk).
(61)

The "autocorrelation width" is

w_a=(int_(-infty)^inftyf*f^_dx)/([f*f^_]_0)
(62)
=(int_(-infty)^inftyfdxint_(-infty)^inftyf^_dx)/(int_(-infty)^inftyff^_dx),
(63)

where f*g denotes the cross-correlation of f and g and f^_ is the complex conjugate.

Any operation on f(x) which leaves its area unchanged leaves F(0) unchanged, since

 int_(-infty)^inftyf(x)dx=F_x[f(x)](0)=F(0).
(64)

The following table summarized some common Fourier transform pairs.

functionf(x)F(k)=F_x[f(x)](k)
Fourier transform--11delta(k)
Fourier transform--cosinecos(2pik_0x)1/2[delta(k-k_0)+delta(k+k_0)]
Fourier transform--delta functiondelta(x-x_0)e^(-2piikx_0)
Fourier transform--exponential functione^(-2pik_0|x|)1/pi(k_0)/(k^2+k_0^2)
Fourier transform--Gaussiane^(-ax^2)sqrt(pi/a)e^(-pi^2k^2/a)
Fourier transform--Heaviside step functionH(x)1/2[delta(k)-i/(pik)]
Fourier transform--inverse function-PV1/(pix)i[1-2H(-k)]
Fourier transform--Lorentzian function1/pi(1/2Gamma)/((x-x_0)^2+(1/2Gamma)^2)e^(-2piikx_0-Gammapi|k|)
Fourier transform--ramp functionR(x)piidelta^'(2pik)-1/(4pi^2k^2)
Fourier transform--sinesin(2pik_0x)1/2i[delta(k+k_0)-delta(k-k_0)]

In two dimensions, the Fourier transform becomes

F(x,y)=int_(-infty)^inftyint_(-infty)^inftyf(k_x,k_y)e^(-2pii(k_xx+k_yy))dk_xdk_y
(65)
f(k_x,k_y)=int_(-infty)^inftyint_(-infty)^inftyF(x,y)e^(2pii(k_xx+k_yy))dxdy.
(66)

Similarly, the n-dimensional Fourier transform can be defined for k, x in R^n by

F(x)=int_(-infty)^infty...int_(-infty)^infty_()_(n)f(k)e^(-2piik·x)d^nk
(67)
f(k)=int_(-infty)^infty...int_(-infty)^infty_()_(n)F(x)e^(2piik·x)d^nx.
(68)
应用:

Fourier Transform


Common Names: Fourier Transform, Spectral Analysis, Frequency Analysis

Brief Description

The Fourier Transform is an important image processing tool which is used to decompose an image into its sine and cosine components. The output of the transformation represents the image in the Fourier or frequency domain, while the input image is the spatial domain equivalent. In the Fourier domain image, each point represents a particular frequency contained in the spatial domain image.

The Fourier Transform is used in a wide range of applications, such as image analysis, image filtering, image reconstruction and image compression.

How It Works

As we are only concerned with digital images, we will restrict this discussion to the Discrete Fourier Transform (DFT).

The DFT is the sampled Fourier Transform and therefore does not contain all frequencies forming an image, but only a set of samples which is large enough to fully describe the spatial domain image. The number of frequencies corresponds to the number of pixels in the spatial domain image, i.e. the image in the spatial and Fourier domain are of the same size.

For a square image of size N×N, the two-dimensional DFT is given by:

Eqn:eqnfour1

where f(a,b) is the image in the spatial domain and the exponential term is the basis function corresponding to each point F(k,l) in the Fourier space. The equation can be interpreted as: the value of each point F(k,l) is obtained by multiplying the spatial image with the corresponding base function and summing the result.

The basis functions are sine and cosine waves with increasing frequencies, i.e. F(0,0) represents the DC-component of the image which corresponds to the average brightness and F(N-1,N-1) represents the highest frequency.

In a similar way, the Fourier image can be re-transformed to the spatial domain. The inverse Fourier transform is given by:

Eqn:eqnfour2

Note the Eqn:oneovern2 normalization term in the inverse transformation. This normalization is sometimes applied to the forward transform instead of the inverse transform, but it should not be used for both.$$

To obtain the result for the above equations, a double sum has to be calculated for each image point. However, because the Fourier Transform is separable, it can be written as

Eqn:eqnfour3
where
Eqn:eqnfour4
Using these two formulas, the spatial domain image is first transformed into an intermediate image using N one-dimensional Fourier Transforms. This intermediate image is then transformed into the final image, again using N one-dimensional Fourier Transforms. Expressing the two-dimensional Fourier Transform in terms of a series of 2N one-dimensional transforms decreases the number of required computations.

Even with these computational savings, the ordinary one-dimensional DFT has Eqn:eqnfour6 complexity. This can be reduced to Eqn:eqnfour5 if we employ the Fast Fourier Transform (FFT) to compute the one-dimensional DFTs. This is a significant improvement, in particular for large images. There are various forms of the FFT and most of them restrict the size of the input image that may be transformed, often to Eqn:eqnfour7 where n is an integer. The mathematical details are well described in the literature.

The Fourier Transform produces a complex number valued output image which can be displayed with two images, either with the real and imaginary part or with magnitude and phase. In image processing, often only the magnitude of the Fourier Transform is displayed, as it contains most of the information of the geometric structure of the spatial domain image. However, if we want to re-transform the Fourier image into the correct spatial domain after some processing in the frequency domain, we must make sure to preserve both magnitude and phase of the Fourier image.

The Fourier domain image has a much greater range than the image in the spatial domain. Hence, to be sufficiently accurate, its values are usually calculated and stored in float values.

Guidelines for Use

The Fourier Transform is used if we want to access the geometric characteristics of a spatial domain image. Because the image in the Fourier domain is decomposed into its sinusoidal components, it is easy to examine or process certain frequencies of the image, thus influencing the geometric structure in the spatial domain.

In most implementations the Fourier image is shifted in such a way that the DC-value (i.e. the image mean) F(0,0) is displayed in the center of the image. The further away from the center an image point is, the higher is its corresponding frequency.

We start off by applying the Fourier Transform of

cln1

The magnitude calculated from the complex result is shown in

cln1fur1

We can see that the DC-value is by far the largest component of the image. However, the dynamic range of the Fourier coefficients (i.e. the intensity values in the Fourier image) is too large to be displayed on the screen, therefore all other values appear as black. If we apply a logarithmic transformation to the image we obtain

cln1fur2

The result shows that the image contains components of all frequencies, but that their magnitude gets smaller for higher frequencies. Hence, low frequencies contain more image information than the higher ones. The transform image also tells us that there are two dominating directions in the Fourier image, one passing vertically and one horizontally through the center. These originate from the regular patterns in the background of the original image.

The phase of the Fourier transform of the same image is shown in

cln1fur3

The value of each point determines the phase of the corresponding frequency. As in the magnitude image, we can identify the vertical and horizontal lines corresponding to the patterns in the original image. The phase image does not yield much new information about the structure of the spatial domain image; therefore, in the following examples, we will restrict ourselves to displaying only the magnitude of the Fourier Transform.

Before we leave the phase image entirely, however, note that if we apply the inverse Fourier Transform to the above magnitude image while ignoring the phase (and then histogram equalize the output) we obtain

cln1fil1

Although this image contains the same frequencies (and amount of frequencies) as the original input image, it is corrupted beyond recognition. This shows that the phase information is crucial to reconstruct the correct image in the spatial domain.

We will now experiment with some simple images to better understand the nature of the transform. The response of the Fourier Transform to periodic patterns in the spatial domain images can be seen very easily in the following artificial images.

The image

stp2

shows 2 pixel wide vertical stripes. The magnitude of the Fourier transform of this image is shown in

stp2fur1

If we look carefully, we can see that it contains 3 main values: the DC-value and, since the Fourier image is symmetrical to its center, two points corresponding to the frequency of the stripes in the original image. Note that the two points lie on a horizontal line through the image center, because the image intensity in the spatial domain changes the most if we go along it horizontally.

The distance of the points to the center can be explained as follows: the maximum frequency which can be represented in the spatial domain are two pixel wide stripe pairs (one white, one black).

Eqn:eqnfour8
Hence, the two pixel wide stripes in the above image represent
Eqn:eqnfour9
Thus, the points in the Fourier image are halfway between the center and the edge of the image, i.e. the represented frequency is half of the maximum.

Further investigation of the Fourier image shows that the magnitude of other frequencies in the image is less than Eqn:eqnfoura of the DC-value, i.e. they don't make any significant contribution to the image. The magnitudes of the two minor points are each two-thirds of the DC-value.

Similar effects as in the above example can be seen when applying the Fourier Transform to

stp1

which consists of diagonal stripes. In

stp1fur1

showing the magnitude of the Fourier Transform, we can see that, again, the main components of the transformed image are the DC-value and the two points corresponding to the frequency of the stripes. However, the logarithmic transform of the Fourier Transform,

stp1fur2

shows that now the image contains many minor frequencies. The main reason is that a diagonal can only be approximated by the square pixels of the image, hence, additional frequencies are needed to compose the image. The logarithmic scaling makes it difficult to tell the influence of single frequencies in the original image. To find the most important frequencies we threshold the original Fourier magnitude image at level 13. The resulting Fourier image,

stp1fur3

shows all frequencies whose magnitude is at least 5% of the main peak. Compared to the original Fourier image, several more points appear. They are all on the same diagonal as the three main components, i.e. they all originate from the periodic stripes. The represented frequencies are all multiples of the basic frequency of the stripes in the spatial domain image. This is because a rectangular signal, like the stripes, with the frequency Eqn:eqnfourb is a composition of sine waves with the frequencies Eqn:eqnfourc, known as the harmonics of Eqn:eqnfourb. All other frequencies disappeared from the Fourier image, i.e. the magnitude of each of them is less than 5% of the DC-value.

A Fourier-Transformed image can be used for frequency filtering. A simple example is illustrated with the above image. If we multiply the (complex) Fourier image obtained above with an image containing a circle (of r = 32 pixels), we can set all frequencies larger than Eqn:eqnfourb to zero as shown in the logarithmic transformed image

stp1fur5

By applying the inverse Fourier Transform we obtain

stp1fil1

The resulting image is a lowpass filtered version of the original spatial domain image. Since all other frequencies have been suppressed, this result is the sum of the constant DC-value and a sine-wave with the frequency Eqn:eqnfourb. Further examples can be seen in the worksheet on frequency filtering.

A property of the Fourier Transform which is used, for example, for the removal of additive noise, is its distributivity over addition. We can illustrate this by adding the complex Fourier images of the two previous example images. To display the result and emphasize the main peaks, we threshold the magnitude of the complex image, as can be seen in

stp1fur4

Applying the inverse Fourier Transform to the complex image yields

stp1fil2

According to the distributivity law, this image is the same as the direct sum of the two original spatial domain images.

Finally, we present an example (i.e. text orientation finding) where the Fourier Transform is used to gain information about the geometric structure of the spatial domain image. Text recognition using image processing techniques is simplified if we can assume that the text lines are in a predefined direction. Here we show how the Fourier Transform can be used to find the initial orientation of the text and then a rotation can be applied to correct the error. We illustrate this technique using

son3

a binary image of English text. The logarithm of the magnitude of its Fourier transform is

son3fur2

and

son3fur4

is the thresholded magnitude of the Fourier image. We can see that the main values lie on a vertical line, indicating that the text lines in the input image are horizontal.
If we proceed in the same way with

son3rot1

which was rotated about 45°, we obtain

son3fur1

and

son3fur3

in the Fourier space. We can see that the line of the main peaks in the Fourier domain is rotated according to rotation of the input image. The second line in the logarithmic image (perpendicular to the main direction) originates from the black corners in the rotated image.

Although we managed to find a threshold which separates the main peaks from the background, we have a reasonable amount of noise in the Fourier image resulting from the irregular pattern of the letters. We could decrease these background values and therefore increase the difference to the main peaks if we were able to form solid blocks out of the text-lines. This could, for example, be done by using a morphological operator.

Common Variants

Another sinusoidal transform (i.e. transform with sinusoidal base functions) related to the DFT is the Discrete Cosine Transform (DCT). For an N×N image, the DCT is given by

Eqn:eqnfourd
with
Eqn:eqnfoure

The main advantages of the DCT are that it yields a real valued output image and that it is a fast transform. A major use of the DCT is in image compression --- i.e. trying to reduce the amount of data needed to store an image. After performing a DCT it is possible to throw away the coefficients that encode high frequency components that the human eye is not very sensitive to. Thus the amount of data can be reduced, without seriously affecting the way an image looks to the human eye.


    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息,谨防诈骗。如发现有害或侵权内容,请点击一键举报。
    转藏 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多