2 of 107

Fourier Transform

You have so far studied the Fourier transform of a 1D or 2D continuous (analog) function.

The functions we deal with in practical signal or image processing are however discrete.

We need an equivalent of the Fourier transform of such discrete signals.

3 of 107

Towards the Discrete Fourier Transform

Given a continuous signal f(t), we will consider signal g(t), obtained by sampling f(t) at equally-spaced time instants.

So g(t) is the pointwise product of f(t) with an impulse train s_ΔT(t) given by:

The impulse functions here are Dirac delta functions.

4 of 107

Towards the Discrete Fourier Transform

This shows that the Fourier transform G(μ) of the sampled signal g(t) is an infinite periodic sequence of copies of F(μ) which is the Fourier transform of f(t). The periodicity of G(μ) is given by the sampling period ΔT. Also though g(t) is a sampled function, its Fourier transform is continuous.

Check the book for the derivation for this (see Appendix at the end of the slides)

5 of 107

G(μ)

6 of 107

Towards the Discrete Fourier Transform

Let us now try to derive G(μ) directly in terms of g(t) instead of F(μ).

f(t)δ(t-nΔT) is the sampled function (a Dirac impulse weighted by f(t)). To get the actual sample values (the strength of the impulse), we do an integration, i.e. of the following form:

7 of 107

Towards the Discrete Fourier Transform

We have seen earlier that the Fourier transform of g(t) is periodic with period 1/ΔT, so it is of interest to characterize one period.

Consider that we take some M equally spaced samples of G(μ) over one period from μ = 0 to 1/ΔT.

This gives us frequencies of the form μ = u/MΔT where u = 0,1,2…,M-1.

8 of 107

Towards the Discrete Fourier Transform

Plugging these values for μ we now have:

This is called the Discrete Fourier Transform of the sequence {f_n} where n = 0,1,…,M-1.

Given the sequence {F_d(u)} where u = 0,1,…M-1, we can recover {f_n} using:

This is the inverse Discrete Fourier Transform (IDFT).

9 of 107

Towards the Discrete Fourier Transform

Consider:

It can be proved that plugging in the expression for f(n) into the expression for F_d(u) yields the identity F_d(u) = F_d(u).

Also plugging in the expression for F_d(u) into the expression for f(n) yields the identity f(n) = f(n).

10 of 107

Towards the Discrete Fourier Transform

In some books, the following expressions are used:

Note that the above expressions can be written in the following matrix vector format:

Vectors of size M by 1

U turns out to be an orthonormal M by M matrix – called the discrete Fourier basis matrix or DFT matrix. The square root of M in the denominator is required for U to be orthonormal, else it would be proportional to an orthonormal matrix.

MATLAB does NOT follow this convention – it follows the one on the earlier slide.

11 of 107

Vectors of size M by 1

U is an orthonormal M by M matrix – called the discrete Fourier basis matrix or DFT matrix. The square root of M is required for U to be orthonormal, else it would be proportional to an orthonormal matrix.

Notice in the last equality how the signal f is being represented as a linear combination of column vectors of the DFT matrix. The coefficients of the linear combination are the discrete Fourier coefficients!

12 of 107

Sampling in time and frequency

Remember that the function g(t) was created by sampling f(t) with a period of ΔT.

And the spacing between the samples in the frequency domain (to get the DFT from G(μ)) is

1/ MΔT since μ = u/MΔT where u = 0,1,2…,M-1.

Likewise the range of frequencies spanned by the DFT is also inversely proportional to ΔT.

13 of 107

DFT properties

Linearity: F(af+bg)(u) = aF(f)(u) + bF(g)(u)

Periodicity: F(u) = F(u + kM) for integer k,

And so is the inverse DFT since f(n) = f(n+kM).

The DFT is in general complex. Hence it has a magnitude |F(u)| and a phase.

14 of 107

Clarification about DFT

We have seen earlier that the Fourier transform G(μ) of the sampled version g(t) of analog signal f(t) is continuous.

We also saw earlier that we take some M equally spaced samples of G(μ) over one period from μ = 0 to 1/ΔT.

This way the DFT and the discrete signal both were vectors of M elements.

Obtaining the DFT given the signal (or the signal given its DFT) is an efficient operation owing to the orthonormality of the discrete fourier matrix (why efficient? Because for an orthonormal matrix, inverse = transpose).

Why don’t we take more than M samples in the frequency domain?

