How to plot a Power Spectrum (log log plot) for an image?

13 views (last 30 days)
Given Problem statement:- to detect blur image
Hint: The power spectrum contains the frequency information of the image, and the slope gives a measure of image blur. A higher slope indicates more lower frequency components, and hence more blur.
what i have been doing is that...
i=imread('actual.png');
i=rgb2gray(i);
i=imresize(i,[256,256]);
figure,imshow(i);
F=fft2(i);
F=fftshift(F);
p=(abs(F).^2);
p=log10(p);
fs=1;
f=fs/2*logspace(-10,10,256);
plot(f,p);
xlabel('frequency(Hz)');
ylabel('power');
title('{\bf LogLogPlot}');
i want to get a single line plot, while it is showing multicolored line plots. Is it because the image is gray image not binary one? or am I just not working in the right direction? Plz help. thanku

Answers (2)

Image Analyst
Image Analyst on 15 Apr 2015
Use imshow() to see the spectrum.
I'm not sure who wrote your hint, but it's misleading. I could be wrong but with a Ph.D. in optics and 37 years of developing imaging algorithms I hope I know at least 0.001% of the imaging concepts out there and a pretty darn good intuitive feel for fourier domain concepts. The first clause of the hint talks about the power spectrum, but what it says after that applies to the SPATIAL domain, not the spectral domain. The slope does get flatter in the spatial domain, like when sharp steep edges get blurred out. It's not right to say that flat or steep slopes in the spectral domain necessarily mean anything. What's more important is the values out in the higher frequencies of the spectral domain. For example, if you have a periodic texture in the spatial domain, you would have steep slopes in some places in the spectral domain because the periodic structure in the spatial domain will give periodic spikes in the spectral domain. Now the author may be thinking of MTF curves where the slope at mid-level frequencies may (or may not) get steeper as the blurring gets worse. Again, it's the value (height) or the amount of power in those frequencies that is the crucial thing, not the slope . As the Full Width Half Max moves inwards as the image becomes more blurred, the slope may get steeper (higher) but you may have a system that has a steeper slope out at a higher frequency and has less blur. I could draw some MTF curves to illustrate if you want.

pfb
pfb on 14 Apr 2015
Edited: pfb on 14 Apr 2015
I think that's because f and p are matrices. In that case all of its columns are plotted. The lines you see are
plot(f(:,1),p(:,1)),
plot(f(:,2),p(:,2)),
...
Try
plot(f(:), p(:))
I guess you probably have to sort the entries of f (and p according to f), in order not to have a very garbled line. Or you could go
plot(f(:), p(:),'.')
that should work, if the abscissae are in the right places.
  2 Comments
dv
dv on 15 Apr 2015
Edited: dv on 15 Apr 2015
Here, as we have size(p)=256 256 and size(f)=1 256 r u saying that the anomaly is because order of p is 256 x 256 not 1 x 256? should i convert p into a vector? plz tell me where am i going wrong. thank you
pfb
pfb on 15 Apr 2015
I'm not an expert in image processing, but what you are doing is plotting a matrix (resulting from a fft) versus a vector with the same length as the column. That is the same as plotting each column of the matrix separately (plot works along columns).
I'm not sure that all of the columns in p should have the same set of abscissae. Maybe so. If this is the case, and you want one plot, you could do the following
% this produces a matrix whose columns are copies of f
ff = f'*ones(size(p,2));
% this "unravels" the matrices into column vectors, and plots them
plot(ff(:),pp(:),'.');
However, for some reason I would think that your "abscissae matrix" should not be made of identical columns. I think that it should be formed from the "wavevectors" corresponding to the modes in the fft. Something like
L = length(p);
x = cos(2*pi/L*(0:(L-1)));
[xx yy] = meshgrid(x,x);
x = xx.^2+yy.^2; % perhaps under sqrt
x = fftshift(f);
plot(x(:),p(:),'.');
or since you have a logarithm
semilogx(x(:),p(:),'.')
But again, I'm no expert in this field.

Sign in to comment.

Categories

Find more on Get Started with Signal Processing Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!