Confession time. When it comes to writing about Fourier transforms on this blog, I'm a chicken. (For my readers not familiar with English slang, chicken in this context means a cowardly or fearful person.)
I'm a chicken because few signal processing topics cause as much confusion, and because few topics seem to inspire as many different, adamant, and sometimes contradictory opinions about the one true way to explain or implement Fourier transform concepts.
The way I see it, the confusion surrounding Fourier transforms stems from four interrelated issues. First, although most people who haven't studied signal processing in depth tend to think there is just one thing called "the Fourier transform," there are actually several. The way I count them, there are four basic types:
- The time domain (or spatial domain for image processing) and the frequency domain are both continuous, infinite domains. There is no explicit or implied periodicity in either domain. This is the one I call the Fourier transform.
- The time domain is continuous and the time-domain functions are periodic. The frequency domain is discrete. I call this the Fourier series.
- The time domain is discrete and infinite, and the frequency domain is continuous. In the frequency domain the transform is periodic. This is the discrete-time Fourier transform (DTFT).
- The time domain and the frequency domain are both discrete and finite. Although finite, the time and frequency domains are both implicitly periodic. This form is the discrete Fourier transform (DFT).
For mathematical analysis of linear time-invariant (or shift-invariant) systems, the Fourier transform and the DTFT are the most useful, depending on whether you are analyzing a continuous-time or discrete-time system. For numerical computation, the DFT is most useful.
The second factor that confuses people is the differing mathematical conventions that are used. The formulas vary in the use of scale factors, the sign in the complex exponential terms, the variables used, and whether sampling period and sampling frequency terms appear explicitly in the equations or are implicit.
The third confusing factor is that the relationships between the different Fourier transform forms can be quite nonintuitive, especially for people not accustomed to thinking in terms of periodic sequences, periodic symmetry, circular shifts, modular indexing, etc.
Finally, there's the confusion caused by the famous FFT acronym. FFT stands for fast Fourier transform, so people expect it to be the Fourier transform, whichever variation they learned. But really it's a fast way to compute one kind of Fourier transform, specifically the discrete Fourier transform.
What finally convinced me to try to write a post involving Fourier transforms was a question received by one of my coauthors of Digital Image Processing Using MATLAB. The questioner wanted to know why the Fourier transform of the image below looked so funny:
I've run out of blog writing time for today, though, so I'll just have to leave you with that teaser. I'll follow up with more about this image and its Fourier transform next time.
Get
the MATLAB code
Published with MATLAB® 7.9
Hey Steve.
Looking forward to your post.
Something I would like to see is how can you sample Sin beneath the sampling rate and observe the Aliasing effect.
By the way, may I suggest a feature for Matlab?
Would you apply the 100% image size into the “image” command and not only “imshow”?
You can see here:
http://www.mathworks.com/matlabcentral/newsreader/view_thread/266485
The trouble it gives me :-).
Thanks.
The image is deceptive, isn’t it? It seems to the eye a binary image with alternative black and white stripes, but in reality the whites stripes are non-uniform within themselves, with a peak at the centre. Morever, the different white stripes also vary in their amplitudes.
Or am I getting something wrong here?
Feel free to not include this if it destroys the surprise, but would the following code be correct to display the transform?
i = imread('http://blogs.mathworks.com/images/steve/2009/lines.png'); i = double(rgb2gray(i)); I = fftshift(fft(fftshift(double(i)))); I = abs(I); I = I./max(I(:)); imshow(I);I had to scale the magnitude because it was all white otherwise, but that was the only “magic” I had to apply. I always get confused when asked to plot a Fourier Transform, since it’s usually (but not *always*) implied that you just want the magnitude. Anyhow, I’m really excited to read the rest of your posts. Fourier transforms are one of those concepts that can never be explained enough.
Drazick—Thanks for the suggestions.
Steve,
Great post - it’s so useful to see things explained from basics - you know, the bit that gets said at the beginning of the lecture and then forgotten forever afterwards.
And Fourier-confusion is my dirty professional secret.
I’ll be looking out for your follow up blogs on this topic.
Thanks!
Looking forward to the answer! I’ll try plotting the FFT now to puzzle about the results :)
@francisco: I think what you want to do is to use fft2 to obtain 2D spectrum like this:
I think what looks surprising is that even though the image consists of single periodic sequence along NW-SE axis at high frequency, the spectrum shows a DC component and peaks at two spatial-frequencies of relatively small values. Looks like a case of aliasing (and perhaps harmonic generation due to squaring). btw, if someone wants to see a ‘point blimfark’ (i.e., aliasing when imaging something rotating) that Steve mentioned in one of his posts, I captured one with my low-speed SLR. http://shalin.wordpress.com/2009/01/11/an-interesting-occurence-of-undersampling/
btw, I have always been terribly puzzled by interplay of fft,ifft,fftshift and ifftshift. When we have real matrix (say,i), it seems that performing fft on ifftshift(i) or i directly gives the same result. Why is that the case?
Shalin—I do not know why you are calling fftshift or ifftshift before calling fft2. Note that fftshift performs a circular shift. One of the properties of the discrete Fourier transform is that a circular shift does not affect the magnitude of the DFT. So you would see exactly the same thing for the magnitude of fft(i) and the magnitude of fft(fftshift(i)).
Steve - When working with real and even functions, obtaining their real and even Fourier/inverse Fourier transforms, I have to use the idiom-
I had determined the above by trial and error. I have fallen into habit of using this idiom even though I didn’t understand how it worked. But the penny drops when you mention that fftshift and ifftshift implement circular shifts :) Here is the post about this on my blog: http://shalin.wordpress.com/2009/12/06/fftifft/. The circular shift doesn’t affect the magnitude but does affect the phase.
Shalin—Yes, I definitely plan to talk more about circular shifts and real/even properties of the DFT.
Shalin, Steve,
The issue with FFT is that is assumes omega=0 is the leftmost element of the array (omega goes from 0 to 2*pi), whereas we are used to studying the Fourier Transform with the origin in the middle (omega goes from -pi to pi). The same is true for the spatial domain — the way we look at it, we want the origin in the middle of the image, not in the top left corner. This is why you need to apply the FFTSHIFT to the image before applying FFT, and then need to FFTSHIFT the result again.
The solution is to either change the way we think about the Fourier Transform to match the definition of the FFT, or we do the sensible thing and define a function that does the shifting “behind the scenes” so we don’t have to think about it. In DIPimage we have a function FT that computes the Fourier Transform of an image using the “correct” coordinate systems. I think it helps didactically.
Cris—There’s no particular need to assume the spatial domain origin of an image is in the middle, and I won’t be teaching it that way here. And there’s no “issue” with the FFT function. It computes the discrete Fourier transform (DFT) according to the standard mathematical definition. As I said in an earlier post, the DFT is one of several related transforms. Much confusion is caused by people not being aware of this and assuming that the DFT is the same as “the” Fourier transform they learned in class.
Steve,
Sorry, yes, poor choice of words on my part. What I meant to say is that the definition of the FFT does not match the definition of the DFT we teach in class. If we were to properly teach the FFT as what it is, there would be less confusion. But we teach about the DFT, and let the students have a go with the FFT without warning them about the subtle difference. It’s not a problem with the FFT, but it is a problem with the software that we give the students to learn about image analysis. This is why I give my students a function FT that uses the FFT under the hood, but produces a result in line with their expectations. I don’t care that much about algorithmic implementation of an idea. The FFT is just a clever way of computing the DFT. This is why I don’t want to bother my students with what the FFT does. I just want them to understand the frequency domain and learn about filtering by examining what happens to the image in that domain.
Cris—Thanks for the clarification. I’m still a bit confused, though. You say that the “FFT is just a clever way of computing the DFT.” I agree with that. That’s what the MATLAB fft function does. But you also say that the FFT does not match the definition of the DFT we teach in class. What definition of DFT (discrete Fourier transform) do you use? For reference, MATLAB uses a definition of the DFT consistent with the formula on this page.
Steve,
I’m from a signal processing background, we use the DFT as the discrete version of the Fourier Transform, and hence use discrete frequencies in the range [-pi,pi). More of a “discretized Fourier transform” than a “discrete Fourier transform”?
The mathematical definition of the DFT as written in Wikipedia is what we were taught is the Fourier series representation of a discrete-time periodic signal. Fourier series decomposes the signal into N basis functions, numbered 0 through N-1. I know this is exactly the same thing as what I understand as the DFT, but it’s just a lot harder to explain filtering, sampling, etc. with that definition. Frequency components are a lot more natural to me than basis functions.
We all have a different view on things! Life truly is diverse. :)
This is unrelated, but interesting: I just learned that Gauss developed an algorithm very similar to the FFT!
Cris—I am also from a signal processing background. From your description, what you call the DFT I call the DTFT (discrete-time Fourier transform). That is the terminology used, for example, in Discrete-Time Signal Processing by Oppenheim and Schafer. Since Ron Schafer led the signal processing research group at Georgia Tech while I was there, it is natural that I am accustomed to that terminology. ;-)
This is the kind of thing I had mind when I said last week that I was afraid of writing about Fourier transforms. In my experience even very knowledgeable signal processing people can have a hard time discussing Fourier transforms with each other, at least until they figure out the terminology roadblocks.
Cris—I made a note to demonstrate how to do 1-D and 2-D DTFT plots.
Also, my viewpoint might not be all that different from yours. I tend to think in terms of the DTFT (what you call DFT) as well, and I think of frequency components rather than basis functions. But I view the DFT (what you call Fourier series of discrete-time periodic signal) as a discrete transform of a finite-length signal, and the resulting DFT coefficients simply sample the DTFT (as long as the DFT length is sufficiently long).
That’s funny. I learned following Oppenheim, Willsky & Young. That’s the same Oppenheim as in Oppenheim & Schafer. I’m guessing my choice of terminology has changed a bit over the years? :)
Cris—My initial introduction to signals and systems was also via Oppenheim, Willsky, and Young.