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.
The power spectrum of a signal indicates the relative magnitudes of the frequency components that combine to make up the signal. The data used to determine the power spectrum must reflect sufficient excitation in the signal.
Power spectrum (PS) of biological time series (of an electroencephalogram recording, for instance) often shows a relationship of decreasing power as a function of frequency (f) according to the general equation: PS(f) = ψ × f-α (Norena et al., 2010).
Computations Using the FFT The power spectrum shows power as the mean squared amplitude at each frequency line but includes no phase information. Because the power spectrum loses phase information, you may want to use the FFT to view both the frequency and the phase information of a signal.
Power spectral density specifies the power levels of the frequency components present in a signal. It is denoted as PSD inshort. The PSD specifies the power of various frequencies present in the signal and we can determine the range of power over which the signal frequencies are operating at.
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:
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)
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With