6

I have begun (a small project) to calculate the power spectrum of an image in the frequency domain.

So, what I have till now is the following:

%// close all; clear all; %// not generally appreciated
img   = imread('ajw_pic.jpg','jpg'); % it is a color image
img = rgb2gray(img); %// change to gray
psd = 10*log10(abs(fftshift(fft2(img))).^2 );
figure(2); clf
mesh(psd)

So far it looks good; I get the mesh plot which resembles the spectra I see in various academic papers.

However, what I am looking for is a graph plot of this power spectra versus the frequency, and I am not entirely sure how to get this frequency vector. I could do for instance:

N=400;        %// the image is 400 x 400
f=-N/2:N/2-1; %// possible frequencies?

but I am not convinced this is entirely correct as this gives rise to negative frequencies.

I'd really appreciate if someone could point me in the right direction to plot log frequency versus the power spectrum.

Adriaan
  • 17,741
  • 7
  • 42
  • 75
AJW
  • 5,569
  • 10
  • 44
  • 57

2 Answers2

2

fft "splits" a signal into frequency "bins", where the maximum frequency you can observe is the Nyquist frequency or half your sampling frequency. This means that for:

Y = fft(X,N); %  (1D case)

the frequencies corresponding to fft-values in Y(1:N/2+1) will be: 

f = [fs/2*linspace(0,1,N/2+1)]; % where fs is your sampling frequency

The other half of Y is just mirrored and comes from the intrinsics of the Fourier transform. If you don't want to understand it fully, I would say that there is no need to bother with it other than what you can find on wikipedia. But for the sake of curiosity you can check out a nice intuitive illustration of the origin of positive and negative frequencies: https://dsp.stackexchange.com/questions/431/what-is-the-physical-significance-of-negative-frequencies/449#449.

Some key differences for your 2D-image case:

  1. Using fftshift you've moved the 0-frequencies to the center of the matrix, rather than having them at the edges as in the above 1D example. So you would actually get f = fs/2 * linspace(-1,1,N) (once again, don't mind the negative frequencies)

  2. The next issue is to obtain the sampling frequency. Spatial frequencies are typically measured in [mm^-1], so in order to obtain it you actually need to know the physical distance between pixel centres (pixel pitch). But you could of course consider the spatial frequencies in [pixels^-1] in which case you are ready to go.

Community
  • 1
  • 1
kamjagin
  • 3,614
  • 1
  • 22
  • 24
  • Thank you very much for your detailed reply. I do seem to understand the negative frequencies bit - but, what I am still not entirely sure on is what should my freq's be after the fft. So, as you have explained, I could use `1/400` (given 400is my px), but I do not have an `fs` since this is an image? Also, when I fft my 400x400 image, and calculate power, the result is also a 400 x 400 matrix - but, I guess what I want is something of size 1x400 - so that I can plot frequencies and the power spectrum. Or am I being stupid here? – AJW Dec 10 '13 at 00:24
  • Well your fs is actually 1 in pixels or 1/(pixel pitch) in mm. If you don't have the pixel pitch you can still use fs=1 (just mind the unit - your nyquist frequency is then half an inverted pixel). You should also consider that you can have different frequencies in x and y direction (e.g. non-quadratic pixels) - that is why you have a 400x400 matrix. So create an fx and fy frequency vector (they are the same since NX=NY=400, and pixel size is unknown). How to visualize it is a different question though (that I don't have a good answer for other than labeling the axis with fx and fy as tics) – kamjagin Dec 10 '13 at 00:52
  • I am trying to do what you have suggested: So, I do: `fs=1;` `f = fs/2 * linspace(-1,1,N); %where N is 400` - this works out to `-0.5 to +0.5` - Does this make sense? Is this a Hz value or radians/sec? :( – AJW Dec 10 '13 at 12:56
  • Just another small thought: I am also using this resource: `http://redwood.berkeley.edu/bruno/npb261b/lab2/lab2.html` and the freq here seems to be `0:N/2` They somehow seem to calculate a `rotavg` - but the function for sure is not from Matlab. Wondering if you had any thoughts. Thanks. – AJW Dec 10 '13 at 13:03
  • Take a look here to get some input on rotational averages: http://www.sbirc.ed.ac.uk/cyril/download/DTP_Fourier_analyses.pdf. The freq range that you mention is consistent with my reply above (i.e. fs=1 and freqs going from 0:N/2 * fs, plus mirrored negative freqs) – kamjagin Dec 10 '13 at 15:02
  • The notes you have pointed out to are simply the best I have seen - thanks so much for this. I have used the method on page 12-15 to obtain the desired plot. Thanks a lot again - I will keep reading on this to learn more! – AJW Dec 11 '13 at 19:32
1

To plot the power spectra versus frequency of the image, one can use a process called 'radial averaging'. This calculates the average value of pixels that are a certain radial distance from the center of the image. Applying this to a power spectral density matrix results in a line plot of power versus frequency.

For more info and a a sample MATLAB code: http://www.mathworks.com/matlabcentral/fileexchange/46468-radialavg-zip/content/radialavg.m

Adriaan
  • 17,741
  • 7
  • 42
  • 75
M. Bedross
  • 109
  • 1
  • 2
  • 13
  • I removed the first paragraph, as this question indeed has a solution, hence the accepted answer. Also, everything is kept on this site, thus future reference is always the case. As for your answer itself: Please include all relevant code from the link (if possible) into the answer itself, as links can rot away. If you are not the author of that code it's customary to add a reference/citation. – Adriaan Jan 11 '16 at 23:55