Extracting Spatial frequency from fourier transform (fft2 on Images)

Steps:

1. fft2 on the Image 2. Extracting Spatial frequency (in Pixels/degree) 3. Plotting magnitude of the fourier transform (power spectral density of the image) Vs Spatial frequency

Now, for step 2 - Rick Rosson's answer on this question might give some directions. http://www.mathworks.com/matlabcentral/answers/13896

But, for step 3 - I have to be able to relate the magnitude of the fourier transform with the corresponding spatial frequencies, to plot them. I am not clear how to extract this correlation from the MATLAB fft2 output. Any help would be appreciated.

PS: For people interested in the Optics part of it, I am trying to plot the Modulation Transfer Function (Figure 7 in the link below). The method is outlined here: http://maxwell.uncc.edu/gboreman/publications/Daniels_OE_4_1995.pdf

Answer 1

[Wayne King]

If you have a simple 2-D complex exponential

x = 0:159;
y = 0:159;
[X,Y] = meshgrid(x,y);
Z = exp(1i*pi/4*(X+Y));
Z = real(Z);
imagesc(Z);

Then note that with a matrix of 160x160, we would expect pi/4 to be shifted from the center (80,80) to (101,101) and (61,61) (if you are using fftshift to move (0,0) to the center of the image)

    Zdft = fftshift(fft2(Z));
    imagesc(abs(Zdft))
    Zdft(101,101)
    Zdft(61,61)

Note they are complex conjugates of each other as expected.

On a matrix 160x160, each step will move you (2*pi)/160 radians/(unit length) in angular spatial frequency, so that's why find pi/4 essentially 20 pixels away from the center of the image (0,0).

Answer 2

[Mona Mahboob Kanafi]

If I get it right, you need a 2D graph of power spectral density versus frequencies. When you do a 2D fft on an image, you get a 2D matrix in matlab representing the fourier transform of your image. You then just need to assign fx and fy in order to plot the 3D graph of FFT. Now, to get power spectral density :

PSD = (a^2/n*m)*(abs(FFT))

where a is pixel width and n & m are number of pixels in each direction.

With this, you obtain the PSD. Now you again obtain a 3D graph of PSD versus fx and fy, But what you want is a 2D graph of PSD-f!

Here it depends on your system. If the system is isotropic then I would do a radial averaging over the 2D FFT which gives me a 2D graph of PSD-fr. Or if it isn't, then I wouldn't use ff2. Instead, I would use fft over parallel lines in the direction I want(say x) and averaged them together in order to obtain a 2D PSD versus fx!

[Gokul Raju]

Thanks very much for the response Mona.

Yes, you got what exactly my problem is. I agree that I can use either fft2 or fft to get the 2D plot. I was actually using fft2, but computed the row-average for my plot. That is clear for me!

My problem solely lies with 'assign fx or fy to plot the graph'. I could still not visualise how to correlate the fft-output-matrix with fx or fy, and get the plot done. It would be really helpful to have a small code snippet or pseudocode/algorithm on this.

Lets say, I have a 1280*720 Image and the field of View being 33.5*19 degrees. I have ~38 pixels/degree on both x and y. Now, can you show me what I do with the fft output, to plot the 2D PSD vs spatial frequency ?.

I appreciate your time on this! Gokul

Answer 3

[Mona Mahboob Kanafi]

qx_1=zeros(m,1);

for k=0:m-1

    qx_1(k+1)=(2*pi/m)*(k);

end

qx_2 = fftshift(qx_1);

% for removing discontinuity in the middle, shift data by 2pi

qx_3 = unwrap(qx_2-2*pi);

qx=qx_3/a;

qy_1=zeros(n,1);

for k=0:n-1

    qy_1(k+1)=(2*pi/n)*(k);

end

qy_2 = fftshift(qy_1);

qy_3 = unwrap(qy_2-2*pi);

qy=qy_3/a;

I understood this by reading steve's blog on this frequency assignment, you can also refer to that:

http://blogs.mathworks.com/steve/2010/06/25/plotting-the-dtft-using-the-output-of-fft/

Just one more hint for you, I believe averaging row data on fft2 does not give you correct result. Just take a look at a 2D fft graph that has been centralized by fftshift. Or at least in my case it doesn't give correct results, so I suggest the methods I explained above.

good luck

[Gokul Raju]

Mona, If I understand this right, qx and qy should contain the values of spatial frequencies in their respective axis, for plotting against the PSD values.

I tried the above snippet with m=1280, n=720 and a=4.65um. All I got in the qx and qy matrices are same values, I mean the entire column matrices in each case are filled with identical elements. Is something wrong ?

btw, I replaces mm,nn with m*m,n*n in the code.

[Mona Mahboob Kanafi]

Sorry, nn and mm are n and m, I fixed them above. after fixing this, what you get is two matrices, one of them 1*1280 and the other one 1*720. The maximum frequency in both fx and fy are the same, because maximum frequency is related to the minimum wavelength and that is pixel width(fmin = 1/a), but according to Nyquist frequency, the maximum f can't exceed 1/2a! To sum up, your maximum f are always 1/2a (or 2pi/2a if you want to work with wavevector rather than frequency).

You minimum f is related to maximum available wavelength in you system which is size of your image in each direction, so min fx is different from min fy. i.e. fyMax = 1/a*n & fxMax =1/a*m