Power Spectral Density Figure

Hi all,

I've spent the better part of two days flailing at this, so I figured it's time to ask. I have a set of images, and I need to calculate the average power at several saptial frequencies such I can derive an average for the entire set of images. The article I am trying to follow said "The spatial frequency content of images used in the experiment was analyzed by computing a periodogram with overlapping Hamming windows (a power spectral density estimate using the Welch’s method in the MATLAB image processing toolbox, averaged over horizontal and vertical dimensions.) Their resulting figure is (note, thre are 4 sets of images in this figure, each a unique color):

I know there are myriad threads on similar questions, but I just can't seem to piece together the process for getting this type of plot. Any help is appreciated.

[Bjorn Gustavsson]

It is difficult to figure out what's going on. Since the plot you show have exactly 10 points between 1 and 10 (presumably in oscillations/image-size), this is most confusing. Considering that there are 52 unique wave-vector amplitudes between 1 and 10, it is impossible to guess what they've done. Either only looked at the purely horizontal and vertical spectral components, or they've binned/averaged the power into 10 bins? Exactly what result they get also ought to depend on the window-length - that determines the absolute value of the smallest wave-number (longest wave-length) they can determine.

[David March]

Bjorn,

Instead of re-creating exactly what they have done, do you have a suggestion on how I can do something analagous? Plot the power by the spatial frequencies?

[Bjorn Gustavsson]

Perhaps you could do something like this (on one image...):

Im = peaks(128)+1;      
Im = Im(1:111,1:end);   
fI = fft2(Im);
asfI =  (abs(fftshift(fI)));
kx = [-64:63];
ky = [-55:55];
[kx,ky] = meshgrid(kx,ky);
plot((kx(:).^2+ky(:).^2).^.5,asfI(:),'.')
set(gca,'YScale','log')
set(gca,'XScale','log')

Here I use the 2-D fft because from your description I havent gotten any reason to use any piece-wise moving-window spectral calculations. (this also means that this lets us resolve longer-wavelength variations (smalle k)). Also I haven't bothered to add up the contributions from each quadrant, that should be done to get the power at each .

[David March]

Thanks for the response, Bjorn. I very much appreciate the time. I'm a MATLAB noob, so in working through it, I had a couple questions. The values, 128 for peaks, and kx and ky matrices, why did you choose those values?

The plot for a single image looks like the below, is this because of what you said in your last sentence about adding up contributions to get power?

:

[Bjorn Gustavsson]

