fftshift of an image

207 views (last 30 days)
Paul T John
Paul T John on 17 Aug 2011
Commented: Walter Roberson on 14 Oct 2023
Hi all,
I'm trying out to find the wavelength of the wave seen in this image: http://www.flickr.com/photos/66468892@N07/6052890522/
I did the fftshift of the red channel [ im(:,:,1) ] and got this image: http://www.flickr.com/photos/66468892@N07/6052339583/
I did the fftshift of the grayscale [ rgb2gray ] and got this image: http://www.flickr.com/photos/66468892@N07/6052365079/
I think the first pair of lateral peaks on either side of the DC peak, in second and third images, corresponds to the grid lines (a grid pattern can be seen in the first image). How could I correlate the co-ordinates of the next pair of lower peaks to the image? Does those peaks correspond to some characteristic of the wave seen in the first image?
Thanks.

Accepted Answer

Rick Rosson
Rick Rosson on 18 Aug 2011
In order to use the DFT to determine the wavelength of the pattern in the image, you will need one key piece of information: some way to relate the size of the pixels in the image to an actual length in a physical unit of measure (for example: inches, centimeters, or millimeters). More specifically, you will need to know (or at least assume) a value for the number of pixels per unit of physical distance.
I will assume for now that you would like to know the wavelength in units of centimeters. Then you would need a value for the number of pixels per centimeter in the x-direction, and likewise in the y-direction (these two values could be the same or different, depending on the scaling of the image).
Let's call these parameters the pixel sampling rates for the X-direction and the Y-direction, respectively, and use the MATLAB variables Fs_x and Fs_y to represent them. The units of measure for these variables are pixels per centimeter. I will assume the following values:
Fs_x = 50; % pixels per centimeter
Fs_y = 40;
We can then compute the pixel sizes in each direction as the reciprocal of the pixel sampling rates:
dx = 1/Fs_x; % centimeters per pixel
dy = 1/Fs_y;
We can then compute a linear distance scale for the X- and Y-axes of the cropped image:
[ M, N, ~ ] = size(cimg); % pixels
x = dx*(0:N-1)'; % centimeters
y = dy*(0:M-1)';
Next, we compute the frequency increments for both X and Y:
dFx = Fs_x/N; % cycles per centimeter
dFy = Fs_y/M;
Finally, we can create a frequency domain array for both X and Y:
Fx = (-Fs_x/2:dFx:Fs_x/2-dFx)'; % cycles per centimeter
Fy = (-Fs_y/2:dFy:Fs_y/2-dFy)';
Almost there (whew!). The last piece is to use fft2 and fftshift to compute the 2D DFT of the image, and then to display the resulting spectral images using Fx and Fy to represent the X- and Y- axes, respectively. From there, you should be able to find the peak values in the spectrum in both directions. The wavelength in each direction is then the reciprocal of the frequencies associated with the peak values.
One other thing: it probably makes sense to remove the DC component of the cropped image prior to computing the DFT. That way, the peaks in the spectrum will show up much more clearly (since they will not have to compete for attention with the DC value).
HTH.
Best,
Rick
  1 Comment
Paul T John
Paul T John on 22 Aug 2011
Thanks a lot Rick. That helped very much.
However, I would like to get two doubts cleared:
1) I did the 2D FFT and cleared the DC component. [Image is seen here: http://www.flickr.com/photos/66468892@N07/6062522739/
Now, how can I relate the peaks to the colour waves and also to the grids? Which peak corresponds to which?
2) Also, once I find out the Fx and Fy of the peak and get the wavelength in each direction, how can I find the resultant wavelength of the propagating wave?

Sign in to comment.

More Answers (5)

