When I dipped my toe into the Fourier transform waters last week, the resulting comments and e-mail indicated there is a lot of interest in tutorial material and explanations. I'm willing, but I'll have to think about how to do it. I haven't taught that material since my professor days, and that was many, many moons ago.
But I can certainly start by following up on the teaser example from last week. A reader of Digital Image Processing Using MATLAB wanted to know why the Fourier transform of the image below looked so "funny." (For the moment I'm going to use the term Fourier transform fairly loosely as many people do.)
url = 'http://blogs.mathworks.com/images/steve/2009/lines.png'; f = imread(url); imshow(f)
I didn't show what the Fourier transform looked like, though, so I suppose it wasn't very clear what problem I was talking about.
The book recommends this procedure for visualizing the Fourier transform of an image:
F = fft2(f); imshow(log(abs(fftshift(F)) + 1), )
The puzzle is why does the Fourier transform look so complicated? The input image, after all, contains only a simple sinusoidal pattern.
The answer is related to the fact that what we're actually computing when we call fft2 is the two-dimensional discrete Fourier transform (DFT). The DFT has an implicit periodicity in both the spatial and the frequency domains. We don't often don't think about the periodicity in either domain, both because we only need one period for computation and because only one period is usually shown.
Let's see what happens when you take the original image and replicate it a few times.
I4 = repmat(I, 2, 2); imshow(I4)
Suddenly the simple sinusoidal pattern doesn't look quite so simple. What you get here (besides eye strain from looking at this image) is a pattern of small vertical and horizontal discontinuities. These discontinuities are producing the rectangular pattern in the frequency-domain visualization above. That's because any kind of discontinuity has frequency-domain energy across a broad range of frequencies.
I see this kind of thing in Fourier transform plots all the time, although often the effect is more subtle. Here's another example.
g = imread('rice.png'); imshow(g)
G = fft2(g); imshow(log(abs(fftshift(G)) + 1), )
What's that vertical line in the frequency domain? There are no horizontal features in the image to account for it. At least, no horizontal features are apparent until you display the periodically replicated image.
imshow(repmat(g, 2, 2))
Since the original image is darker at the bottom than at the top, there is a strong horizontal discontinuity at the periodic boundary. That's what causes the vertical line to appear in the frequency domain.
One way to remove these frequency-domain effects is to taper the original image values toward 0 at the image boundaries. That removes the periodic discontinuities. The tapering can be done by elementwise multiplication of the original image by a matrix equal to 1 in the center and trending toward 0 at the edges. This is called windowing.
Let's try that here using a raised cosine function.
[M, N] = size(f); w1 = cos(linspace(-pi/2, pi/2, M)); w2 = cos(linspace(-pi/2, pi/2, N)); w = w1' * w2; f2 = im2uint8(im2single(f) .* w); imshow(f2)
F2 = fft2(f2); imshow(log(abs(fftshift(F2)) + 1), )
Although there are still some visible artifacts, this frequency-domain plot is much closer to showing what we might have expected: the Fourier transform of a simple sinusoidal pattern.
Here's the same procedure applied to the rice image.
[M, N] = size(g); w1 = cos(linspace(-pi/2, pi/2, M)); w2 = cos(linspace(-pi/2, pi/2, N)); w = w1' * w2; g2 = im2uint8(im2single(g) .* w); imshow(g2)
G2 = fft2(g2); imshow(log(abs(fftshift(G2)) + 1), )
That vertical line in the previous frequency-domain plot is no longer visible.
Since people seem definitely interested in understanding Fourier transforms better, I'll work on putting together more material for the blog. You can help me by continuing to post your comments, including the things that puzzle you.
Get the MATLAB code
Published with MATLAB® 7.9
32 CommentsOldest to Newest
Steve, maybe you can address this in next episode, but what are the artifacts reamining in the diagonal sinusoid windowed FFT plot? Is this a result of the windowning function?? This is an excellent subject for a blog…keep em coming!
John—I haven’t looked at it very closely, but I can give you a “hand-wavy” explanation. One source of artifacts for apparently “ideal” functions like sinusoids is that we’re always looking at a signal over some finite-extent domain. A real sinusoid goes on forever; anything less than that and you’ll see artifacts.
Can the artifacts of the ‘ideal’ functions be thought of as the convolution of the window freq domain with the ideal freq domain?
Ketan—It depends on how you look at it. You can also think of them as being artifacts remaining from the original frequency-domain view that have been only partially removed.
I wanted to know why you calculated the fourier transform as you did. This is waaay off from my guess at it last post, and I want to know what I’m missing.
More specifically, why do you add 1 to the magnitude of the fourier transform, and why do you take the log afterwards? Should we always plot Fourier transforms on a log scale?
Francisco—Taking the log is just a way to make small details in the magnitude of the Fourier transform more visible. Adding 1 first is to avoid taking the log of 0.
the F-T image seems to has AWGN noise. Maybe using Huffman compression we can compress the F-T image more efficiently than we can in spacial domain. So my finding is that the amount of information that can represented by an image on spacial domain has insight a surplus amount of information. So when transform an image in another domain we achieve get over of this surplus information. Is that true? But now i am thinking that if we transform again the image in frequency domain assuming that this image is in spacial domain maybe we can achieve minimize information. Maybe we can break the theorem of Shanon. Maybe i opened a new section of conversation with title:lossless image compression using transforms…i am intresting a lot about it.
Dimitris—No, the Fourier transform image does not contain noise. Image compression is off-topic, and I have no plans to write about it.
JohnA, the cross artifacts are due to the windowing. Because Steve used a cos window in each direction, you’re going to get that. It’s kind of like multiplying by a rect but smoother, so Fourier Transform is convolved with the spectrum of the window, which is sort of a smooth-sided rect function. This will give you the spectrum of the window at every spike. Now you know a sharp-sided rect function would give you a 2D sinc function (which looks like a lumpy cross) but since the edges are smoothed by the cosine, you get kind of that shape but smoother, hence the smooth cross. If he had uses a circularly symmetric window, the artifacts would be like a smoothed Bessel function (you get the Bessel function for a sharp sided perfectly circular disk). This would appear more like rings than a cross. That should explain it but it won’t really sink in until you see the equations:
Let’s say that the FT(I) = F (fourier transform of I, let’s call that F). Now multiply I in the spatial domain by a window. FT(I * W) = FT(I) ** FT(W) because multiplication in the spatial domain equals convolution in the Fourier domain. (Notation: * = multiplication, and ** = convolution) This is why you see the spectrum of the window everywhere. It IS EVERYWHERE actually but you really only notice it when the value is high, such as at a spike.
You have used windowing here. How would the result differ if we use zero padding. Will we achieve similar results with zero padding in these cases? If results are same which method would you suggest?
Rakesh—The only effect of zero-padding is to sample the discrete-time Fourier transform. It won’t change the overall frequency-domain shape, and at this resolution I wouldn’t expect to much visual difference at all. Zero-padding is one of the topics I will definitely discuss in a later post.
Great, great post!!
I am trying to push this a little bit further to filter in the fourier domain. I am using these folk’s elliptical butterworth for that:
This is the first time I work with images so I am still a bit confused by formats. I get an error message when I try to use ifft2. Any suggestions/recommendations?
Thank you Steve
Here is my code:
url = 'https://5243270857495011345-a-1802744773732722657-s-sites.googlegroups.com/site/arnigolidatapost/home/images/Image004_crop.png?attachauth=ANoY7cpXS_3ULcnNIvsBkDZjdNDwP9q0e7cc-2IvA1QuWRklvZ_KD4h3kvwCWP6XwqTEtjyTqVYQyxMJiqTPZRn2Spt29w_0grqE41SgyZmh4W4CMpMSzTq6ycsizDHNj5HNXrShYTh9t_V16HDsGcnOltRdgLRurrnrauyK97dutLR9sCEfmm9qM1HI-qqbyO1QdMiWahC4KcBufbbY4CcDfzU2Kpkx1NA49zgGAmD-HrpUDSqBDJw%3D&attredirects=0'; f= imread(url); f=rgb2gray(f); imshow(f); F = fft2(f); imshow(log(abs(fftshift(F)) + 1), ); ef1=(1-EFILTER(sz, 0.99, 0.01,35)); ef2=(1-EFILTER(sz, 0.99, 0.01,35,pi/2)); eF2=eF.*ef2; imshow(eF2); invf=iff2(eF2);
And this is the message
??? Undefined function or method ‘iff2′ for input arguments of type ‘double’.
Matteo—If you read your error message very carefully, you might achieve enlightenment. ;-) Hint: There is no MATLAB function called “iff2″.
Oh, great! Should wear an an arlequin costum rather than Santa nex week, eh? ;-) Sorry for the trouble.
But now that I reflected on my code, I notice that I am not applying the filters properly – i.e. I am applying them on F wher I intended to apply them to the Fourier with the low frequencies in the centre. Any advice on how to adjust my code? This is not covered in too great detail in the user manual. Thank you again Steve. This blog is awesome.
Matteo—Possibly later in my Fourier transform series you will see some tips on frequency-based filtering.
Will wait patiently and eagerly. I hope you will do so. Tricky topic. Keep up, this is very exiting. Thank you
Your article is interesting in that it shows the benefit of windowing. But how about the different types of windows? Does it matter?
Matt—For power spectrum estimation, the choice of window is definitely important. Different windows are often compared in the frequency domain, where the key characteristics are the width of the main lobe and the height of the side lobes.
Spectrum estimation, though, isn’t what people are usually doing with images. For frequency-domain visualization, which is what I was showing here, I doubt that the choice of window function matters that much.
I have another question for you. I do some work with filtering of potential filed data, particularly gravity. In Matlab, I treat these maps as images. In the past to separate regional (low frequency) from local (high frequency) effects I used with good results either upward continuation, or best-fit polynomial trend analysis. But I’ve seen many authors use bandpass filters and I decided to tackle it. In a recent paper I’ve seen the wave numbers converted to kilometer scale, which is very intuitive because then you can relate the size of the bandpass filter to that of the geological features you are trying to filter out. Any hint as to how this could be done? I have seen next to zero literature on the math of it for matlab or even in general. Thank you, Matteo
Matteo—I’ll see if I can get this far in my series. At the least I’ll try to show how to do a frequency-domain visualization with a kilometer scale. If you have seen such a visualization that you like, in a paper maybe, can you send it to me? The frequency-domain processing chapter of Digital Image Processing Using MATLAB has a lot of information about bandpass filtering.
hi,can you tell me how I can get the fequency from the frequency domain?
Lillyne—Your question is kind of vague, but perhaps you are asking about the scaling of the frequency axis for different kinds of Fourier transforms. Please take a look at the posts in my Fourier transform series for more information.
Reading your problem statement and what we observe after running the code seems almost like an detective story with lots of mysterious effects but very nice explanations in the end.
One question that came to my mind is related to the “rectangular pattern in the frequency-domain visualization” for the first pattern. Is there any general reason for the appearance of rectangular pattern ?
Jayram—The quick answer is that is that you’re computing the Fourier transform of something that’s nonzero over a rectangular domain. Does that help?
Hi Steve,I know F(0,0) is often called the "DC component" and with zero frequency,but what is the frequency of F(m,n)? After F(0,0) is moved to the center with commond fftshift,then what is the requency distribution now? I'm lost.
Can you give me any suggestion as to how I could create a sine wave at an angle like in your first figure? I could easily figure out how to create a vertical 2D sinusoid, but introducing an angle?
frq=0.02; x=0:1:499; [xx,yy]=meshgrid(x,x); csw=cos(12*pi*frq*xx); imagesc(x,x,csw); colormap(bone); axis equal axis tight axis off
please ignore my previous question. It was easy to figure out with a little more thinking, and I should not have asked. This worked well:
% a 2D sinusoidal at 45 deg frq=0.02; x=0:1:511; [xx,yy]=meshgrid(x,x); csw45=cos(12*pi*frq*xx+12*pi*frq*yy);
But I have a more interesting question. Suppose I now create a second sinusoid at right angle from the one above one and want to combine the two. My guess is I could just Fourier transform both, sum those up, then just inverse Fourier transform the result. What do you think?
Thank you so much,
Arni—Since the Fourier transform is linear, the procedure you describe is equivalent to just adding the signals together.
Hi again Steve
Thanks for your comments, very much appreciated.
Thanks for the tutorial. I have a doubt about the illustration of rice image. When rice image is replicated there is a discontinuity in horizontal direction then why is there a a vertical line in its fft? How can directional information be infered from fft of an image. Can you please comment on this.
Rakesh—Good question about the horizontal/vertical line issue. I might do a blog post on it.