Yup, fft2 gives you the Fourier-coefficients for wavenumbers in all four quadrants. This is something you should see as 4 or more points at every unique wave-number amplitude (alopng x-axis). You migh (or not you have to understand what's done in the paper) want to add (the power of) those coefficients toghether. This is just the finding the points on a discrete grid with the same radius, and add those. On the other hand: maybe they use only the vertical and horizontal components, but that seems very strange to me.

[David March]

They state that they "averaged over horizontal and vertical dimensions", so perhaps not adding the coefficients, but taking an average at each amplitude? If so, I'm guessing I can take the average of each column. Though currently X and Y are matrices, and they would need to be only single dimensions. And if I do that, is there a way to make this plot a line like the original?

Granted, I'm not stuck trying to recreate the exact original, just close enough to their process.

[Bjorn Gustavsson]

If the only take the power for purely horizontal and vertical wave-vectors, then that seems utterly stupid to me, unless they have some very good justification for that. The reason this rubs me the wrong way is that it is not rotation-invariant. If you have one dominant harmonic wave-pattern in one direction, then if the camera is oriented such that the rows or columns align with that direction we'd get a large contribution to the power at the corresponding wave-number. If on the other hand the camera is rotated random degrees off that direction the contribution will be zero or very small. I very much hope they have a good enough justification to ignore all other contributions, otherwise this is upsetting...

[David March]

Well I think it's their motivation. This comes from psychology research, and the goal is often to ensure that sets of images are equal in some ways. One of those ways is power, and often our concerns are limitedto things like total or averages, easily digestable to a non-vision scientist or photographer. In this sense, I need to get an average of the power at each frequency for each image in each set (say a set of 10 positive images, a set of 10 negative images). Then I can say, on average, these sets have the same power at high and low frequencies. So I think they are less concerned with direction or quadrant, and more want to capture the entire image.

[Bjorn Gustavsson]

That is an upsetting motivation. It might be marginally tolerable, in the event that they look at machine-generated images. That however, puts some doubt on the validity of the research since biological creatures typically have visual systems that to some degree has to operate with "quite some orientation invariance agility". The reason this is poor analysis is most easily shown with a simple example.

left panels, 2 images with simple harmonic patterns, upper row aligned with rows and columns, middle panel log of the magnitude of the 2-D fft of these images, the upper middle panel have practically 4 delta-spikes (period perfectly matching the image-size so no spectral leakage), the lower panel with the peaks in the diagonal directions from the DC component, spectral leakage due to the implicit periodicity assumed when applying the FFT. right panels are the horizontal and vertical fft-components. Notice how much lower the amplitudes are at the peaks in the lower panel.

Can you share the reference?

[David March]

That's a great example, and one non-vision scientists could likely benefit from being aware of. I've attached a pdf of the article.

Also I want to say that I very much appreciate your concerns with their method. That said, is a way to adapt what you have provided to do something similar, just to appease those who are expecting something similar?

[Bjorn Gustavsson]

Haven't read the paper, but for consumption of the supplemental figure 1, you should have a look at slide 18-19 in this PDF:

Some Image processing fundamentals. There you'll see that the Fourier-amplitudes decrease very similarly for most natural images, and what happens when one swaps the phase of one 2-D fft to another - the images almost swap (admittedly not very well).

[David March]

That's a cool reference, and a very helpful piece of information for interpreting the supplemental figure 1. I'm sure there are indeed issues with that process, as you pointed out earlier. That being said, can I replicate it to a degree with my image sets?

[Bjorn Gustavsson]

Sure, this is what I did to make my illustration:

x = 1:128;
[x,y] = meshgrid(x);
kx = [-64:63];
[kx,ky] = meshgrid(kx);
I1 = sin(x/128*2*pi*16) + sin(y/128*2*pi*8);
I2 = sin((x-y)/2^.5/128*2*pi*16) + sin((y+x)/2^.5/128*2*pi*8);
fI2 = fft2(I2);
fI1 = fft2(I1);
asfI2 =  (abs(fftshift(fI2)));
asfI1 =  (abs(fftshift(fI1)));
subplot(2,3,1)
imagesc(I1)
subplot(2,3,4)
imagesc(I2)
subplot(2,3,2)
pcolor(kx-1/2,ky-1/2,log10(asfI1)),shading flat,set(gca,'TickDir','out')
subplot(2,3,5)
pcolor(kx-1/2,ky-1/2,log10(asfI2)),shading flat,set(gca,'TickDir','out')
subplot(2,3,3)
plot(kx(1,65:end),asfI1(65,65:end),'.-')
hold on 
plot(ky(65:end,1),asfI1(65:end,65),'.-')
set(gca,'YScale','log','XScale','log')
subplot(2,3,6)
plot(kx(1,65:end),asfI2(65,65:end),'.-')
hold on 
plot(ky(65:end,1),asfI2(65:end,65),'.-')

You can just add the amplitudes (or amplitudes squared) of the FFTs from the images together and plot the horizontal and vertical components - the DC ones. If you don't do the fftshift operation you'll find them on the first row and column of the fft:

fI = fft2(Im);
A_hor = abs(fI(1,1:end/2));
A_ver = abs(fI(1:end/2,1));

[David March]

I'll give that a shot. Given I don't understand much of that code, I'll likely come back with another question, or three :)