Power spectrum of an image

AJW picture AJW · Dec 9, 2013 · Viewed 16.5k times · Source

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.

Answer

kamjagin picture kamjagin · Dec 10, 2013

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.