Interplay of fft, ifft, fftshift, and ifftshift in MATLAB

Contents

• The problem – Fourier transform of a real and even function.
• Let us take a concrete example – circular pupil
• First attempt at obtaining the point spread function.
• Second attempt – change the origin before ifft2
• The prescribed use
• Where was the confusion?
• The confusion resolved
• After-thoughts

Hello world,

I routinely simulate forward imaging in optical systems to understand how they function. One of the beautiful phenomena of optics is that a lens is a fourier transformer. It produces 2D fourier transform of field between planes located one focal length away on either side of itself. MATLAB is a ‘matrix laboratory’, so it implements discrete Fourier transform (DFT). The standard way of implementing DFT is FFT (Fast fourier transform) algorithm.

Although I have been using MATLAB since 6 years, I never quite grasped the interplay of fft, ifft, fftshift, and ifftshift functions found in MATLAB. Well, the penny dropped today after a comment from Steve on his blog that fftshift and ifftshift cause circular shifts and nothing more. By the way, Steve’s image processing blog is one of the most useful blogs and these days he is clarifying how Fourier transforms work in MATLAB. He just described that FFT algorithm inherently assumes that the function is periodic in both space and frequency. This occurs because whenever one of the domain (space or spatial-frequency) is sampled, its Fourier counterpart becomes periodic. Multiplication with impulse train (i.e., discretization) in one domain leads to convolution with impulse train in the other domain (i.e., periodicity). In DFT, both space and spatial-frequency are discrete (actually digital), and therefore, both domains are periodic.

The problem – Fourier transform of a real and even function.

Typical functions that occur in optics (e.g., pupil of an imaging system which is simply a circle) are real and even around origin (i.e., optical axis). We expect the Fourier transform or inverse Fourier transform of a real and even function to be real and even itself. However, I had trouble achieving this seemingly basic task. After good amount of experimentation, I had decided that the the idiom spec = fftshift(fft2(ifftshift(signal))) gives real and even spectrum for real and even signal. However, I understood why that is the case just recently.

Phase-space optics 2: Instantaneous correlation and pseudo-symbolic computation in MATLAB

As I noted in the previous post about FiO, I wish to share some of the simulations and basics about phase-space signal analysis. While doing simulations, I sometimes find it tricky to interpret analog equations and convert them to a sequence of programmed operations on discrete data. Also, one needs to be careful about not losing information when going from analog to discrete domain. This all-important step is known as A to D (analog to *discrete*) conversion and is responsible for Nyquist and Shannon being so famous.

I rely on two tools for my modeling needs: MATLAB (TM) most of the time and Mathematica (TM) sometimes. Mathematica is excellent for allowing quick visualization of the model due to its ability to work directly with symbols, whereas MATLAB is good for fast and maintainable implementation. It may be possible to write maintainable code in Mathematica but that requires going through steep learning curve, especially for folks who had learned classic ‘C’ as their first language. Is there a way of marrying symbolic ease of Mathematica with programming ease of MATLAB? I seem to have found one using anonymous functions. The code used in this post has been posted here on MATLAB file-exchange, which will be available within a day after moderation.

Let us consider the problem of calculating instantaneous correlation (which is a primary quantity required to compute and understand the Wigner distribution and Ambiguity function). For a signal it is defined as,

We know that the auto-correlation of a function is defined as Thus, we see that instantaneous correlation is in some sense the quantity which when integrated along gives us the auto-correlation function.

As we shall see in upcoming posts, Wigner distribution is the Fourier transform of the instantaneous correlation (Icorr) along the delay () dimension, whereas the Ambiguity function is Icorr’s Fourier transform along the normal () dimension. Therefore, there are two ways of interpreting what the instantaneous correlation is:

1. is the signal which is shifted by and The conjugate of the delayed version () when multiplied with the advanced version () gives us Icorr.
2. is the signal which is first dilated to give us This dilated signal’s conjugate is flipped to give The signal is advanced to give and the signal is delayed to give us

As is often the case – hindsight is 20/20. I could write above paragraph so clearly because I have already gone through simulating the problem with two tools. But things should be presented logically and not chronologically, so I explained the model first. Next are the simulations.

Phase-space optics 1: Introduction

As noted in the last post about Frontiers in Optics (FiO), one of the interesting things that will happen at the conference is the special symposium on phase-space optics. I have been intrigued by these ideas since sometime now. In upcoming posts, I will recount some basics of the phase-space optics. I hope those who are getting started in this area will find these posts informative or at least interesting, and those (particularly FiO attendees) who already have some insights will share them via comments.