If we did, the aforementioned computational efficiency would be lost. The inverse transform would require a matrix pseudo-inverse which is costly.

Also the columns of the orthonormal matrix U^T(size M x M) constitute a basis: i.e. any vector in M-dimensional space can be uniquely represented using linear combination of the columns of that matrix. If you took more than M samples in the frequency domain, that uniqueness would be lost as U^Twould now have size M’ x M where M’ > M.

15 of 107

Convolution of discrete signals

Discrete equivalent of the convolution is:

Due to the periodic nature of the DFT and IDFT, the convolution will also be periodic.

The above formula represents one period of that periodic resultant signal.

16 of 107

Convolution of discrete signals

The convolution theorem from continuous signals has an extension to the discrete case as well:

Therefore discrete convolution can be implemented using product of DFTs followed by an IDFT.

But in doing so, one has to account for periodicity issues to avoid wrap-around artifacts (see next slide).

17 of 107

Convolution of discrete signals

Consider the discrete convolution:

To convolve f with h, you need to (1) flip h about the origin, (2) translate the flipped signal by an amount n, and (3) compute the sum in the above formula.
Steps 2 and 3 are repeated for each value of n.
The variable n ranges over all integers required to completely slide h over f.
If h and f have size M, then n has to range up to 2M-1.

18 of 107

Convolution of discrete signals

In other words, the resultant signal must have length of 2M-1.

The convolution can be implemented in the time/spatial domain – and MATLAB has a routine called conv which does the job for you!

Now imagine you tried to implement the convolution using a DFT of M-point signals followed by an M-point IDFT.

Owing to the assumed periodicity, you would get an undesirable wrap-around effect. See next slide for an illustration.

20 of 107

Convolution of discrete signals

How does one deal with this conundrum?

If f has length M and h has length K, then you should zero-pad both sequences so that they now have length at least M+K-1.

Then compute M+K-1 point DFT for both, multiply the DFTs and compute the M+K-1 point IDFT.

21 of 107

Convolution of discrete signals: 2 methods

Consider you want to convolve signal f having N elements with signal h having K elements.

You can use the conv routine in MATLAB – which by default produces a signal of N+K-1 elements for you. It takes care of zero-padding internally for you.

The other (equivalent) alternative is to:
Append f with K-1 zeros.
Append h with N-1 zeros.
Compute the N+K-1 point DFT of f and h (using fft).
Multiply the two DFTs point-wise.
Compute the IDFT of the result – this gives you a signal with N+K-1 elements.

In both cases, you may wish to extract the first N elements of the resultant signal (note that the trailing K-1 elements may not be zeros!).

22 of 107

Why convolution using Fourier transforms?

The time complexity of computing the convolution of a signal of length M with another of length M is O(M²).

With a DFT computed naively, it would remain O(M²) – the complexity of multiplying a M x M matrix with a M x 1 vector.

23 of 107

Why convolution using Fourier transforms?

But there’s a smarter way of doing the same which computes the DFT in O(M log M) time.

That is called the Fast Fourier Transform, discovered by Cooley and Tukey in 1965.

It is based on a divide and conquer algorithmic strategy.

The same strategy can be used to compute the IDFT in O(M log M) time.

Hence convolution can now be computed in O(M log M) time.

24 of 107

Towards the Fast Fourier Transform

Even indices

Odd indices

25 of 107

Towards the Fast Fourier Transform

The M-point DFT computation for F is split up into two halves. The first half requires two M/2-point DFT computations – one for F_even, one for F_odd. The second half follows directly without any additional transform evaluations once F_evenand F_oddare computed. Now F_evenand F_oddcan be further split up recursively.

26 of 107

Towards the Fast Fourier Transform

This gives:

There is also a fast inverse Fourier transform which works quite similarly in O(M log M) time. The speedup achieved by the fast fourier transform over a naïve DFT computation is rather dramatic for large M!

T(M)

27 of 107

MATLAB implementation

The Fast Fourier Transform is implemented in MATLAB directly – there are the routines fft and ifft for the inverse.

28 of 107

2D-DFT

Given a 2D discrete signal (image) f(x,y) of size W₁ by W₂, its DFT is given as:

29 of 107

2D-DFT computation

Given a 2D discrete signal (image) f(x,y) of size W₁ by W₂, its DFT is given as:

Compute row-wise FFT, followed by a column-wise FFT. Overall complexity is W₁W₂ log(W₁W₂).