Paul T John
Paul T John on 18 Aug 2011
I did 2D-FFT prior to fftshift. Also, the FFT-image is actually a 2D spectrum, and it appears as 1D because I used 'view(0,0)' command. I am adding the code below:
close all;
clear all;
clc
img = imread('f75p0.png'); %importing pixel information
figure(1); imshow(img);
cimg = img(1:485,1:216,1:3); %cropping the pixels
figure(2); imshow(cimg);
imgred = cimg(:,:,1); %just the red channel
imggrey = rgb2gray(cimg); %greyscale of the cropped image
figure(3); imshow(imgred);
dred = double(imgred);
dgrey = double (imggrey);
figure(4); mesh(dred);
imgfftred = fftshift(fft2(dred));
figure(5); surf(abs(imgfftred));
view(0,0);
shading interp;
imgfftgrey = fftshift(fft2(dgrey));
figure(6); surf(abs(imgfftgrey));
view(0,0);
shading interp;
I hope now you would be able to reproduce the code. Sorry for the confusion.
  2 Comments
Mehri Mehrnia
Mehri Mehrnia on 14 Oct 2023
why fftshift is necessary?
Walter Roberson
Walter Roberson on 14 Oct 2023
fft() returns back the 2-sided fft. By convention, the order of frequency bins returns back from fft() of real input signals is like
0, 1, 2, 3, 4, ... n-1, n, conj(-n), conj(-(n-1)), ... conj(-(4)), conj(-(3)), conj(-(2)), conj(-1)
but people often prefer to visualize it in the order
conj(-(n-1)), ... conj(-(4)), conj(-(3)), conj(-(2)), conj(-1), 0, 1, 2, 3, 4, ... n-1, n
The change is not mathematically significant: it is just easier for people to look at the plots.

Sign in to comment.


Joshua Canfield
Joshua Canfield on 8 May 2019
Edited: Joshua Canfield on 8 May 2019
How do I plot the 2D DFT against Fx or Fy to obtain the usefull spacial frequency information
Could I sum columns and rows and then plot against Fx or Fy respetively?

Rick Rosson
Rick Rosson on 18 Aug 2011
The result of the DFT of an image should also be an image, but the graphs are showing a 1D spectrum. How did you compute the DFT of the image? Prior to calling fftshift, did you use the fft or fft2 function (or some other method)? Can you show a more complete set of MATLAB code that is reproducible by others?
  1 Comment
Paul T John
Paul T John on 18 Aug 2011
Hi Rick, I'm attaching the code (as another answer). Please see below.

Sign in to comment.


Rick Rosson
Rick Rosson on 18 Oct 2011
Hi Joseph,
Walter is correct. A single line or line segment, as in the linked image, contains an infinite number of frequencies (or wavelengths).
To expand a bit on this concept:
The DFT of an image is a transform from the space domain (in mm or cm) to the spatial frequency domain (in radians per cm or radians per mm), also known as the wave number. If we use x and y to denote the space domain, then the spatial frequency would be kx and ky respectively. A line segment or rectangle in space does transform into a sinc function in the spatial frequency domain.
It is a bit misleading or confusing to characterize it as sinc( x ) versus x. Rather, it should be sinc( kx ) versus kx. This distinction is crucial in this discussion, because there is a one-to-one relationship between kx and the x-component of the wavelength (and likewise for y). That is,
Lx = 2*pi/kx;
Ly = 2*pi/ky;
So it is actually a function of wavelength, and it includes an infinite number of wavelengths, as Walter points out.
HTH.
Rick
  1 Comment
Gova ReDDy
Gova ReDDy on 20 Oct 2011
Thankyou walter and Rick...

Sign in to comment.


dustin Gallegos
dustin Gallegos on 1 Nov 2011
Hi there - Can you help me do something similar: I am trying to determine the wavelengths of the different colors in a given image. Ultimately, I would like to plot lambda(nm) vs Intensity from the picture.
Thanks !
  1 Comment
Walter Roberson
Walter Roberson on 1 Nov 2011
I suggest you start a new Question with that, and that you post a sample image; See http://www.mathworks.com/matlabcentral/answers/7924-where-can-i-upload-images-and-files-for-use-on-matlab-answers
Also maybe you would get the information you need from http://www.mathworks.com/matlabcentral/answers/17011-color-wave-length-and-hue

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!