I have found books by Leon Cohen [1] and a compilation of selected papers as part of SPIE milestone series [2] useful. Cohen has contributed greatly to the basic ideas of phase-space representations and showed that all phase-space distributions are special cases of a particular distribution – which has come to be known as the Cohen class. The phase-space representations have been used in optics (and signal processing in general) for analysis and representation of signals and systems in a way that matches our intuition more closely. However, once someone has been used to the Fourier tools for few years , the phase-space representations may not seem that intuitive, as happened in my case:-).

Joint distributions were first invented in quantum mechanics. They were used by Eugene Wigner [2, pp. 30] to represent probability of a particle possessing given position and given momentum. The distribution represented how particle may be `distributed’ as a function of position and momentum. Gabor [2, pp. 120] and Ville [2, pp. 149] introduced joint distributions to signal analysis to represent temporal signals in a way that matches human intuition. Until the works by Gabor and Ville, signal analysis was performed either in time-domain or in frequency-domain. However, when a human analyzes signal, he/she is usually interested in how signal’s frequency changes over time, e.g., how pitch of some one’s voice changes over time. The last sentence of the abstract of Ville’s paper reads “These notions of instantaneous frequency and of the instantaneous spectrum are introduced to furnish a firm theoretical basis for studies of frequency modulation, …, and in a general way,of all problems for which classical harmonic analysis furnishes a description which departs too far from physical reality.”

An interesting occurence of undersampling

Hello world,

I was checking out the burst mode of my digital SLR on a ceiling fan, and incidentally captured a nice demonstration of the effects of undersampling. Those who are from the digital signal processing area will find nothing new in the following video of a ceiling fan picking up speed.

(Updated video is now embedded.)

Those who are still thinking why the fan seems to start rotating in one direction and then reverses the direction, read on:

The world around us is analog (i.e. continuous in space-time),  but most of the storage equipment digital (i.e. discrete space-time + discrete values). In going from the reality of continuous space-time to stored bits, we perform analog-to-digital conversion.  The first step in this conversion is the discretization or ’sampling’ of the analog information. In this particular case, the DSLR is sampling the motion of the fan in time.  We are clearly losing lots of information in going from continous time to discrete time. But how much information can we lose, and still be able to observe the features of interest? i.e., what is the slowest sampling frequency that will let us have faithful representation of the signal that represents the event. Harry Nyquist discovered that for avoiding loss of information during discretization of continous signals, one needs to sample the value of the signal atleast two times faster than the rate at which the signal changes.  This idea came to be known as Nyquist sampling theorem.

When I switched on the fan and started capturing the burst of images using the camera, the camera was fast enough to capture the movement of the fan. Or I should say, the fan was slow enough such that its movement could be captured faithfully by the camera’s shutter speed. Nyquist theorem was satisfied and we meausred ‘true’ direction in which the fan was rotating.  However, when the fan picked up the speed, the camera shutter was no longer fast enough. More specifically, at some point after the fan was turned on, the fan’s speed of rotation (measured in rotations per second) exceeded twice the value of the speed  of camera’s burst mode (measured in number of exposures per second). This led to ‘false’ measurement of the direction in which the fan was rotating. This effect is called aliasing of information. This was incidental but intriguing occurrence of undersampling.

Why use sine-waves to represent everything?

This issue was triggered off by my friend Bhavesh who asked:

why do we always want to represent any wave/ any signal in terms of sine wave?

We represent the whole spectrum in terms of sine wave frequencies.. all our modulation techniques are based on sine wave..

why sine wave is so prominent? What is its basic charactistic that helps us?

———————————————————————————————

The answer to the key question (what is the basic characteristic that helps us?) is this:

We usually model our systems as LTI (Linear Time invariant) system.

You get the general picture if you think about basic equations which arise in *mathematical representation of real life*.

Whenever we model something, we assume (or we are taught to assume) that system is linear (so that we can easily predict the output for unknown input from the output for known inputs) and time-invariant (who likes system that will work one way today and the other way tomorrow!!).

Now here comes linear algebra. Consider input and output to be signals and the system to be an operator on those signals.
Spectral theory of linear algebra says that if we express signals as linear combination of eigenfunctions of an operator;  analysis, interpretation and computation of determining output for a given input becomes very simple.

(This is because eigenfunction is a signal that remains unchanged by an operator except for some scaling or pahse-shift (or equivalent effects in other domains).)

It just so happens that complex exponentials and sinusoids are eigenfunctions of LTI systems (which are our models of real situations). We want to express our signals in terms of these eignefunctions (sinusoids) to simplify analysis, interpretation and computation.

By similar argument, if you choose some other model for your system or phenomenon sine waves are not that suitable.