30 of 107

MATLAB implementation

The FFT for 2D arrays or images is implemented in MATLAB routines fft2 and ifft2.

31 of 107

2D-DFT properties

Linearity: F(af+bg)(u,v) = aF(f)(u,v) + bF(g)(u,v) where a and b are scalars.

Periodicity: F(u,v) = F(u + k₁W₁,v) = F(u,v + k₂W₂)= F(u + k₁W₁,v + k₂W₂) for integer k₁ and k₂.

And so is the inverse DFT since f(x,y) = f(x + k₁W₁,y) = f(x,y + k₂W₂)= f(x + k₁W₁,y + k₂W₂) for integer k₁ and k₂.

32 of 107

2D-DFT: Magnitude and phase

The phase carries very critical information.

If you retain the magnitude of the Fourier transform of an image, but change its phase, the appearance drastically changes.

33 of 107

IFT (Magnitude of FT of Salman Khan times phase of P V Sindhu)

IFT (Magnitude of FT of P V Sindhu times phase of Salman Khan)

34 of 107

2D-DFT properties

Translation:

Note that translation of a signal does not change the magnitude of the Fourier transform – only its phase.

35 of 107

2D-DFT properties

The Fourier transform of a rotated image yields a rotated version of its Fourier transform

If you consider conversion to polar coordinates, then we have:

36 of 107

2D-DFT visualization in MATLAB

In the formula for the DFT, the frequency range from u = 0 to u = M-1 actually represents two half periods back to back meeting at M/2 (see next slide).
It is more convenient to instead view a complete period of the DFT instead.
For that purpose we visualize F(u-M/2) instead of F(u) in 1D, and F(u-W₁/2,v-W₂/2) instead of F(u,v) in 2D.
Thereby frequency u = 0 or u = v = 0 now occurs at the center of the displayed Fourier transform.
Note that F(u-W₁/2,v-W₂/2) = F(f(x,y)(-1)^x⁺^y)(u,v), so the centering operation is easy to implement.

37 of 107

One full period

Two half periods – back to back

39 of 107

2D-DFT visualization in MATLAB

To visualize the DFT of an image, we visualize the magnitude of the DFT.

The Fourier transform is first centered as mentioned on the previous slide.

Due to the large range of magnitude values, the DFT magnitudes are visualized on a logarithmic scale, i.e. we view log(|F(u,v)|+1) where the 1 is added for stability.

40 of 107

fim2 = fft2(im2);

absfim2 = log(abs(fim2)+1);

imshow(absfim1,[-1 18]);

colormap (jet); colorbar;

The (0,0) frequency lies at the

top left corner

Without centering the fft

fim2 = fftshift(fft2(im2));

absfim2 = log(abs(fim2)+1);

imshow(absfim1,[-1 18]);

colormap (jet); colorbar;

The (0,0) frequency lies in the

Center.

With centering the fft

41 of 107

“Viewing” the rotation property of the DFT

42 of 107

Visualizing DFT bases in 2D

http://en.wikipedia.org/wiki/JPEG

Real part of DFT

Discrete cosine transform (for comparison only)

Each image represents a column of the 2D DFT matrix reshaped to form an image. Only the real component is displayed here. Our theory says that any discrete image can be accurately represented as a linear combination of such basis images. The coefficients of the linear combination may be positive, negative or zero, and real or complex in value. Note that pixel (x,y) of the image for frequency (u,v) represents: Exp(j 2*pi*(ux + vy)/N)

43 of 107

2D convolution

In MATLAB, 2D convolution can be implemented using the routine conv2.

This can be very expensive if the signals you wish to convolve another with, are of large size.

Hence one resorts to the convolution theorem which holds in 2D as well.

44 of 107

2D convolution

The convolution theorem applies to 2D-DFTs as well:

Consider an image f of size W₁ x W₂ which you want to convolve with a kernel k of size K₁ x K₂ using the DFT method.
Then you should symmetrically zero-pad f and k so that they acquire size (W₁+K₁-1) x (W₂+K₂-1).
Compute the DFTs of the zero-padded images using the FFT algorithm, point-wise multiply them and obtain the IDFT of the result.
Extract the central W₁ x W₂ portion of the result – for the final answer.

45 of 107

Image Filtering: Frequency domain

You have studied image filters of various types: mean filter, Gaussian filter, bilateral filter, patch-based filter.

The former two are linear filters and the latter two aren’t.

Linear filters are represented using convolutions and hence have a frequency domain interpretation as seen on the previous slides.

46 of 107

Image Filtering: Frequency domain

Hence such filters can also be designed in the frequency domain.

In certain applications, it is in fact more intuitive to design filters directly in the frequency domain.

Why? Because you get to design directly which frequency components to weaken/eliminate and which ones to boost, and by how much.

47 of 107

Low pass filters

Edges and fine textures in images contribute significantly to the high frequency content of an image.

When you smooth/blur an image, these edges and textures are weakened (or removed).

Such filters allow only the low frequencies to remain intact and are called as low pass filters.

In the frequency domain, an ideal low pass filter can be represented as follows:

Note: we are assuming (0,0) to be the center (lowest) frequency. Frequencies outside a radius of D from the center frequency are eliminated.

48 of 107

Low pass filters

49 of 107

Low pass filters

To apply such a filter to an image f to create a filtered image g, we do as follows:

D is a design parameter of the filter – often called the cutoff frequency.

This is called the ideal low pass filter as it completely eliminates frequencies outside the chosen radius.

50 of 107

Low pass filters

Notice the ringing artifacts around the edges!

51 of 107

Low pass filters

The ringing artifacts can be explained by the convolution theorem.

The corresponding spatial domain filter is called the jinc function.

The jinc is a 2D and circularly symmetric version of the sinc function – see also the next slide.

Any cross section of the jinc is basically a sinc function.

https://en.wikipedia.org/wiki/Sombrero_function

52 of 107

An image can be regarded as a weighted sum of Kronecker delta functions (impulses). Convolving with a sinc function means placing copies of the sinc function at each impulse. The central lobe of the sinc function causes the blurring and the smaller lobes give rise to the “ripple artifacts” or ringing.

53 of 107

Low pass filters: other types

These ringing artifacts can be quite undesirable.

Hence the ideal low pass filter is replaced by other types of low pass filters which weaken but do not totally eliminate the higher frequencies.

For example:

Butterworth filter:

D = cutoff frequency, n = order of the filter

Gaussian filter

55 of 107

Low pass filters: Butterworth

As D increases, the Butterworth filter allows more frequencies to pass through without weakening (see previous slide).

As the order n increases, the Butterworth filter weakens the higher frequencies more aggressively (why?) – which actually increases ringing artifacts (see previous slide).

57 of 107

Low pass filters: Gaussian

The spatial domain representation for the Gaussian LPF is basically a Gaussian.

In fact, the following are Fourier transform pairs:

http://www.cse.yorku.ca/~kosta/CompVis_Notes/fourier_transform_Gaussian.pdf

Notice from the equations that a spatial domain Gaussian with high standard deviation corresponds to a frequency domain Gaussian with low cut off frequency! The extreme case of this is a constant intensity signal (Gaussian with very high standard deviation – infinite as such) whose Fourier transform is an impulse at the origin of the Frequency plane). Another extreme case is a spatial domain signal which is just an impulse – its Fourier transform has constant amplitude everywhere.

58 of 107

High pass filter

It is a filter that allows high frequencies and eliminates or weakens the lower frequencies.

The equation for the frequency response of an ideal high pass filter is given as:

In fact, given an LPF, an HPF can be constructed from it using:

61 of 107

Ringing artifacts for an Ideal HPF!

Ringing artifacts gone!

62 of 107

Ringing artifacts for an Ideal HPF!

Ringing artifacts gone!

63 of 107

High pass filters

Note that any gradient-based operation on images is essentially a type of high-pass filter.

This includes first order derivative filters or the Laplacian filter.

Why? Consider the first order derivative filter in y direction.

64 of 107

High pass filters

Consider the Laplacian filter:

65 of 107

Boosting high frequencies

In some applications, you neither wish to weaken the higher frequencies nor the lower frequencies.

For example, in image sharpening applications you wish to boost the edges or textures (basically higher frequencies).

In that case you wish to perform operations of the following form:

67 of 107

Notch filters

There are applications where images are presented with periodic noise patterns – e.g. images scanned from old newspapers (see next slide).

Usually, the amplitudes of the Fourier components of natural images are seen to decay with frequency.

The Fourier transform of images with periodic patterns show unnatural peaks – unlike ordinary natural images.

Such images can be restored using notch filters, which basically weaken or eliminate these unnatural frequency components (or any frequency components specified by the user).

68 of 107

The yellow ellipses indicate unnatural peaks in the Fourier transform of the image with the superimposed interfering pattern. These are unnatural because typical images do not have such Fourier transforms. Therefore one can apply a notch filter to remove the interfering pattern.

70 of 107

Courtesy: http://www.robots.ox.ac.uk/~az/lectures/ia/lect2.pdf

71 of 107

Courtesy: http://www.robots.ox.ac.uk/~az/lectures/ia/lect2.pdf

72 of 107

Notch filters

The expression for a notch filter frequency response is given as follows (Q = number of frequency components to suppress):

Ideal notch reject filter

Gaussian notch reject filter

73 of 107

Algorithm for frequency domain filtering

Consider you have to filter an image f of size H x W using one of the filters described in these slides.
Zero pad the image symmetrically to create a new image f’ of size 2H x 2W.
Compute the Fourier transform of f’ and center it using fftshift. Let the Fourier transform be F’(u,v).
Design a frequency domain filter H(u,v) of size 2H x 2W with the zero frequency at index (H,W) of the 2D filter array.
Compute the product F’(u,v)H(u,v).
Compute the inverse Fourier transform of the product after applying ifftshift (to undo the effect of fftshift).
Consider only the real part of the inverse Fourier transform and extract the central portion of size H x W.
That gives you the final filtered image.
Note that the zero padding is important. To understand the difference, see the results on the next slide with and without zero-padding.

74 of 107

Effect of Gaussian LPF with appropriate zero padding

Original image

Effect of Gaussian LPF without appropriate zero padding – zero the border artifacts!

75 of 107

Interpreting the DFT of Natural Images

In a few cases, one can guess some structural properties of an image by looking at its Fourier transform.

Most natural images have stronger low frequency components as compared to higher frequency components. This is generally true, although its not a strictly monotonic relationship.

For example, how does the DFT of a vertical edge image look like?

76 of 107

There is a sinc function along the u axis!

77 of 107

Strong horizontal and vertical edges = strong responses along u and v axes in the Fourier plane!

A strong edge in direction θ in XY plane = a strong response in the direction perpendicular to θ in UV plane

78 of 107

Interpreting the DFT

For an image (or any 2D array) with Fourier transform F(u,v), F(0,0) is proportional to the average value of the image intensity – why?

The DFT of a non-zero constant intensity image is basically an impulse at the zero frequency. What is the height of the impulse?

79 of 107

DFT limitations

A strong response at frequency (u,v) (i.e. a large magnitude for F(u,v)) indicates the presence of a sinusoid with frequency (u,v) and with large amplitude somewhere in the image.

In general, the Fourier transform cannot and does not tell us where in the image the sinusoid is located (read section 1 of http://users.rowan.edu/~polikar/WAVELETS/WTpart1.html ).

For that, you need to compute separate Fourier transforms over smaller regions of the image (Short time Fourier transform) – that will tell you which region(s) contained a particular sinusoidal component.

80 of 107

Fun with Fourier: Hybrid images

Hybrid images are a superposition of an image J₁ with strong low frequency components and weak higher frequency components, and an image J₂ with stronger higher frequency components and weak lower frequency components.

J₁ = apply LPF with some cutoff frequency D on an image.
J₂ = apply HPF with some cutoff frequency D on another image.

The appearance of the hybrid image changes with viewing distance!

81 of 107

Tiger – when seen up-close, Cheetah – when seen from far

Einstein from nearby,

Marilyn Monroe from far away (or if you squint)!

Refer to:

http://cvcl.mit.edu/hybrid_gallery/monroe_einstein.html

http://cvcl.mit.edu/publications/OlivaTorralb_Hybrid_Siggraph06.pdf (paper from a graphics journal!)

82 of 107

Fourier in action: tomographic reconstruction

Tomography refers to the process of estimation of the internal structure of an object (tomos in Greek means interior).

This is different from photographic imaging which images only the object surface.

Tomography is typically accomplished via X Rays and the so called Radon transform.

83 of 107

Radon transform

Imagine a line was drawn through the 2D image in a certain direction α, and you integrated the intensity values along that line.

Now you repeat this for lines parallel to the original one but at different offsets.

Each such summation produces a “bin” of the “tomographic projection”.

The collection of bins form a 1D array which is called the tomographic projection or the Radon transform of the object in the direction α.

84 of 107

Reconstruction from the Radon Transform

The Radon transform computation can be similarly repeated for multiple angles.

The aim of tomographic reconstruction is to reconstruct the 2D image from a collection of its 1D tomographic projections in different angles.

85 of 107

Fourier transform of the Radon Transform

The Radon transform is given as:

Its 1D Fourier transform w.r.t. ρ (keeping ϴ fixed to some value) is given by:

G(μ, ϴ) is the Fourier transform of the projection of f(x,y) along some direction ϴ.

86 of 107

Fourier transform of the Radon Transform

Its 1D Fourier transform w.r.t. ρ (keeping ϴ fixed) is given by:

The RHS of this equation is a slice of the 2D Fourier transform of f(x,y), i.e. F(u,v), along the angle ϴ in the frequency plane, and passing through the origin

This equation above is called the Projection Slice Theorem or the Fourier Slice Theorem. It states that the Fourier transform of a projection of the 2D object along some direction ϴ (i.e. G(μ, ϴ)) is equal to a slice of the 2D Fourier transform of the object along the same direction ϴ (in the frequency plane), passing through the origin.

87 of 107

Reconstruction from the Radon transfer

The projection slice theorem thus tells you that you can fill up the 2D Fourier transform of the original image – since the 1D Fourier transform of the projection at angle θ is equal to a slice of the 2D Fourier transform of the image at the angle θ through the origin of the (u,v) plane.

Thus given different projections, each at different angles, you can fill up the Fourier transform of the original image. The final inverse Fourier transform gives you a reconstruction of the image.

88 of 107

How does the Radon transform show up in tomographic acquisition?

Let I₀ be the intensity of the X-ray beam sent through the object f at angle θ.
The materials in the object absorb the X-rays partially.
The remaining energy is sent out and is measured by a detector.
As per Beer’s law, the energy measured at the detector is given by:

Measured projection vector (many bins) at angle θ

Radon transform of object f at angle θ

89 of 107

Fourier transforms in actual acquisition

Computed tomography (in an indirect way via the Fourier slice theorem)
Magnetic resonance imaging – the machine measures the Fourier transform of the object placed inside!
Optics: in diffraction, only the squared magnitudes of the Fourier transform are measured – phase retrieval problem.

90 of 107

Applications of the Fourier transform in image processing

Filtering: efficient implementation of the convolution
Filtering: frequency domain filters
Others: solving partial differential equations (not covered).

91 of 107

Appendix: Fourier transform of an impulse train

Impulse train is given by:

It is a periodic signal with period ∆T, and hence it can be expressed by a Fourier series:

92 of 107

Appendix: Fourier transform of an impulse train

Simplifying the Fourier series coefficients, we have:

Thus the Fourier transform of an impulse train with period ∆T is also an impulse train with period 1/∆T. (Compare: The Fourier transform of a Gaussian is also a Gaussian with an inverted standard deviation). There are a few other functions that are their own Fourier transform.

93 of 107

Appendix: Sampling Theorem

This is a very important theorem in signal processing.

Different parts of it were independently discovered by Shannon, Nyquist, Whittaker and Kotelnikov.

It is often called Shannon’s sampling theorem or the Shannon-Nyquist sampling theorem, but it may well be called the Shannon-Nyquist-Whittaker-Kotelnikov sampling theorem.

94 of 107

Appendix: Sampling Theorem

The theorem states that any band-limited signal with maximum frequency B can be accurately reconstructed from its equi-spaced time-domain samples provided the sampling rate is more than 2B (equivalently the sampling period is less than 1/2B).

The reconstruction of the original (analog) signal from its time-domain samples proceeds by sinc interpolation – we will see it later.

An intuitive proof of the sampling theorem can be seen from the next four slides (pictorially).

95 of 107

Appendix: Sampling Theorem

f(t) is an underlying analog signal which is bandlimited, and whose equi-spaced time-domain samples we observe (second figure). The aim is to recover f(t) accurately from its sampled version. The gap between the samples is ∆T.

96 of 107

Appendix: Sampling Theorem

F(µ) is the Fourier transform of f(t) (observe the bandlimitedness). In the second row, we see the Fourier transform of its sampled version – we have earlier proved that it would be equal to the summation of infinitely many copies of F(µ). If ∆T is too large, then the different periods will merge and cause aliasing. Hence ∆T has to be “small enough”.

97 of 107

Appendix: Sampling Theorem

98 of 107

Appendix: Sampling Theorem

99 of 107

Appendix: Sampling Theorem

Refer to the previous slide. We have:

This is called the Whittaker-Shannon interpolation formula. In practice, it requires infinite summation, and hence methods like bilinear or bicubic interpolation are preferred. It can be implemented in truncated form, with some error

Sampled version of f

The inverse fourier transform of a rect is a sinc

100 of 107

Appendix: Sampling theorem in 2D

A 2D signal is band-limited if its Fourier transform is 0 outside a rectangle established by intervals [-µ_max,+µ_max] and [-ν_max,+ν_max].

The theorem states that such a continuous band-limited 2D signal can be accurately reconstructed from its equi-spaced spatial-domain samples provided the sampling rate is more than 2µ_max in the X-direction and more than 2ν_max in the Y-direction (equivalently the sampling period is less than 1/2µ_max and 1/2 ν_max).

100

102 of 107

Appendix: Sampling Theorem and Aliasing

If a signal is sampled at a rate less than the sampling rate, its analog version cannot be accurately reconstructed.

We get an erroneous version of that signal from its samples.

This phenomenon is called aliasing, because the (genuinely) higher frequency components “masquerade as” (i.e. pretend to be) lower frequency components.

Aliasing is illustrated in figure 4.6 of the book (also on slide 96 here).

A certain amount of aliasing is inevitable in practical signal or image acquisition, as in practice the signals can only be of finite length. And we know that time-limited or space-limited signals cannot be bandlimited.

102

103 of 107

Appendix: Sampling Theorem and Aliasing in 1D

103

Analog signal (to be reconstructed from its samples) shown in solid line:

f(t) = sin (πt),

Time period = 2 second

As per the sampling theorem, we need to sample this signal at a rate strictly greater than twice the largest frequency, i.e. greater than 2 x ½ = 1 sample per second. Equivalently the sampling rate must be strictly less than 1/(2 x ½) = 1 second. If you sample at exactly twice the largest frequency, i.e. 1 sample per second, you will collect values of sin (πt) for integer t, i.e. all zeroes. This will not allow for good reconstruction.

104 of 107

Appendix: Sampling Theorem and Aliasing in 2D (example 1)

104

Consider a perfect imaging system that can faithfully capture images of size 96 x 96 pixels (one pixel is the resolution or sampling period of the system in both X and Y directions). Such a system can faithfully capture a checkerboard pattern in which each square has size 1 pixel or more. See top row where the square width is 16 and 6 pixels respectively. But there is a problem if the size of the squares is less than 1 pixel: see bottom row where the squares have width 0.9174 and 0.4798 respectively. These produce misleading patterns which can masquerade as normal images.

105 of 107

Appendix: Sampling Theorem and Aliasing in 2D (example 2)

105

“original” shirt texture

Aliased shirt texture (see appearance of slanting patterns from top right to bottom keft near the midde shirt button)

106 of 107

Appendix: Sampling Theorem and Aliasing in 2D (example 3)

106

Aliasing can be controlled by a small amount of blurring – this is of course not without loss of information, but at least it does not produce fake textures.

1 of 107

2 of 107

3 of 107

4 of 107

5 of 107

6 of 107

7 of 107

8 of 107

9 of 107

10 of 107

11 of 107

12 of 107

13 of 107

14 of 107

15 of 107

16 of 107

17 of 107

18 of 107

19 of 107

20 of 107

21 of 107

22 of 107

23 of 107

24 of 107

25 of 107

26 of 107

27 of 107

28 of 107

29 of 107

30 of 107

31 of 107

32 of 107

33 of 107

34 of 107

35 of 107

36 of 107

37 of 107

38 of 107

39 of 107

40 of 107

41 of 107

42 of 107

43 of 107

44 of 107

45 of 107

46 of 107

47 of 107

48 of 107

49 of 107

50 of 107

51 of 107

52 of 107

53 of 107

54 of 107

55 of 107

56 of 107

57 of 107

58 of 107

59 of 107

60 of 107

61 of 107

62 of 107

63 of 107

64 of 107

65 of 107

66 of 107

67 of 107

68 of 107

69 of 107

70 of 107

71 of 107

72 of 107

73 of 107

74 of 107

75 of 107

76 of 107

77 of 107

78 of 107

79 of 107

80 of 107